Impact on splicing in Saccharomyces cerevisiae of random 50-base sequences inserted into an intron

Generation of a library of introns with random 50-nt elements

We sought to analyze splicing in S. cerevisiae by assessing the effect of a large number of random sequence elements inserted into a common intron. We constructed a library of over 750,000 introns, each containing a unique N50 element. The intron background for the library was based on the native S. cerevisiae ACT1 intron sequence of 309 nt. We shortened the ACT1 intron to 115 nt by removing the centermost nucleotides between the 5′ splice site and branchpoint sequence; this shorter ACT1 intron served as a size control (Fig. 1A). For the library, 50 nt in the region between the 5′ splice site and the branchpoint, the majority of the intron sequence in this region, were replaced by synthetic oligos of 50 random bases (Fig. 1A). However, the 5 nt just downstream from the 5′ splice site sequence and the 5 nt just upstream of the branchpoint sequence were maintained to avoid effects on splicing efficiency due to changes that can influence the direct binding of the spliceosome to these consensus sites (McGrail et al. 2006). The N50 oligos had a generally uniform distribution of all 4 nt (Supplemental Fig. S1A).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 1.

Experimental design and analysis of splicing in controls. (A) Design of constructs. The background for the library was created by removing the centermost nucleotides between the 5′ splice site and branchpoint sequences of the native S. cerevisiae ACT1 intron. A random 50-nt (N50) sequence replaces bases 12 to 61 of the intron between the 5′ splice site and the branchpoint of the shortened ACT1 intron. The intron variants were inserted into the GFP gene 10 nt downstream from the start codon with a unique barcode added 38 nt upstream of the GFP coding sequence, which was under control of the CYC1 promoter and terminator. Numbers indicate the location relative to the GFP open reading frame start site. The different sequence contexts used for analysis are the “nearby” context spanning from 15 bases upstream to 15 bases downstream from the N50 (80 bases total), the “barcode” context spanning from the start of the barcode through 15 bases downstream from the N50 (124 bases total), the “upstream” context spanning from the transcriptional start site (TSS) through 15 bases downstream from the N50 (174 bases total), and a “broad” context spanning from the transcriptional start site through 15 bases downstream from the intron (259 bases total). (B) Analysis of splicing of the full-length and shortened ACT1 introns by semiquantitative RT-PCR using RNA from exponentially growing liquid yeast cultures. PCR using DNA (left) and cDNA (right), unspliced (upper band) and spliced (lower band). (C) Expression analysis determined by GFP intensity measured via flow cytometry using cells containing plasmids with GFP and either no intron, the full-length, or the shortened ACT1 intron and cells without GFP. (D) Splicing of the shortened ACT1 intron paired with 1905 unique barcode variants as assayed by sequencing. (E) Experimental workflow. The N50 intron library was transformed into yeast and RNA was isolated. The splicing efficiency, or the ratio of spliced versus total transcripts, was measured using targeted sequencing of the intron, and variants were identified by their unique barcode.

Each intron was placed in a green fluorescent protein (GFP) reporter gene at a position 10 nt downstream from the start codon, as in the native ACT1 gene (Fig. 1A). Transcription of the reporter gene was under the control of the S. cerevisiae CYC1 promoter and terminator (Fig. 1A). Each intron variant was paired with a unique 15-nt barcode in the 5′ untranslated region (UTR), 38 nt upstream of the start codon, such that the barcode remains in the mRNA after splicing to identify the pre-mRNA N50 variant for each transcript (Fig. 1A). The constructs were sequenced to match barcodes to N50 elements. The barcodes contained an approximately even distribution of the 3 nt other than thymine, which was excluded to avoid the creation of an upstream start site (Supplemental Fig. S1A). Four sequence contexts were used for analysis, including a “nearby” context spanning 15 bases each side of the N50 (80 bases total); a “barcode” context, spanning from the start of the barcode through 15 bases 3′ of the N50 (124 bases); an “upstream” context spanning from the transcriptional start site through 15 bases downstream from the N50 (174 bases); and a “broad” context spanning from the transcriptional start site through 15 bases downstream from the intron (259 bases) (Fig. 1A).

