The almost invariant consensus sequence for mRNA splicing in animals and plants is gu_ag, where gu is the splice donor sequence and ag is the splice acceptor sequence. A longer splice donor consensus sequence in most mammals is guragu, where r is either g or a (Mount, 1982; Black, 2003). Usually, an expression of “gu_ag” means that only the 5′ and 3′ terminal two nucleotides of the sequence are invariable as gu and ag, respectively, and that a sequence represented by the underscore can be any sequences. However, here we use this expression to indicate that the sequence represented by the underscore can be any sequences except for sequences that do not match any of the other consensus sequences.

The splice acceptor consensus sequence is preceded by a branch point sequence, which contains an adenine, which is ligated to the 5′ splice site ribonucleotide to form the intron lariat, and a polypyrimidine tract (c or u), which is between the branch point and the splice acceptor sequence. While the short gu_ag consensus sequence of introns is clearly not sufficient to differentiate amongst the multitude of alternative splicing events, surprisingly little is known about what other sequence information is required to regulate alternative RNA splicing (Ladd and Cooper, 2002; Barash et al., 2010; Witten and Ule, 2011). The flanking one or two nucleotides on either side of the intron are also often conserved, and they are included in our supplementary tables, but they will not be discussed further in this paper so that we can focus our analyses on consensus sequences at the ends of the introns.

Alternative RNA splicing occurs in almost all human genes and vastly increases the number of proteins and transcripts that an organism can produce (Pan et al., 2008; Wang et al., 2008). Exons that are involved in all types of RNA splicing can be classified into five major categories: (1) exons containing alternative 5′ splice sites (A5), (2) exons containing alternative 3′ splice sites (A3), (3) retained or invariant exons (R), (4) skipped exons (S), and (5) mutually exclusive exons (ME) (Ast, 2004; Sugnet et al., 2004). In addition to these five types of exons, we also include in our bioinformatics analysis exons that contain an alternative promoter (APr) and exons that contain an alternative poly A (APA) site to make a total of 7 exon types. Since an intron is flanked by two exons, APr can only be at the 5′ end, and APA can only be at the 3′ end, there are 36 possible pair-wise types of introns that are distinguished by the combinations of 7 types of flanking exons. Here we present bioinformatics evidence to support our hypothesis that there is an alternative splicing code based on the 36 different types of introns.