To establish the splicing efficiency—defined as the fraction of reporter gene transcripts that were spliced correctly—initially of control constructs, we transformed yeast with plasmids containing either the full-length or the shortened ACT1 intron within the GFP reporter gene. We analyzed transcript abundance using semiquantitative reverse transcription PCR. For the full-length ACT1 intron construct, only the spliced mature transcript was visible by gel electrophoresis (Fig. 1B), comparable to the complete splicing of this intron in other contexts (Vijayraghavan et al. 1986; Agarwal and Ansari 2016). In the shortened intron construct, ∼82% of the transcripts were visible as the spliced form (Fig. 1B). We also assessed protein expression from the GFP reporters using flow cytometry to measure the fluorescence of exponentially growing yeast. The full-length or shortened intron-containing GFP constructs displayed higher average GFP fluorescence than the “no intron” control (Fig. 1C), in line with yeast studies that show that genes with an intron have higher mRNA stability and transcription (Furger et al. 2002; Parenteau et al. 2008; Moabbi et al. 2012). Yeast with the shortened intron displayed 66% of the fluorescence relative to the full-length intron (Fig. 1C), establishing that GFP protein production and RNA splicing efficiency measured by gel electrophoresis were largely concordant.

To determine whether the upstream barcode affected splicing, we paired the shortened ACT1 control with ∼2000 barcodes and measured splicing efficiency by targeted cDNA synthesis, amplification and sequencing across the 3′ splice junction. The shortened intron variant spliced with 97.0% efficiency in this assessment. The average splicing efficiency of all barcode variants of the shortened intron was 96.2%, with only 0.7% variance between different barcodes (Fig. 1D). Although this efficiency suggests that the barcodes had minimal impact on splicing, 4% of the barcoded controls displayed a splicing efficiency of 80% or less (Fig. 1D). Features of the barcodes may affect splicing, although not their GC content, as splicing was similar among transcripts carrying barcodes of different GC content (Supplemental Fig. S1B). The analyses of RNA and protein expression demonstrate that the reporter construct is a high-splicing background; however, as it is not spliced as efficiently as the full-length ACT1 intron, it can be used to assess insertions that cause either positive or negative effects on splicing (Fig. 1E).

Intron variants with random N50 elements display a wide range of splicing efficiencies