We hypothesize that there are 36 different types of introns with unique consensus sequences based on their flanking exon types. To test this hypothesis, we analyzed the 36 different types of introns separately rather than combining all introns into one pool, as is usually done in bioinformatics studies of introns (e.g., Mount, 1982). It is possible, for instance, that there are many different types of macromolecular splicesosome complexes, in addition to the canonical U2-type splicesosomes and the non-canonical U12-type splicesosomes (Padgett, 2012), that utilize the numerous U2-variant small RNA sequences in the genome (O'Reilly et al., 2013). To test this hypothesis, we determined whether the 36 types of introns are enriched for a particular paired consensus sequence(s) that is derived from both ends of the intron and flanking exon regions. Supplementary Table 1 has a list of all 36 types of introns and the number of introns in each class in H. sapiens, M. musculis, D. melanogaster, C. elegans, and A. thaliana.

To begin our bioinformatics analysis of introns, we first generated a table of paired splice donor and acceptor consensus sequences, from the most common to the least common (see Materials and Methods). For statistical reasons, we selected an arbitrary cutoff of each paired consensus sequence being represented by at least 100 introns. Using a modification of our program SnpEff, which classifies sequences in any sequenced genome (Cingolani et al., 2012), we analyzed the genomes and compared them with annotated full-length RefSeq cDNAs of 50 different plant and animal species. The total number of different types of paired consensus sequences ranged from one in baker's yeast, Saccharomyces cerevisiae, guaugu_ag, which has only 282 introns, all of which are always flanked by invariant exons (R-R), to 95 different consensus sequences in the marmoset, which has 184,882 introns (Supplementary Table 2). The average number of introns in the 50 species that we analyzed was 116,288 with a standard deviation of 45,266. Almost half of the animals' genomes that we analyzed have between 40 and 50 different types of paired consensus sequences with at least 100 introns in each type. Table 1 shows the 42 consensus sequences in humans, which are in at least 100 introns, in rank order from most common to least common.

Second, the 42 different paired consensus sequences in humans (Table 1) were analyzed individually to determine whether they are enriched or depleted for any of 36 types of introns (see Materials and Methods). The 36 types of introns are written in the form X_a–X_b, where X is one of the seven types of exon, and the X_a exon precedes the X_b exon in the same gene. For example, the class A5-A3 is an intron that is flanked by an upstream exon with an alternative 5′ splice site and a downstream exon with an alternative 3′ splice site. In Figure 1A, all 4 types of splicing, indicated with dashed lines, would generate introns in the A5-A3 class. Figures 1B–D show R-R, S-S, and A3-S classes of introns and the consensus sequences that are most significantly enriched for these classes of introns. Notice that there are only 36 possibly combinations for the 7 types of exons rather than 49 (i.e., 7² = 49) because alternative poly A (APA) is never first (X_a) and alternative promoter (APr) is never second (X_b) in the X_a-X_b nomenclature system.

In an attempt to concisely summarize our analyses, Table 2 shows an “alternative splicing code” for human introns in a format similar to the genetic code. The left column indicates the exon type upstream on the intron (APr, 3S, 5S, ME, R, and S) and the top row indicates the exon type downstream of the intron (3S, 5S, ME, R, S, and APA). The numbers in the table show the intron consensus sequences that are most significantly enriched (upper cell) and depleted (lower cell) for the indicated class of intron. Supplementary Figure 1 shows splicing diagrams of the 18 classes of introns that have consensus motifs that are over-represented for that intron class, and Supplementary Figure 2 shows splicing diagrams of the 16 classes of introns that have consensus motifs that are under-represented for that intron class. Notice that many of the classes of introns have more than one over-represented motif and more than one under-represented motif (Table 2), indicating that the alternative splicing code, much like the genetic code, is degenerate. Also, notice that almost half of the intron types are not associated with consensus sequences, such as A3–A5, possibly because macromolecular complexes presumably do not exist that recognize certain rare types of alternative RNA splicing event in humans.

Six of the 36 types of introns are enriched in the non-canonical core consensus sequence gc_ag, including APr-A3 (gc_ag, gcaag_ag, and gcaagu_ag), APr-R (gc_ag), A5-A3 (gc_ag and gcagg_ag), A5-R (gc_ag, gcaag_ag, gcaagu_ag, and gcagg_ag), A5-S (gcaug_ag), and S-R (gcaug_ag) (Table 2). However, the only intron type that is depleted in a non-canonical core consensus sequence is R-R, which is depleted for gc_ag, gcaag_ag, gcaagu_ag, and gcagg_ag (Table 2). We interpret this as indicating that many types of alternative splicing events utilize the non-canonical core consensus sequence gc_ag, but that invariant splicing almost always uses the canonical core consensus sequence gu_ag. These findings are consistent, with slight variations indicated below, in all 50 species studied. A complete list of paired consensus sequences is available upon request, as well as enrichment (up) and depletion (down) p-values for the 36 types of introns for each paired consensus sequence, for all 50 organisms' genomes that we analyzed (see Supplementary Table 2 for a partial list).

In our analyses of the 50 species, the most frequent intron class in most of the species is R-R, (e.g., 79% for H. sapiens, 29% for M. musculis, 62% for D. melanogaster, 80% for C. elegans and 91% for A. thaliana) which means a invariant exon is followed by another invariant exon (Figure 1B). The second most common intron class, in most of the 50 species analyzed, is S-S (e.g., 5% for H. sapiens, 22% for M. musculus, 4% for D. melanogaster, and 3% for C. elegans), which means that two consecutive exons are skipped, either together or individually, in mature RNA (Figure 1C). Other studies, have suggested that exon skipping is the most frequently occurring alternative splicing event. For example, it was found that over one third of exons can be skipped (~38%) (Ast, 2004; Sugnet et al., 2004). “Pathological” exon skipping is commonly seen in diseases with multiple disrupted alternative splicing events, especially in cancer (Watson and Watson, 2010).

There are many practical uses for understanding the alternative RNA splicing code. For example, many diseases, including cancer, have mutations that cause changes in alternative RNA splicing that contribute to pathogenesis (Watson and Watson, 2010). It is estimated that at least 15–50% of mutations that cause human diseases affect splice-site selection (Wang and Cooper, 2007; Singh and Cooper, 2012). Here we show how the alternative RNA splicing code in Table 2 helps to interpret human genetic diseases that are caused by mutations near splice donor and acceptor sites that could not be adequately explained without this code. Using the databases of disease-causing mutations at spliced 3′ and 5′ splice sites, dbass5 and dbass3 (http://www.dbass.org.uk/dbass5/viewlist.aspx; Singh and Cooper, 2012), we analyzed all intron mutations at intron positions +3, +4, +5, and +6 (the first intron nucleotide at the splice donor is +1) and successfully correlated the alternative RNA splicing code to 96 different mutations in 56 genes (Supplementary Table 3). For example, Menkes disease (MD), which has several alleles in the ATP7A gene that are associated with alternative splicing defects, is a lethal disorder of copper metabolism that lead to severe neurological degeneration (Møller et al., 2000). Occipital horn syndrome (OHS) is a milder allelic form that is caused by partial loss of function of the ATP7A gene (Møller et al., 2000). Both MD and OHS are caused by mutations in the intronic sequences of the ATP7A gene, which encodes an ATPase that is responsible for copper efflux from cells (Figure 2A; Nissim-Rafinia and Kerem, 2002).

In the ATP7A gene, two splice-site mutations (IVS6+1G>A, IVS6+5G>A) for MD and one (IVS6+6T>A) for OHS were identified in a previous study (Figure 2; Møller et al., 2000). The main biological effect of the mutation in the first position of the splice donor site of intron 6 (gu to au) is cryptic downstream splice donor usage followed by exon 7 skipping (Figure 2B; Møller et al., 2000). Exon skipping and cryptic splice site activation are typical results of mutations in any of the four core consensus bases, gu_ag, and can be explained without the alternative splicing code. However, why the ATP7A mutation in position 5 of intron 6 (IVS6+5G>A) has such a severe effect on alternative splicing was previously not understood since this is outside of the canonical gu_ag consensus sequence (Figure 2C; Møller et al., 2000).

Using the alternative RNA splicing code, we can now better explain the alternative splicing phenotypes caused by the mutations the 5th position of the 5′ splice site of intron 6 of ATP7A. The wild-type sequence guaagu_ag corresponds to a paired consensus sequence that is over-represented for R-R, which means that there is little or no alternative splicing in the wild-type ATP7A gene for this intron (Figure 1A). However, the 5th position mutation (Figure 2C) corresponds to the guaa_ag paired consensus sequence that is over-represented for the intron class S-S (Figure 1C). Therefore, the alternative RNA splicing code helps explain why two adjacent exons, exons 6 and 7, are skipped as the result of the mutation in the 5th position (Figure 2C). A similar argument can be made for the milder ATP7A mutation in OHS, (IVS6+6T>A), which leads to a motif change to guaag_ag, which is an over-represented motif for the intron class A3-S, and leads to incomplete exon 6 and/or exon 7 skipping and cryptic splice site usage 50 nucleotides downstream of the normal 5′ splice site in intron 6, at a second guaag_ag sequence (Figure 2D). In the OHS allele, exon 6 becomes an A3 exon because the 5′ splice site of exon 5 can join with the normal 3′ splice site or exon 6 or the alternative 3′ splice sites of exon 7 or exon 8 (Figure 2D).

We note that the above analysis for ATP7A intron 6 is an over simplification of what is required to predict the effect of an intron mutation because multiple consensus sequences are often enriched or depleted in several of the 36 types on introns. For example, the wild type ATP7A intron 6 consensus sequence, guaagu_ag, corresponds to a consensus sequence that is enriched for R-R, R-S, and S-APA. Therefore, in order to predict the outcome of a mutation in a consensus sequence, one must determine which intron classes are uniquely enriched when a mutation is present that was not enriched in the wild-type sequence. The sixth position mutation in OHS has the intron sequence guaag_ag which is enriched in A3-S and R-R. This might explain why both A3-S and R-R splicing events are induced by the OHS mutation (Figure 2D). Similarly, the fifth position mutation in MD2 has the sequence guaa_ag, which is only enriched in the intron type S-S. This might explain why S-S splicing events are induced by the MD2 mutation (Figure 2C). The alternative RNA splicing code can also be used to explain +3 to +6 intron mutations in neurofibromatosis type 1 (NF1), one of the most prevalent inherited disorders in human (Hastings and Krainer, 2001), beta thalassemia (HBB) (Felber et al., 1982), and many other human diseases.

In addition to the canonical splicing pathway, which uses the gu_ag consensus sequence, there are non-canonical (a.k.a., minor) splicing pathways that sometimes do not use the gu_ag consensus (Padgett, 2012). The canonical splicing pathway generally uses the U1 and U2 small RNAs in their splicing mechanism, always at gu_ag introns, while the non-canonical pathway uses U11 and U12 small RNAs, at both gu_ag and au_ac introns. The U12-like introns also have several conserved nucleotides that flank the splice donor and splice acceptor sequences (Padgett, 2012). When we searched for U12-like consensus sequences in the lists of intron consensus sequences, we found that human and mouse share the top three U12-like sequence matches: (1) guauccuuu_ag (Rank 29, Table 1), (2) auauccuu_ac (Rank 37, Table 1) and (3) guauccuu_ag (Rank 36, Table 1). The U12-like motif guauccuuu_ag is also the best match with the U12-like splicing pathway in A. thaliana (Supplementary Table 2). Curiously, both D. melanogaster and C. elegans have the weakest matches to the U12-like splicing sequence, gugggu_cag, and guucguuuuu_uuucag, respectively, even though they are presumably evolutionarily closer to humans than plants (Supplementary Table 2).

As we showed with mutations that affect the major splicing machinery, mutations that affect the minor splicing machinery can also be better interpreted with the paired consensus sequence motifs that we identified. One example involves a tumor suppressor gene, LKB1, whose splice acceptor mutation in the second intron is thought to cause Peutz-Jeghers syndrome (PJS) (Hastings et al., 2005). This mutation changes the splice junction sequence from auauccuu_ac to guauccuu_ac, and causes aberrant splicing, even though the mutation is changing a non-canonical “au” splice donor to a canonical “gu” splice donor (Figure 3A). Perusing the alternative RNA splicing code, we noticed that the wild-type LKB1, auauccuu_ac, is present, but the sequence found in PJS, guauccuu_ac, is not present on the paired RNA splicing consensus sequence table in humans (Table 1). Therefore, even though the consensus sequence table indicates that the splice donor sequence guauccuu is a good minor splice donor sequence, the paired-sequence analyses indicate that the “gu” core splice donor sequence must be paired with another canonical splice acceptor sequence, “ag,” even in U12-type introns. In other words, our analyses suggest that there are at least two distinct classes of U12-type introns in humans; one with the core sequence gu_ag and the other with au_ac, and the machinery that recognizes the two ends of the introns in the U12-type splicesosomes cannot be swapped. This hypothesis might also help explain the unusual splicing reactions at the 3′ splice site to be multiple cryptic dinucleotide termini (such as cg, au, ug, and gg) observed from different patients since no “ag” is present in vicinity of the splice acceptor site (Figure 3B; Hastings et al., 2005).

Here, we admit that interaction through direct pairing of both the 5′ and 3′ splice-site recognition complexes is one of the possible explanations that affect splicing events. In our hypothesis for developing the alternative splicing code, we emphasized that the 5′ and 3′ splice sites should be considered as a pair. However, mostly these sites can be recognized by independent complexes, 5′ splice sites by U1 and 3′ splice sites by U2AF in the first step of canonical splicing. In addition, in most of the consensus sequences we identified, the 3′ splice sites are far less variable than the 5′ splice sites. Other possibilities that could explain the alteration of splicing events which need to be further investigated include interaction of splicing enhancers or silencers.

In addition to the long consensus sequences identified in the U12-type non-canonical introns by rank-order analyses (Table 1), we found several uncommonly long non-U12-type consensus sequences (>8 nucleotides) in humans that are also conserved in mouse, but not C. elegans, D. melanogaster, and A. thaliana. Splice junctions that are targeted by antisense oligonucleotides can often block alternative splicing, such as “morpholino” oligonucleotides that are used to inhibit splicing (Morcos, 2007). Therefore, one possibility is that the long consensus sequences might be targets for long-noncoding RNA or RNA-binding proteins that regulate alternative RNA splicing in these introns. To test this idea, we used the Mirbase database (http://www.mirbase.org) to search for complimentary miRNA sequences for the long human and mouse alternative mRNA splicing motifs. We found that the long paired alternative RNA splicing motif gu_uuugagacag (Rank 42, Table 1) is an excellent match with human hs-miR-6510-3p as well as mouse mmu-miR-706 (Figure 4A). This bioinformatics result suggests that this uncommonly long splice site consensus sequence might be negatively regulated by one of these microRNAs, but this hypothesis requires genetic validation (Figure 4B).

We also compared alternative splicing consensus sequences among the 50 species to determine whether they are evolutionarily conserved and whether they are enriched in the same classes of introns. The top two consensus sequences that are shared by the greatest number of species are gugagu_ag and gu_ag (Figure 5). The motif gugagu_ag is enriched for the intron class A5-A3 in 10 of the 50 species (Figure 5A), and the motif gu_ag is enriched in the intron class A5-S in 11 of the 50 species and depleted in the intron class A5-S in 4 of the 50 species (Figure 5B). Humans and mouse share 80% of all alternative RNA splicing motifs. When we looked for a possible reason why D. melanogaster, C. elegans and A. thaliana share only small portion of significant motifs with human (14, 26, and 29%), we found that, although the canonical sequence gu_ag is the most highly conserved (98%), the third base after the splicing donor “gu” varies among taxa. The base adenine was hardly ever observed in the third position of the intron donor sequence in D. melanogaster or C. elegans (<1%), while adenine is the most common nucleotide in the third position in the splice donor, i.e., gua_ag, for both human (58.4%) and mouse (58.3%).

Next, we used the program High-Throughput GoMiner (HTGM) (Zeeberg et al., 2005) to analyze and interpret the functional significance of conserved gene sets that have at least one intron with identical paired consensus sequences between humans and A. thaliana. In this analysis, significant Gene Ontology (GO) categories were determined for genes with a single particular type of paired consensus sequence in an intron. Not surprisingly, since they are separated by hundreds of millions of years of evolution, not many paired alternative splicing motifs are conserved between humans and plants. However, the paired motif gugagu_ag is the third most common motif in A. thaliana (i.e., Rank 3) and the most common motif in humans (Rank 1). Similarly, gua_ag is the most common motif in A. thaliana and the 6th most common motif in humans. Genes that have these two sets of motifs in A. thaliana and humans cluster together in the HTGM analysis with the GO categories “protein phosphorylation,” “protein tyrosine kinase,” “tRNA metabolic processes,” “tRNA aminoacetylation,” and “amino acid activation” (Figure 6A). The identification of GO categories in protein phosphorylation is interesting because species-specific alternative exons were analyzed in a recent study, and the two most enriched keywords are “alternative splicing,” which reflects the abundant alternative RNA splicing in this set of genes, and “phosphorylation” (Merkin et al., 2012). Alternative RNA splicing is used to alter protein phosphorylation, which can alter protein stability, subcellular localization, and activity (Merkin et al., 2012). Our results suggest that using alternative RNA splicing to regulate phosphorylation of proteins might be an evolutionarily conserved mechanism in plants and humans and is regulated by these two shared alternative RNA splicing motifs (Figure 6A).

We also did HTGM analyses of the gene sets with identical paired alternative RNA splicing sequences in H. sapiens, M. musculis, C. elegans, and Drosophila, but excluding A. thaliana because it is too divergent from the other species. Multiple-species cluster analyses of Gene Ontology categories determined by HTGM indicates that “cell morphogenesis” and “cell development” and related GO categories are conserved across the four animal species for multiple similar consensus sequences (Figure 6B; for a full HTGM analysis of the four animals see Supplementary Figure 3). Our HTGM result in the four animal species is consistent with previous studies demonstrating alternative splicing is often evolutionarily conserved in a tissue and developmental stage-specific manner (Chen et al., 2012; Merkin et al., 2012).

Finally, we wanted to determine whether genes with multiple introns tend to have the same or different intron consensus sequences. If, as we hypothesize, the 36 types of introns in the alternative splicing code table (Table 2) utilize many different types of splicesosomes, then it is possible that genes with multiple introns will utilize as few different types of splicesomes as possible. It might not be practical for different types of splicesomes to splice each of the different introns in a gene. Rather, since RNA splicing occurs concurrently with transcription, it might be more efficient for a particular type of spicesosome to move on to the next intron once it has completed the splicing reaction of the upstream intron.

The most significant cluster of intron types with genes with multiple introns corresponds to U12-like introns (rank 36, 51–52 in mouse) and the second most significant cluster corresponds to U2-like introns that correspond to the top 12 ranked consensus sequences (Figure 7). Introns in ranks 13–26 form a third cluster and introns in ranks 27–50 (excluding rank 36, which is a U12-type intron) form a fourth cluster, but this is not well segregated from the third cluster (Figure 7). We interpret this as suggesting, using the logic in the previous paragraph, that there might be as many as 3 or 4 different types of splicesosomes in the mouse. We performed similar cluster analyses in Arabidopsis, C. elegans, Drosophila, and humans and identified 3 or 4 similar clusters of intron consensus sequences in genes with multiple introns (Supplementary Figures 4–7). In humans, rank 1–42 introns all have remarkably similar distributions in genes with 1–14 introns. For example, over 50% of all genes have 6–10 introns, and over 50% of all the rank 1–42 introns are in genes with 6–10 introns (Supplementary Figure 8).

The alternative mRNA splicing code in Table 2 provides a new toolkit for characterizing cis-acting sequences that are important for generating the enormous diversity of processed mRNAs in animals and plants. Other groups have attempted to decipher the alternative mRNA splicing code by characterizing the binding sites for a large number of RNA binding proteins that are known to affect tissue-specific alternative mRNA splicing, such as MBNL, PTB, RBFOX, STAR, and TIA families of splicing factors (Ladd and Cooper, 2002). The research group of Xiong et al. (2015) have extracted human DNA sequences as input and used a Bayesian model to predict the percentage of exons that are spliced in (PSI) based on sequences within the transcripts. For example, they found among intronic variants that are known to cause disease alter splicing nine times as often as the common variant when they are more than 30 nucleotides away from any splice site. We believe that our study complements the study by Xiong et al. because their program can be used to calculate PSI-values, whereas, our method will help to predict the type of alternative splicing that SNPs with altered PSI-values produce. We did not investigate branchpoints in our paper, but Mercer et al. (2015) analyzed RNA-seq data with novel bioinformatics methods to identify 60,000 high-confidence human branchpoints.

Combining our alternative mRNA splicing consensus sequence information with the RNA-binding protein datasets should allow the construction of better RNA splicing maps (Witten and Ule, 2011) that can be used to better understand the mechanism of tissue-specific alternative mRNA splicing events. The alternative RNA splicing code can also be used to better understand how human germline disease mutations and somatic mutations in cancer affect alternative RNA splicing and lead to disease etiology. Future biochemical experiments are needed to test the hypothesis that the many classes of paired alternative RNA splicing events in humans with paired consensus sequences have unique macromolecular complexes that regulate RNA maturation. Future bioinformatics analyses are needed to predict how a particular splice site mutation in any of the first or last few nucleotides in an intron precisely affects alternative splicing. The alternative splicing code should help inform both of these endeavors.

DR is the PI of the lab and directed the analyses and the writing of the manuscript. WQ helped write the paper, made the figures, and performed the human SNP analyses. PC wrote the software and came up with the original idea, BZ did the multiple species GO analyses.