We initially characterized the overall splicing efficiency of the intron library. DNA and RNA were extracted from two independent cultures of yeast cells containing the library. To distinguish each transcript, we added unique molecular identifiers (UMIs) to the gene-specific primer used during reverse transcription. Using massively parallel sequencing, we identified each N50 element by its corresponding barcode and determined the abundance of different splicing isoforms by aligning sequencing reads across the 3′ splice junction to either the intron (i.e., unspliced) or the first exon (i.e., spliced). The UMIs were deduplicated to remove PCR copies, and unique transcripts were counted as the number of UMIs per variant. To obtain higher confidence data, we included in further analyses only variants detected with at least three RNA transcripts and 20 DNA reads (152,514 variants in Rep 1; 193,200 in Rep 2; 141,710 overlap between both replicates; Supplemental Table S1). The splicing efficiency of N50 variants in the two biological replicates was well correlated (Pearson's R = 0.96; Fig. 2A). Using averages of the two replicates in the overlap set, we classified variants into three groups: high (>80%), moderate (20% to 80%), or low (<20%) splicing efficiency (Fig. 2A). The library has a large fraction of variants that showed reduced splicing efficiency compared to the control (Figs. 1D, 2A,B). Variants displayed the full range of splicing efficiencies from 0% to 100% (Fig. 2A,B), with 72% of the variants splicing at lower than 80% (Fig. 2B). This reduced splicing efficiency indicates that features of the N50 elements have a substantial impact on the splicing process.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 2.

Sequencing analysis of splicing efficiency and transcript abundance in the N50 library. (A) Splicing efficiency of 141,710 intron variants present in two independent biological replicates of the intron library. Variants that displayed low (under 20%), moderate (20% to 80%), or high (above 80%) splicing efficiency in both replicates are boxed. (B–D) Averages of the two replicates are shown without overabundant outliers (>2000 reads/transcript, n = 141,702). (B) Histogram of RNA splicing efficiency displaying the distribution and range of splicing efficiencies of the N50 library. (C,D) RNA abundance was determined by the log₁₀ of the total number of transcripts relative to total DNA for each variant. (C) RNA abundance versus splicing efficiency of each variant averaged across each splicing group (i.e., low, moderate or high). Brackets [ ] indicate inclusive values and parentheses ( ) indicate exclusive values. Horizontal brackets on top indicate significant differences between compared positions using a two-sample t-test, with three asterisks denoting P-values <5⁻¹⁶. The gold diamonds indicate the mean for each group. (D) RNA abundance versus splicing efficiency of each variant individually. Pearson's correlation coefficients are denoted by R.

As splicing efficiency and steady-state RNA levels are correlated (Clancy and Hannah 2002; Ding and Elowitz 2019; Schirman et al. 2021), we examined the abundance of transcripts containing the N50 elements. To account for variation in the frequencies of the N50 elements in the population, we normalized each variant's total RNA to its DNA and termed this value “RNA abundance,” though overall, RNA and DNA levels moderately correlated (R = 0.57, Supplemental Fig. S1C). Intron variants that displayed moderate splicing efficiency had on average 33% lower RNA abundance, and those with low splicing efficiency a striking 65% lower, than those with high efficiency (Fig. 2C). Furthermore, unspliced transcripts were 50% less abundant than spliced transcripts (UMIs per variant, Supplemental Table S1). Transcripts with incorrect splicing can be degraded by nonsense-mediated decay (for review, see He and Jacobson 2015), thereby reducing the accumulation of unspliced transcripts and protecting cells against potentially deleterious proteins produced from these transcripts. Overall, the intron library displayed a positive correlation between splicing efficiency and RNA abundance (R = 0.38, Fig. 2D), supporting the idea that unspliced transcript isoforms have reduced accumulation.

Hexamers enriched in GU reduce splicing efficiency

To identify short sequences in the N50 elements that may contribute to the broad range of splicing efficiencies observed in the intron library, we calculated the average splicing efficiency of library members containing every possible 6-nt-long sequence. We ranked each hexamer from highest (rank 1) to lowest (rank 4096) splicing efficiency, based on the average splicing of introns containing them, and inspected hexamers more than 2.5 standard deviations away from the mean (Fig. 3A; Supplemental Table S2). The hexamers with the highest associated splicing efficiency were AU-rich (Supplemental Table S2). Similarly, hexamers associated with introns with the lowest splicing efficiency rank were generally GC-rich. However, 20 of the 25 hexamers with the lowest associated average splicing efficiency contained at least one GU or UG dinucleotide, and 11 contained multiple repeats, with GUGUGU ranked lowest of all hexamers (Fig. 3A; Supplemental Table S2). The effect of this GU-rich element on splicing efficiency was not dependent on its specific location within the N50 (Fig. 3B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 3.

Analysis of all possible hexamers within the N50 library variants. (A–D) Values averaged between two biological replicates, n = 141,710. (A) Average splicing efficiencies of N50s containing each hexamer ranked highest (#1) to lowest (#4096). (B) Average splicing efficiencies of N50s containing GUGUGU and starting at each possible position within the N50. (C) Average RNA abundance of N50s containing each hexamer ranked highest (#1) to lowest (#4096). (D) Average minimum free energy (MFE) of the main predicted structure for the library variants using the N50 sequence and its “upstream” context (entire upstream RNA sequence and 15 bases downstream, 174 bases total) versus the “nearby” context (15 bases each side of N50, 80 bases total). The red numbered dots indicate the specific hexamers listed.

In total, the 25 hexamers with the largest reductions on splicing efficiency were found in ∼26% of the N50s, due in part to a slight bias for both guanine and thymine in library synthesis (Supplemental Fig. S1A). We compared the average splicing efficiencies of N50 elements that contain either contiguous GU or UG repeats (e.g., GUGU) to those with noncontiguous repeats (e.g., GUAAGU). In all cases, N50s with contiguous dinucleotide repeats had on average lower splicing efficiency than those with noncontiguous dinucleotide repeats (Supplemental Table S2). Although GUGUGU has not been reported to have a drastic effect on splicing when present in native yeast introns, we analyzed the abundance of short sequences in native introns, which might reflect evolutionary selection on sequence content (due to the limited number of native introns, we analyzed all 1024 possible pentamers rather than hexamers). GUGUG is enriched over 896 other pentamers and UGUGU over 919 other pentamers, indicating that GU repeats are not preferentially excluded from S. cerevisiae introns (Supplemental Table S3). Conversely, CG-rich pentamers tend to have lower enrichment over a shuffled pentamer (Supplemental Table S3). Taken altogether, these results suggest that GC content has a general negative effect on splicing, whereas the effect of the GU repeats is likely due to the specific context of the reporter gene.

We also determined average RNA abundances across N50 elements containing each hexamer and again ranked each hexamer from highest to lowest. N50s containing hexamers with GU repeats were not among those ranked lowest in average RNA abundance, with 29 and 84 sequences ranked lower than GUGUGU and UGUGUG, respectively (Supplemental Table S2). Rather, the lowest-ranked hexamers for RNA abundance were GC-rich, with CCGCCG corresponding to the least abundant RNA on average (Fig. 3C; Supplemental Table S2).

Secondary structure and intron GC content affect splicing efficiency

We considered whether the GU-rich N50 elements might be forming some secondary structure within the reporter gene pre-mRNA that inhibits splicing. Because of wobble-pairing of guanine, GU-rich sequences could potentially hybridize to a large range of complementary sequences within the pre-mRNA. However, a single AC-rich sequence, ACACACAC, is present in the 5′ UTR of the reporter gene, beginning 11 bases upstream of the barcode. To determine if this sequence may be related to the reduced splicing efficiency of introns with GU-rich N50 elements, we calculated the lowest MFE structure of regions containing the N50 element and variable additional sequence (Supplemental Fig. S2; Supplemental Table S4). The correlation of MFE to splicing efficiency was greatest (R = 0.23) using an “upstream” context (Fig. 1A) that contains the AC repeat region (Supplemental Fig. S2; Supplemental Table S4). Exclusion of the full 5′ UTR sequence using “nearby” or “barcode” contexts (Fig. 1A) resulted in much lower correlations between MFE and splicing efficiency (Supplemental Fig. S2; Supplemental Table S4). Additionally, some hexamers had a much stronger effect on predicted secondary structure when 5′ UTR sequences were included than when they were not (Fig. 3D; Supplemental Table S2). These outliers included GUGUGU and similar GU-rich sequences, as observed in the lowest-ranked hexamers for associated splicing efficiency (Fig. 3D; Supplemental Table S2). These results support the idea that the GU-rich sequences in the N50 elements form a secondary structure with the 5′ UTR that reduces splicing, suggesting that context is an important consideration in splicing.

We also examined energy of the MFE structure and its frequency in the ensemble of all structures, the free energy and the diversity of the thermodynamic ensemble of structures and the base-pairing probability between each pair of nucleotides, across the N50. The MFE of the main structure (R = 0.21) and of the ensemble of structures (R = 0.22) were most strongly correlated with splicing efficiency among the predicted secondary structure features using the full N50 sequences (P-value < 2.2 × 10⁻¹⁶, Supplemental Table S4). The lower the MFE of the ensemble, the lower the splicing efficiency, suggesting that formation of strong RNA structures reduces splicing efficiency. This relationship is congruent with results that show that specific secondary structures of the pre-mRNA near or involving the three consensus intron sequences can reduce splicing efficiency and alter splice site selection by hindering spliceosome access to these sites (Eperon et al. 1988; Halfter and Gallwitz 1988; Goguel and Rosbash 1993; Mougin et al. 1996; Singh et al. 2007; Barrass et al. 2015; Zafrir and Tuller 2015). Across the N50, the base-pair probability increased in variants with moderate versus high splicing efficiency and was greatest in the low splicing group (Supplemental Fig. S3A). The base-pair probability was also consistently higher across the AC-repeat region of the 5′ UTR for variants with low splicing relative to the other groups (Supplemental Fig. S3B). Congruently, variants that contain the GUGUGU hexamer displayed a much higher base-pair probability in this region than any of the other splicing groups (Supplemental Fig. S3B) as well as the other transcript regions (Supplemental Fig. S3A–C). This correlation further supports the idea that secondary structure between the intron and 5′ UTR can strongly impact splicing. For the 5′ splice sequence, there was no consistent trend for base-pair probability for the three splicing groups, and the site-specific differences across the 5′ splice site were greater than between any of the splicing groups, indicating minimal interaction of the 5′ splice site with the N50 for splicing (Supplemental Fig. S3C). The other features of the predicted secondary structures (the frequency of the main structure in the ensemble of all predicted structures and the ensemble diversity) had little to no correlation with the splicing efficiency of the variants (R < |0.02|, Supplemental Table S4).

We also analyzed GC content across the entire N50 as well as smaller windows of 10 nt, and the predicted RNA secondary structure across a “broad” context (259 bases) of the transcript spanning from the transcriptional start site through 15 bases downstream from the intron (Fig. 1A). Generally, introns with higher GC content displayed lower splicing efficiency, with a correlation of R = −0.07 averaged across all 10-nt windows and R = −0.14 for GC across N50s (Supplemental Fig. S3D,E), similar to prior studies (Galante et al. 2004; Wong et al. 2013; Schirman et al. 2021; Gnan et al. 2022). These correlations are in line with the low splicing efficiency rank of hexamers with high GC content (Supplemental Table S2). The differences in both measures of GC content among the low, moderate and high splicing efficiency groups were greatest for regions at the edges of the N50 elements relative to the middle (Supplemental Fig. S3D).

Motifs associated with reduced splicing efficiency

We analyzed the N50 elements of RNAs that displayed high, moderate, or low splicing efficiency in both biological replicates (Fig. 2A) for sequence motifs that were significantly enriched. We identified three motifs enriched in the low splicing group, five in the moderate splicing group and one in the high splicing group, when variants were grouped by unique barcodes (Fig. 4A). One of the motifs from the low splicing group and four from the moderate splicing group were similarly found when variants were grouped by unique barcode and N50 (Supplemental Fig. S4A). To validate these results, we chose representative motifs from each group: (i) UUAGUG, (ii) AGGCUUCGGG and (iii) GAUGUUUGAAUA from the low-splicing group; (iv) UUUGAAUAAUGAC, (v) AAAUGUUCGA, (vi) UUUAGAUAUAUGUC, (vii) GACUUUAGAU and (viii) ACUUAUUU from the moderate-splicing group; and (ix) CCGUGUAG from the high-splicing group (Supplemental Fig. S4A). Validation controls of a silencer motif (UUUGUGUA, designated Y1) identified from native yeast introns and a mutated motif that acts as a silencer (UUUAUGCU, designated Y2) as well as a splicing enhancer (GUACAUGU, designated Y3) were chosen from Yofe et al. (2014) (Fig. 4A).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 4.

Overview and sequencing analysis of motifs enriched from splicing efficiency status in different intron backgrounds. (A–F) Results for the average of the two replicates. (A) Motifs that were enriched among the variants with unique barcodes that displayed low (under 20%), moderate (20% to 80%), or high (above 80%) splicing efficiency in both replicates relative to a randomization of those sequences with P-values <0.05. Also listed are motifs previously reported to reduce (Y1 and Y2) or enhance (Y3) splicing in Yofe et al. (2014). (B) Overview of the different sequence backgrounds and positions of motifs tested. Numbers indicate the position in the different 50-nt test regions, and regions of the introns are denoted (ss) splice site, (bp) branchpoint. (C) Average splicing efficiency for all motif variants and positions as well as the no motif controls, denoted by a gold rectangle, for each background sequence condition, n = 336. (D) The difference in splicing efficiency for each motif variant from the no motif control with the corresponding background sequence, averaged across all positions and backgrounds, n = 28. (E) The difference in splicing efficiency for each motif position from no motif control with the corresponding background sequence, averaged across the five low splicing motifs and all backgrounds, n = 35. Brackets indicate significant differences between compared positions using a two-sample t-test, with one asterisk denoting P-values <0.05 and three asterisks denoting a P-value <0.0005. (F) The difference in splicing efficiency for motifs in each background from no motif control with the corresponding background sequence, averaged across the five low splicing motifs and all positions, n = 20. Single asterisks on the top of the graphs indicate significant differences from the no motif control within each group when using a one-sample t-test with P-values <0.05.

The enriched motifs were introduced into seven backgrounds and at three locations within the N50 to test their efficacy in a broad range of settings (Fig. 4B). The seven backgrounds were: one N50 from each of the three splicing efficiency groups that did not contain any of the motifs; the shortened ACT1 intron sequence with a medium GC content (40%), as well as low (20%) and high (60%) GC content versions constructed by replacing bases as needed; and the shortened ACT1 intron paired with a mutation of the 5′ splice site that substantially reduces splicing (Fig. 4B; Vijayraghavan et al. 1986). To assess positional dependency, we placed the motifs at the start of the N50 (the initial nucleotide of the motif placed at position 1 of the N50), the middle (the central nucleotide of the motif placed at position 25) and the end (the last nucleotide of the motif placed at position 50) (Fig. 4B). We also assessed constructs that had the same motif at all three positions as well as each background containing no motif (Fig. 4B). In all, we tested 12 motifs in four positional configurations and seven backgrounds, as well as seven no motif control backgrounds, for a total of 343 unique introns.

To test whether the enriched motifs influenced splicing efficiency independent of context, we first compared the splicing efficiency of the seven control backgrounds (no motif). The backgrounds chosen from the low, moderate and high splicing groups displayed the same overall trends in the validation experiment as in the library experiment (Fig. 4C). For backgrounds with higher GC content, there was a corresponding decrease in splicing efficiency (Fig. 4C), again as in the library experiment (Supplemental Fig. S3D,E). The control background with the 5′ splice site mutation had a splicing efficiency of 52% (Fig. 4C). To compare the effect of the motifs in these different backgrounds, we calculated the difference in splicing efficiency of each intron compared to the no motif control for each background. A significant decrease in splicing was observed for the three motifs (1, 2, and 3) from the low splicing group, as well as for the two previously reported splicing-reducing motifs Y1 and Y2 (Fig. 4D). No differences in splicing efficiency were detected for motifs 4, 6, 7, or 8 from the moderate splicing group or from the previously reported splicing-enhancing motif Y3 (Fig. 4D). Motif 5, from the moderate splicing group, and motif 9, from the high-splicing group, were highly variable, but showed a significant reduction in splicing across backgrounds (Fig. 4D). These data suggest that only motifs that strongly reduce splicing maintain their effect independent of the background.

We assessed the position and background specific effects for all low splicing motifs. The average splicing for all low splicing motifs, including the three from the low splicing group and the two from Yofe et al. (2014), were observed in all positions and were stronger when the motifs were present in triplicate (Fig. 4E). In analyzing background-specific effects, we found that these five low splicing motifs displayed lower average splicing efficiency in all backgrounds except the low and high splicing backgrounds (Fig. 4F). These results indicate that the motifs were not as effective when present in a background that had a strong enhancing or reducing impact on splicing. In contrast, the motifs further reduced splicing efficiency in the poor splicing background created by the 5′ splice site mutation (Fig. 4F), suggesting that the 5′ splice site and more 3′ regions of the intron impact splicing in different ways. Overall, the results demonstrate that multiple motifs can affect splicing efficiency in context-dependent and context-independent manner.

N50 variants and cryptic splicing

Cryptic splice sites are noncanonical transcript locations that undergo splicing aberrantly, thereby competing with and reducing the amount of splicing at the canonical site (Nelson and Green 1990; Cunningham et al. 1991; Lesser and Guthrie 1993; Rosenberg et al. 2015). To examine 5′ and 3′ cryptic splice sites in the variants, we reexamined sequences that did not match the expected junction sequences of canonically spliced or unspliced variants. We identified 5′ cryptically spliced variants if their junction sequence exactly matched a sequence anywhere downstream from the barcode, upstream of the branchpoint, and not at the canonical site (Fig. 5A). Cryptic 3′ splice sites were identified by junction sequences that contained combinations of the last six bases of exon 1, assuming these transcripts are spliced at the canonical 5′ splice site, and at least three bases of 3′ unspliced sequence (Fig. 5A). Since we sequenced the intron–exon junction using a reverse primer that anneals to the start of exon 2, no downstream 3′ cryptic splice sites could be identified (Fig. 5A). Approximately 1% of the N50 variants (1017 variants) had at least one transcript that was cryptically spliced. Most of the 5′ cryptic splices fall within the first 25 nt of the N50 elements (Fig. 5B), likely because there is a minimum length required upstream of the branchpoint to form a lariat (Wieringa et al. 1984; Smith and Nadal-Ginard 1989) and 5′ cryptic splices in the last 25 nt of the N50 would not meet this minimum.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 5.

Overview and sequencing analysis of cryptically spliced N50 variants. (A) Diagram of the different types of splicing in library variants. Cryptic 5′ and 3′ splicing are depicted by a red dashed line and a blue dash-dot line, respectively. Splicing at a cryptic 5′ or 3′ splice site (ss) results in the inclusion of the 5′ or 3′ end of the intron, correspondingly. The sequencing primer to identify the differentially spliced variants anneals to the junction between the intron and the second exon. (B,C) Results for the average of both replicates. (B) The number of 5′ and 3′ cryptically spliced variants at each position of the N50. (C) The splicing efficiency versus the portion of cryptically spliced transcripts in 5′ and 3′ cryptically spliced variants. (D) Portion of cryptic splicing per variant found in both independent biological replicates of the intron library (n = 237). Pearson's correlation coefficient is denoted by R. (E) The significantly enriched sequence motif among the 5′ cryptically spliced variants is listed, P-value <0.05. (F) The consensus motif of the actual cryptic 5′ splice sites identified in the library. (E,G) The number of variants containing the enriched motif at each starting position found in 5′ cryptically spliced variants with (E) only 5′ cryptically spliced variants or (G) the entire library. (H) The four most common sites used for 3′ cryptic splicing among the intron variants.

Cryptic 3′ splice sites were similarly constrained, found only within 30 nt of the canonical 3′ splice site and only at specific locations (Fig. 5B). In the control shortened ACT1 intron sequence, no 5′ cryptic splicing was observed (Supplemental Fig. S4B). However, six of the barcode variants displayed 3′ cryptic splicing at the same top two locations as the 3′ cryptically spliced library variants, indicating that only a few positions present a viable splice sequence (Fig. 5B; Supplemental Fig. S4B). No relationship was found between canonical splicing efficiency and the proportion of transcripts cryptically spliced (Fig. 5C). For 25% of the N50 elements associated with cryptic splicing, the cryptically spliced transcripts accounted for the majority of transcripts (Fig. 5C). Only 23% of the 5′ and none of the 3′ cryptically spliced variants were found in both biological replicates. However, the proportion of cryptic splicing for these 23% (237 variants) was well correlated (R = 0.78; Fig. 5D). The low abundance and lack of replicable 3′ cryptic splicing compared to 5′ cryptic splicing may indicate that the N50 region had comparatively little impact on 3′ splice site selection and is likely not a general feature of cryptic splicing. For N50 elements in which 5′ cryptic splice sites were identified, we searched for enriched sequence motifs. The only motif identified contains the canonical 5′ splice site, GUAUGU, found in 73% of the 5′ cryptically spliced variants (Fig. 5E); this site provided the primary 5′ splice site used in the variants with 5′ cryptic splicing (Fig. 5F). The enriched motif is primarily located in the first half of the 5′ cryptically spliced variants (Fig. 5E), but evenly distributed in the canonically spliced library variants (Fig. 5G). For the 3′ cryptic splice sites, we found that the three bases most commonly used as the cryptic splice site did not match the canonical 3′ splice sequence (5′-[C/U]AG-3′, Fig. 5H). It is possible that cryptic splicing occurs when another site is present that can be bound and cleaved by the spliceosome, even if it is weaker than the canonical splice site.

Predicting splicing efficiency from intron sequence using multiple models

To determine sequence patterns in the N50 element that influence splicing efficiency, we trained five convolutional neural network (CNN) models with different architectures. For 140,017 sequence variants (the set of overlapped variants with unique barcodes), we trained on 80% of the data, validated on another 10% and tested performance on the remaining 10%. While these models made predictions that are significantly different from random, none were able to explain the majority of the variance in splicing efficiency (max Pearson's R = 0.26, Supplemental Table S5). Down-sampling of the number of sequences in the training set suggested that the CNN models require many more training sequences to learn the patterns that are controlling splicing efficiency.

We also applied linear regression, random forest regression and elastic-net models to relate splicing efficiency to 14 known features of each N50 alone or to these features along with k-mers in the N50 (see Materials and Methods for feature details; Table 1; Supplemental Table S6). The simple linear regression outperformed the random forest and elastic-net models. The best linear regression model, which included the 14 features as well as the presence of trimers in the N50, performed almost on par with the more complex CNN models, R = 0.24 versus R = 0.26, respectively (Table 1; Supplemental Table S5). The most important feature in determining splicing efficiency in this linear model was the predicted MFE of the structural ensemble, followed by the GC-content of the N50, the presence of an additional 5′ splice site, the presence of hexamers ranked highest in splicing efficiency, and the presence of hexamers ranked lowest in splicing efficiency (Table 1). Due to its inherent sparsity, the elastic-net model trained on sequence k-mers within the N50 was also used to determine which sequence patterns have the strongest influence on splicing efficiency. The two most predictive features were the occurrence of GUG and GCG elements, predicted to negatively affect splicing efficiency (Table 1). However, structure predictions did not indicate the formation of specific base-pairing between the N50 and other regions of the pre-mRNA, suggesting either that the predictions lack accuracy or that the formation of secondary structure is insufficient to explain the splicing efficiency.