Role and convergent evolution of competing RNA secondary structures in mutually exclusive splicing

ABSTRACT

Exon or cassette duplication is an important means of expanding protein and functional diversity through mutually exclusive splicing. However, the mechanistic basis of this process in non-arthropod species remains poorly understood. Here, we demonstrate that MRP1 genes underwent tandem exon duplication in Nematoda, Platyhelminthes, Annelida, Mollusca, Arthropoda, Echinodermata, and early-diverging Chordata but not in late-diverging vertebrates. Interestingly, these events were of independent origin in different phyla, suggesting convergent evolution of alternative splicing. Furthermore, we showed that multiple sets of clade-conserved RNA pairings evolved to guide species-specific mutually exclusive splicing in Arthropoda. Importantly, we also identified a similar structural code in MRP exon clusters of the annelid, Capitella teleta, and chordate, Branchiostoma belcheri, suggesting an evolutionarily conserved competing pairing-guided mechanism in bilaterians. Taken together, these data reveal the molecular determinants and RNA pairing-guided evolution of species-specific mutually exclusive splicing spanning more than 600 million years of bilaterian evolution. These findings have a significant impact on our understanding of the evolution of and mechanism underpinning isoform diversity and complex gene structure.

KEYWORDS:

• Competing RNA secondary structures
• convergent evolution
• MRP1
• Mutually exclusive splicing

Introduction

Alternative pre-mRNA splicing is a major mechanism for enhancing transcriptomic and proteomic diversity and regulating complex gene expression patterns during evolution.^1-4 Exon duplication has been an important source of novel alternative splicing and gene function, and has tended to be involved in mutually exclusive alternative splicing.⁵ Mutually exclusive splicing is a very strictly regulated form of alternative splicing, which is characterized by the inclusion of only one of 2 or more candidate exons in the mature mRNAs.⁶ The most extreme example of mutually exclusive splicing is the insect Down syndrome cell adhesion molecule (Dscam) gene, which contains multiple clusters of variable exons and can potentially produce 38,016 different isoforms in Drosophila melanogaster.⁷ Mutually exclusive splicing plays important regulatory roles in a myriad of cellular and developmental processes.^8-10 Mutations in variable exons and related splicing elements may contribute to diverse diseases, including Timothy syndrome,¹¹ cancer,¹² spinal analgaesia,¹³ cardiomyopathy,¹⁴ and heart disease.¹⁵

Several mechanisms have previously been proposed to guarantee the mutually exclusive use of variable exons that are arrayed in tandem, including steric interference between splice sites,¹⁶ spliceosome incompatibility,¹⁷ nonsense-mediated decay (NMD),⁶ and competing RNA secondary structures.^18-20 The former 3 mechanisms may explain the mutually exclusive usage of no more than 2 variable exons, such as mammalian α-tropomyosin,¹⁶ Caenorhabditis elegans 14–3–3ξ,²¹ and Dscam exon 17.²² The competing RNA secondary structures function in controlling the mutually exclusively splicing of 2 or more variable exons, including Drosophila Dscam exon 4, 6, and 9 cluster, 14–3–3ξ, mhc.^18-20,23 Recently, an analogous model based on the competition of upstream and downstream intronic RNA secondary structures that evolved in Drosophila srp and hymenopteran Dscam exon 4 and lepidopteran Dscam exon 9 clusters has been proposed.²⁴ Additionally, bioinformatics studies have suggested that RNA pairing may be involved in the mutually exclusive splicing of mammalian dynamin 1.²⁵ However, these RNA pairing systems have not been identified in non-arthropod species, and it is unclear whether different phyla rely on similar strategies. Moreover, no vertebrate genes with more than 2 mutually exclusive exons have been identified to date,¹⁸ suggesting a specific evolutionary constraint of mutually exclusive splicing in vertebrates.

The ABC gene family is now recognized as one of the largest transporter families in all kingdoms of life.²⁶ Multidrug resistance proteins (MRP), a large group of ABC transporters, are evolutionarily conserved throughout the animal kingdom.²⁷ In addition to having canonical roles in the cellular export and removal of toxic metabolites and compounds, the MRPs likely play critical roles in a diverse array of cell physiologic and disease processes.²⁷ MRPs are also thought to be essential parts of the ancient immune defense system, working in “cellular antitoxic” mechanisms.²⁸ Although the duplication and loss of MRPs during metazoan evolution have been well studied, less is known about the evolution of alternatively spliced isoforms. In D. melanogaster, MRP1 contains 10 alternatively spliced exons organized into 2 clusters; these can potentially generate 16 different isoforms via mutually exclusive splicing.²⁹ The C. elegans MRP homolog contains 4 mutually exclusive exons, which could regulate larval development.³⁰ These MRP splice variants are largely lineage-specific, possibly representing a strong and flexible system capable of adapting rapidly to multiple environmental challenges. Therefore, the evolutionary dynamics of the MRP exon cluster across different taxa are of great interest in terms of the evolution, mechanism, and function of complex mutually exclusive splicing.

In this study, we found that tandem exon duplication was deeply conserved in bilaterian MRPs. Interestingly, these events were of independent origin in flies (D. melanogaster), worms (C. elegans), and fish (Danio rerio), suggesting the convergent evolution of alternative splicing. However, mammalian MRP1 contains only one exon copy, although MRP1 underwent bursting exon duplication in the early-diverging chordate Branchiostoma lineage. Moreover, we showed that dual sets of evolutionarily conserved RNA pairings direct mutually exclusive splicing of the MRP1 pre-mRNA in insects. Importantly, we also discovered and verified a similar structural code in MRP exon clusters of the annelid, C. teleta, and the chordate, B. belcheri, suggesting a similar strategy for regulating mutually exclusive splicing across bilaterians. Taken together, these data provide a remarkably detailed picture of RNA pairing-guided evolution of species-specific mutually exclusive splicing spanning more than 600 million years of bilaterian evolution. Our findings also provide new evolutionary and mechanistic insights into isoform diversity and complex gene structure.

Results

Tandem exon arrays are deeply conserved in metazoan MRPs

To examine the origin and evolution of tandem exon arrays in metazoan MRPs, we analyzed exon duplication in major phyla, including Spongia, Cnidaria, Nematoda, Platyhelminthes, Annelida, Mollusca, Arthropoda, Echinodermata, and Chordata. Similar to Drosophila and Caenorhabditis, MRP1 underwent tandem exon duplication in roughly identical regions in most bilaterian species, including Nematoda, Platyhelminthes, Annelida, Mollusca, Arthropoda, Echinodermata, and Chordata. The number of exon copies varied considerably across species, ranging from 2 to a maximum of 11 (Fig. 1A; Figs. S1 and S2). In Drosophila, the hypervariable regions encoded by the duplicated exons roughly corresponded to the fourteenth and fifteenth transmembrane regions (Fig. 1B and C). Alternative splicing of the duplicated exons could modify the transport properties of the protein.³¹ However; the MRP genes do not contain duplicated exons in early-diverging metazoans, such as Hydra vulgaris, Nematostella vectensis, and Amphimedon queenslandica. This suggests that exon duplication occurred before the divergence of protostomes and deuterostomes, which occurred ∼600 million years ago.³²

Figure 1. Exon duplication is conserved across invertebrate and vertebrate MRPs. (A) Origin and evolution of exon duplication in MRPs. This phylogenetic tree illustrates the relative positions of the organisms examined here and is provided as a quick reference for taxonomic relationships.⁴⁶ Exons that are duplicated in tandem are shown in different colors, based on the size of the duplicated exon. Exons that are the same color are homologous. Lost exons are indicated by a dashed box. Independent exon duplication events are indicated by solid circles. Phylogenetic analysis indicates that exon duplication may have occurred before the split of the Bilaterians, over 600 million years ago. (B) Schematic illustration of MRP domain. The regions in yellow are encoded by alternative exons. Brown squares indicate the Walker A and Walker B motifs of the ABC superfamily. Pentagrams indicate the positions of thiolester motifs. (C) Tertiary structure model of MRP. The hypervariable regions (shaded in yellow) are encoded by alternative exons. (D) The independent origin of MRP exon duplications in the basal chordate Branchiostoma, nematodes, and insects. See Table S1 for explanation of abbreviations.

Interestingly, both the number and structure of MRP1 duplicated exons diverged considerably during chordate evolution. We found that 7 and 10 duplicated exons were present in MRP1 homologues of the basal chordates B. belcheri and Branchiostoma floridae, respectively (Fig. 1; Fig. S1). Multi-exon duplication events have been detected in some Cionidae species, which underwent pairwise mutually exclusive splicing (Fig. 1; Fig. S1). This gene structure is conserved in Ciona intestinalis, Ciona savignyi, and Botryllus schlosseri. These sequences, which are encoded by a single alternative exon in most species, have been split into 2 exons in Cionidae and then further duplicated through exon expansion. Additionally, 2 tandemly arrayed exons are present in the Pisces MRP1, such as zebrafish. However, MRP1 contains only a single copy of the exon in late-diverging vertebrates, including humans, lizards, chickens, and Xenopus tropicalis (Fig. S1). Phylogenetic analysis indicated that the second duplicated exon ortholog of zebrafish MRP1 was lost in Tetrapoda MRP1 at the divergence of the Amphibia and Pisces classes.

Evolutionary convergence of exon duplication in different phyla

Next, we analyzed the relationships of MRP variable exons in various phyla. Comparison of sequences of the duplicated exons indicated that the 2 copies from insect MRPs resembled each other more closely than they resembled any of the duplicates from the Branchiostoma or Caenorhabditis MRPs (Fig. 1D). Moreover, although these variable exons belonged to a generally similar domain, their intron/exon structures differed, albeit with the same phase. The 5′ junctions of the variable exon in Branchiostoma MRP corresponded to those of Drosophila MRP, whereas their 3′ junctions corresponded to those of Caenorhabditis MRP. These 2 lines of evidence indicate that exon duplication events occurred independently in various phyla of arthropods and nematodes as well as of chordates. This suggests that these alternative splicing outcomes may have intrinsic functional benefits. Rescue experiments also revealed that both the b-isoform and the c-isoform, but not the e-isoform, can rescue the mutant phenotype, suggesting that this may be due to a discrepancy in substrate specificity among the variant isoforms.³⁰ Therefore, these splice isoforms provide novel functional properties, probably in the recognition of or affinity for different substrates.

Dual RNA pairings control the exon 8 selection of lepidopteran MRP1

To elucidate the mechanism underlying the mutually exclusively splicing of duplicated exons, we first investigated the evolutionary dynamics of variable exons in arthropods, particularly in insects. Two species of the order Crustacea (Daphnia pulex and Daphnia magna) and 68 species representing 4 major clades of the order Insecta were examined. Together, these organisms allowed us to examine MRP1 evolution over ∼450–500 million years, with both short-term resolution (comparison with Drosophila, Anopheles, Lepidoptera, or Hymenoptera) and long-term resolution (comparison among the 4 insect clades and Daphnia).

We first searched for sequences that were conserved among lepidopteran species, because all contain only 2 exon 8 variants. Interestingly, comparative sequence analysis identified 4 highly conserved intronic elements, one in the intron downstream of exon 7, 2 between exon 8.1 and exon 8.2, and one upstream of exon 9 (Fig. 2A). These elements could potentially form 2 pairs of RNA secondary structures between introns in each species (Fig. 2B; Fig. S3). To examine the effects of dual RNA pairings on exon 8 selection of lepidopteran MRP1, we tested disrupting and compensatory mutations in the dsRNA stem of the silkworm MRP1 pre-mRNA (Fig. 2B). Mutating the upstream docking sequence (Ud) or the downstream selector (Ds) could abolish inclusion of exon 8.2 (Fig. 2C–E). Importantly, the reduced efficiency of exon 8.2 inclusion due to individual mutations reverted to wild-type (WT) levels through compensatory double mutations. In addition, the branch points are predicted to be located downstream of selector sites. Together, these data indicated that inclusion of exon 8.2 is activated by upstream RNA pairing (stem I) in the silkworm.

Figure 2. Dual RNA pairing controls alternative exon 8 inclusion in Bombyx mori MRP1. (A) Genomic organization of B. mori (Bmo) MRP1. The upstream docking site (Ud), upstream selector sequence (Us), downstream docking site (Dd), and downstream selector sequence (Ds) are indicated. Above are sequences of consensus intronic elements in Lepidoptera species (for abbreviations used, see Table S1), with spacings as indicated (nt) below. Color shading indicates the most common nucleotide at each position. (B) Predicted RNA pairing of B. mori MRP1 pre-mRNA. Mutations introduced into the dsRNA stem are indicated on the left or right as mutated sequences (M1–M4). The RNA pairings for other Lepidopteran species are shown in Fig. S3. The gray dashed arrow represents potential RNA-RNA interaction. The green arrow depicts activating the inclusion of the proximal exon while the dashed arrow indicates the potential activation. (C) The effects of mutations on the exon 8 inclusion, including disruptive single mutations (M1–M4) and compensatory double mutations (M12 and M34). WT, wild type. RNA for RT-PCR was isolated from silkworm cells that were transiently transfected with the plasmids indicated. The band marked with an asterisk (*) is a non-specific RT-PCR product. (D) Effects of mutations on exon 8.1 and exon 8.2 selection. The PCR lanes represent a single band cut from panel C, containing exon 8. Exon 8.2-specific restriction digestion was performed to determine the frequency of exon 8.2 utilization. (E) Effects of dual RNA pairings on exon 8 selection in B. mori. The exon 8.2 selection frequency was quantified based on the data in panels C and D (exon 8.2 inclusions/total inclusions). Data are expressed as a percentage based on the mean ± SD of 3 independent experiments. *P < 0.05 (Student's t-test).

We next attempted to test the effect of the dsRNA stem II of the silkworm MRP1 on exon 8.1 inclusion using a similar strategy (Fig. 2B). As silkworm exon 8.1 was infrequently included in minigene constructs with BmN cells, we have not demonstrated the proposed effect of dsRNA stem II on exon 8.1 inclusion. However, mutating the upstream selector (Us) or the downstream docking site (Dd) markedly enhanced exon 8.2 inclusion (Fig. 2C–E). The increased efficiency of exon 8.2 inclusion due to individual mutations was reverted to WT levels through compensatory double mutations (Fig. 2C–E). These results were concomitant with repression of the in-loop exon by RNA pairing³³ and thus strongly verified the presence of RNA stem II.

Additionally, we observed that the intronic distance between the selector element and the 3′ splice site was consistently small (45–75 nt) among all these species except Plutella xylostella, whereas the distance between the 5′ splice site and the selector element was large and varied among species (38–1506 nt) (Fig. 2A). As the spacing between the splice site and the selector is important for exon activation,²⁰ we speculated that the remarkable distance discrepancy may have contributed to the different inclusion patterns of exon 8.1 and exon 8.2. To confirm this suggestion, we examined the role of the RNA pairing of P. xylostella MRP1, in which the distance between the 5′ splice site and the upstream selector is similar to that between the downstream selector and the 3′ splice site (Fig. S4A and B). Disrupting and compensatory mutations revealed that exon 8.2 inclusion was activated by upstream RNA pairing, whereas exon 8.1 inclusion was enhanced by downstream RNA pairing in P. xylostella (Fig. S4C–E). Taken together, these results demonstrate that dual RNA pairing evolved to control exon 8 inclusion in Lepidoptera. Since these dual RNA pairings are structurally similar to the srp gene,²⁴ an analogous regulatory mechanism may control mutually exclusive splicing.

RNA pairing-guided evolution of mutually exclusive splicing in Hymenoptera

Next, we examined whether the exquisite architecture found in Lepidoptera is conserved in Hymenoptera. We found that only 2 mutually exclusive exons exist in bee species, whereas 3 copies are present in closely related ants, and there are 4 copies in distantly related Pteromalus puparum and Athalia rosae (Fig. 3A; Fig. S2). Phylogenetic analysis indicated that exon 8.1 and exon 8.3 were sequentially lost in ancestral MRP1 during hymenopteran evolution, nearly 100 and 200 million years ago, respectively (Fig. 3A).³⁴ How did MRP1 evolve to guarantee species-specific mutually exclusive behavior during hymenopteran speciation? Remarkably, 4 intronic elements within exon cluster 8 were found to be conserved among 15 hymenopteran species, with estimated divergence times ranging from 10 million to 200 million years ago (Fig. 3A). Similar to the observations for Lepidoptera, these Hymenoptera-specific intronic sequences potentially form 2 pairs of RNA secondary structures (Fig. 3B; Fig. S5). Therefore, we concluded that these intronic elements are clade-specific, but with the potential to form conserved RNA architectures.

Figure 3. The evolution of mutually exclusive splicing in Hymenopteran MRPs is guided by RNA pairing. (A) Phylogenetic arrangement of the cis intronic elements flanking mutually exclusive exons during Hymenoptera evolution. Symbols are as described in Fig. 2. A schematic diagram of the partial pre-mRNA with constitutive exons depicted as black boxes, alternative exons as colored boxes, and introns as lines. Above are sequences of consensus intronic elements for different taxa (for abbreviations used, see Table S1) with spacings as indicated (nt). The curved lines represent RNA-RNA interactions. (B) The predicted intronic RNA pairing of A. mellifera (Ame) pre-mRNA. Mutations introduced into the dsRNA stem are indicated (M1–M4). The RNA pairings for other Hymenopteran species are shown in Fig. S5. (C) Effects of mutations on the exon 8 inclusion are indicated for disruptive single mutations (M1–M4) and compensatory double mutations (M12 and M34). WT, wild type. The band marked * is a non-specific RT-PCR product. (D) Effects of mutations on the exon 8 inclusion. (E) Predicted intronic RNA pairing of P. puparum (Ppu) MRP1 pre-mRNA. Mutations introduced into dsRNA stem are indicated (M5, M6). (F,G) Effects on exon 8 inclusion (F) exon 8 selection (G) of disruptive single mutations (M5 and M6) and compensatory double mutations (M56). Data are expressed as a percentage representing the mean ± SD from 3 independent experiments. *P < 0.05 (Student's t-test).

Similar to the silkworm, structure-restoring double mutations in dsRNA stem I restored the efficiency of Apis mellifera exon 8 inclusion to a degree similar to that observed with the WT minigene (Fig. 3B–D). These combined data indicated that the inclusion of A. mellifera exon 8.2, which is orthologous to silkworm exon 8.2, was activated by upstream RNA pairing (stem I) in A. mellifera, which is similar to the case in silkworms. Although A. mellifera exon 8.1 is not orthologous to silkworm exon 8.1, a pattern similar to that of silkworm stem II was obtained by disruptive and compensatory mutations for A. mellifera stem II (Fig. 3C and D). These results indicated that, as in Lepidoptera, dual RNA pairings control exon 8 inclusion of A. mellifera MRP1.

As the ant exons 8.2 and 8.3 are orthologous to A. mellifera exons 8.1 and 8.2, we speculated that dual RNA pairings (Stem I and Stem II) control their inclusion in ants. Next, we examined how ant and P. puparum exon 8.1 splicing is controlled. We hypothesised that if the selector sequences required to ensure mutually exclusive splicing of the duplicated exon duplicated simultaneously with exon expansion via staggered homologous recombination,³⁵ downstream RNA pairing may direct its paralogue exon 8.1, as it does for exon 8.2. However, our analyses failed to find a selector with reverse complementarity to Dd1. Interestingly, a conserved intronic element was found in the intron downstream of exon 8.1 among 10 non-bee hymenopteran species, which was complementary to a conserved intronic element (Dd2) about 100 bp downstream of Dd1 (Fig. 3A). These predicted RNA pairings (Stem III) are highly conserved with a common core among 9 non-bee hymenopteran species (Fig. S5B). Further disruptive and compensatory mutations experimentally confirmed this proposed Stem III (Fig. 3E–G; Fig. S6). Moreover, Dd2 was absent in A. mellifera MRP1, concomitant with the loss of its ortholog of ant and P. puparum exon 8.1 (Fig. 3A). Taken together, these results indicated that the choice of alternative exon 8 was modulated by multiple sets of species-specific RNA pairings during hymenopteran evolution.

Bidirectional expansion of dual RNA pairings in Drosophila

Next, we explored how RNA pairing acts in Drosophila MRP1, which shows the largest exon 8 cluster expansion during insect evolution. The variable exon 8 number ranges from 4 in Mayetiola destructor to 11 in Drosophila virilis (Fig. S2). We found that 4 highly conserved intronic elements were restricted to the Sophophora clade of fruit flies, one in the intron downstream of exon 7, one in the intron downstream of exon 8.4, one in the intron upstream of exon 8.8, and one upstream of exon 9. These intronic sequences potentially form 2 pairs of RNA secondary structures (Fig. 4B; Fig. S7). We have not convincingly predicted selector sequences for other variable exons. However, since the selector sequences that ensure mutually exclusive splicing of the duplicated exon might be duplicated during exon expansion,³⁵ other variable exons are likely regulated via a similar dsRNA-directed mechanism. Indeed, such architectures and secondary structures are conserved between Anopheles, M. destructor, and Musca domestica, although their sequences are different (Figs. S8 and S9). To examine whether upstream RNA pairing acts in Drosophila MRP1 as it does in the silkworm and A. mellifera, we tested disrupting and compensatory mutations in the dsRNA stem downstream of exon 8.4 (Fig. 4B). Mutation in the Dd element led to a remarkable decrease in exon 8.4, but exon 8.8 inclusion was increased (Fig. 4C and D). Similarly, Us mutation specifically decreased the inclusion of exon 8.4. Curiously, a structure-restoring double mutation (Us and Dd) increased the efficiency of exon 8.4 inclusion to approximately the WT level (Fig. 4C and D). Taken together, these data indicated that exon 8.4 inclusion is dependent on base pair interactions between the downstream docking site and the upstream selector sequence in Drosophila.

Figure 4. Mutually exclusive splicing of the Sophophora MRP1 pre-mRNA is directed by bi-directional RNA pairing. (A) Genomic organization of the D. melanogaster MRP1 gene. Symbols used are the same as in Fig. 2, and the introns are not drawn to scale. (B) The predicted RNA pairings of MRP1 pre-mRNA. The RNA pairings for other Sophophora subgenus species are shown in Fig. S7. Mutations introduced into dsRNA are indicated on the left or right as mutated sequences (M1–M6). (C) Effects on exon 8 inclusion of disruptive single mutations (M1–M6) and compensatory double mutations (M12, M34, and M56). WT, wild type. The RNA was used as a template for RT-PCR and variant-specific restriction digestion to analyze the frequency of exon 8 utilization. (D) Effects of competing intronic RNA pairings on exon 8 inclusion are validated by disruptive single or combined mutations (M1–M4) and compensatory double mutations (M12 and M34). Data are expressed as a percentage representing the mean ± SD of 3 independent experiments. ***P < 0.001 (Student's t-test). (E) Effects of upstream RNA pairing on exon 8 inclusion are validated by disruptive single or combined mutations (M5 and M6) and compensatory double mutations (M56). **P < 0.01, ***P < 0.001 (Student's t-test) (F) Comparison of the selection frequency of exon 8.8 in various tissues of WT and Mrp^ΔLeft8.8. Overview of naturally-occurring MRP1 deletion mutants in D. melanogaster populations. (G) The inclusion of exon 8.8 was inhibited in natural Mrp^ΔLeft8.8 mutant isolates.

We used a similar strategy to verify the dsRNA stem upstream of the D. melanogaster exon 8.8, which is orthologous to the silkworm exon 8.2. Mutating the upstream docking sequence (Ud) or the downstream selector (Ds) markedly reduced exon 8.8 inclusion (Fig. 4C–E). This suggested that either Ud or Ds mutation, concomitant with a failure in inter-intronic RNA pairing, could specifically decrease exon 8.8 inclusion. Importantly, the reduced efficiency of exon 8.8 inclusion caused by individual mutations was recovered to WT level by compensatory double mutations (Fig. 4C–E). Interestingly, the results of mutational studies were consistent with the frequency analysis of the alternative exon 8 inclusion in naturally occurring deletion mutants (Fig. 4F and G). We found that 2 D. melanogaster mutant isolates lacked exon 8.7 and the selector sequence upstream of exon 8.8 (Mrp^ΔLeft8.8), which was consistent with the minimal inclusion of the last exon 8 variant in all tissues examined (Fig. 4F). In contrast, the last exon 8 variant generally occupied a considerable portion of all examined tissues of WT isolates. Taken together, both artificial mutation experiments and naturally existing mutant analysis indicated that the last exon 8 variant is dependent on upstream RNA pairing.

Conservation and divergence of RNA pairing in insects

To trace the origin and evolution of such dual RNA structures, we systematically analyzed a major taxonomic group of insects and crustaceans, which last shared a common ancestor approximately 420 million years ago.³⁶ We found clear evidence that the conserved docking site is located in the intron upstream of constitutive exon 9 in the coleopteran MRP1, which could pair with selector sequences upstream of its alternative exons 8 (Fig. S10). Notably, a similar structural architecture was found in the exon 8 clusters of water flea species (D. pulex and D. magna) (Fig. S11), suggesting that a dual RNA pairing system to direct mutually exclusive splicing was established in the Pancrustacea ancestor gene. We concluded that these RNA architectures and bidirectional RNA pairing interactions are evolutionarily conserved components required for splicing throughout insect or arthropod evolution.

However, upstream and downstream RNA pairings did not evolve symmetrically. We found that downstream docking site-dependent RNA pairing was mainly responsible for exon cluster expansion during insect evolution. However, we also found strong evidence for species-specific expansion of upstream docking site-dependent RNA pairing with exon 8 duplication in M. domestica (Fig. S8). Moreover, downstream docking sites were generally more conserved than upstream docking sites, particularly in Coleoptera and Daphnia species (Figs. S10 and S11). Differences in sequence conservation between the upstream and downstream docking sites are likely to reflect adaptable expansions of the guided exon cluster, with more conserved docking sites purifying selected by pairing with more selector sequences. Overall, evolutionary history and expansion processes for the bidirectional RNA pairing system shaped a regulatory network of species-specific mutually exclusive splicing.

In addition to RNA pairing interactions between docking site and selector sequences, the positions were also subject to purifying selection, particularly for the selector sequence (Fig. S12). For example, selector sequences are located ∼25 bp from the 5′ alternative splice site in the coleopteran species investigated, whereas their intron sizes differed markedly (Fig. S12). Similar patterns have been observed for Stem III in non-bee hymenopteran species. In contrast, the short distance between the selector element and the 3′ splice site was highly conserved among Lepidoptera (48–75 nt) and bee species (291–377 nt), whereas the 5′ splice site and the selector element were divergently large (Fig. S12). In the case of RNA pairings near the splice site, they mainly tended to act as an approximation activator for the most proximal exon, such as stems I that direct exon 8.2 in lepidopterans and bees (Figs. 2 and 3). When RNA pairings are divergently distant from the exon outside the loop, they have a primarily suppressive function for the exon within the loop, as shown by the suppression of stems II for exon 8.2 inclusion in lepidopterans and bees. Different evolutionary constraints of position may reflect the tuning and even switching between different functions of RNA pairing. This result also strongly suggests an unusual mode of regulating the selection of alternative exons by tuning the positions of these cis-elements.

RNA competing pairings are conserved in chordate MRPs

Finally, we examined whether other cells used a splicing mechanism of commitment to the one-only choice from a tandem array of exons that was similar to that used by arthropod MRP. Branchiostoma belcheri and B. floridae MRPs contained 7 and 10 duplicated exons, respectively, which were spliced in a mutually exclusive manner. We discovered similar arrangements of intronic elements and one-to-more RNA pairings within exon 8 clusters (Fig. 5A–C). Given the large size of the Branchiostoma MRP exon 8 cluster and the concomitant technical difficulties, we generated deletion constructs to test the effects of RNA pairing on splicing in Drosophila S2 cells (Fig. 5A and B).

Figure 5. The RNA-pairing-directed mechanism is conserved in Chordates. (A) Overview of partial genomic organization and minigene constructs of the Branchiostoma MRP1 gene. Constitutive exons (black boxes), alternative exons (blue boxes), conserved intronic elements IE1 and IE2 (blue), and introns (line) are indicated. Below are sequences of identical intronic elements (IE) in B. belcheri (Bbe) and B. floridae (Bfl) with spacings as indicated (nt). The dashed arrow represents RNA-RNA interactions. (B) Predicted RNA pairing of B. belcheri MRP pre-mRNA. Mutations introduced into the dsRNA stem are indicated on the left or right as mutated sequences (M1–M2). (C) Predicted mutually exclusive RNA pairings of B. floridae MRP pre-mRNA. (D) Effects on exon 8 inclusion of disruptive single mutations (M1 and M2) and compensatory double mutations (M12). (E) Effects of mutations on alternative exon 8 selection. Data are expressed as a percentage representing the mean ± SD of 3 independent experiments. *P < 0.05 (Student's t-test).

The reduced efficiency of exon 8.6 inclusion caused by individual mutations (M1, M2) was restored by compensatory double mutation (M12) in the B. belcheri minigene (Fig. 5D and E). These observations suggested that RNA pairings may control the mutually exclusive choice of exon variant in Branchiostoma MRP. Moreover, we also found similar arrangements and architectures of the cis elements of MRPs from the colonial chordate, B. schlosseri, and the annelid, C. teleta (Fig. S13), although this has not yet been confirmed experimentally due to the lack of an effective experimental system. Therefore, we concluded that the dsRNA-directed mechanism for directing mutually exclusive splicing is evolutionarily conserved in chordates and even in bilaterians.

Discussion

We demonstrate that MRP1 genes independently underwent tandem exon duplication in Nematoda, Platyhelminthes, Annelida, Mollusca, Arthropoda, Echinodermata, and early-diverging Chordata. Our findings also indicate the convergent evolution of exclusive spliced exons for several types of ion channels in arthropods and vertebrates.³⁷ MRPs are transmembrane proteins, and thus they exhibit the general adaptive advantages of mutually exclusive splicing suggested by Copley. However, MRP1 genes in late-diverging vertebrates do not contain tandemly duplicated exons. Similarly, no exon duplication has occurred in vertebrate Dscam genes while insect Dscam1 genes underwent multiple exon duplication events,³⁸ suggesting that there is an adaptive benefit to maintaining at least some of these genes during metazoan evolution. One possibility is that identical expression patterns of isoforms will be maintained, which is arguably less likely with paralogues. Another possibility is that, because the choice of exon is ultimately made at the post-transcriptional stage, the properties of the channel may be modified on a shorter timescale and in a more precise manner than would be possible by repressing transcription of one gene and initiating transcription of another. Intriguingly, all of the ion channel families discussed here undergo RNA editing in Drosophila,³⁷ providing another means of post-transcriptional modification of ion-channel functions, and again suggesting an adaptive benefit that could not be achieved with 2 gene copies.

Traditional models of regulated alternative splicing involve cis-regulatory elements (i.e. splicing enhancers and silencers) in addition to the canonical 5′ and 3′ splice site signals, which either enhance or repress the ability of the spliceosome to recognize particular splice sites through binding to trans-acting factors.^3,4 Mutually exclusive splicing is one of the most elaborately regulated forms of alternative splicing, which may require extra components. Recent studies have indicated that RNA secondary structures play crucial roles in the higher-order regulation of complex mutually exclusive splicing, particularly for more than 2 variable exons (e.g. Dscam, which has 38,016 isoforms).^18-20,23,24 Additionally, such splicing control systems generally include enhancers upstream of the docking site, silencers, weak splice sites, and their binding proteins.^39,40 The competing RNA pairing systems described here are involved in directing mutually exclusive splicing not only in Arthropoda but also in other phyla, such as Chordata and Annelida, suggesting that it may be an evolutionarily conserved mechanism in bilaterians.

RNA pairing systems evolved with the expansion of exon clusters under intense purifying selection. By integrating the genetic and molecular data from 121 metazoan species, we developed a credible long- and short-term phylogeny of mutually exclusive splicing, directed by an RNA pairing controller, in MRP1 (Fig. 6). Although the RNA pairings are an evolutionarily conserved mechanism in bilaterian species, these competing pairing systems originated independently in various phyla, suggesting the mechanistic convergence of mutually exclusive splicing. In arthropod MRP1, mutually exclusive splicing seems to have initially evolved in the common ancestor of the Pancrustacea and Myriapoda, possibly beginning with the unidirectional downstream pairing system. During subsequent evolution, the RNA structural system could have been made more complex in 2 ways. One, the simultaneous duplication of selector sequences with exon expansion through non-allelic homologous recombination,³⁵ would form one-to-many competitive pairing systems, such as those in coleopterans and P. puparum (Fig. S10; Fig. 3). Alternatively, the emergence of upstream RNA pairing would generate another competitive system between upstream and downstream RNA pairing, such as those in lepidopterans and hymenopterans (Figs. 2 and 3). During dipteran speciation, these 2 sets of pairing systems evolved in concert with exon cluster expansion (i.e., M. domestica and Anopheles gambiae, Figs. S8 and S9). These pairing systems could potentially evolve, through additional rounds of recombination, to flexibly direct mutually exclusive splicing for an ever-changing exon cluster during speciation.

Figure 6. RNA-pairing-guided evolution of species-specific mutually exclusive splicing. Extant MRP1 pre-mRNA structures and proposed ancestor molecules are shown associated with a cladogram indicating the phylogenetic relationships of organisms analyzed in this study. Intronic RNA pairings in D. melanogaster (Dme) were confirmed in Drosophila S2 cells (Fig. 4), and B. mori RNA pairing was confirmed in BmN cells (Fig. 2). Intronic RNA pairing in A. mellifera, P. puparum, and B. floridae was validated using disruptive and compensatory mutations in heterologous Drosophila S2 cells, respectively (Figs. 3 and 5). Because the minigene constructs for T. castaneum were not effective in transfection experiments with Drosophila S2 and BmN cells, we did not experimentally validate their proposed structures. Abbreviations used are listed in Table S1. Nodes denoting the independent ancestral origins of exon duplication are indicated by solid circles. Dashed lines indicate hypothesized relationships. Lost exons are shown in the dashed circles. Constitutive MRP1 mRNAs are depicted as unstructured (black).

Materials and methods

Materials

Insects and other organisms used in this study are listed in Table S1. Fruit flies (D. melanogaster), silkworms (Bombyx mori), red flour beetles (Triboleum castaneum), honeybees (A. mellifera) and other species were obtained as reported previously.²⁰ Natural isolates of D. melanogaster species were donated by the Kunming Institute of Zoology, Chinese Academy of Sciences and the College of Life Sciences, Hubei University, China. Buccinum belcheri samples were donated by Dr. Wang (School of Life Sciences, Xiamen University, China), and other insect species (i.e., P. puparum), were donated by the Institute of Insect Sciences, Zhejiang University.

Annotation and identification of MRP1s

The sequences of MRP1 genes in some arthropod species were obtained from FlyBase and NCBI (http://flybase.org/blast/ and http://blast.ncbi.nlm.nih.gov/blast, Table S1). The genomic sequences of the MRP1 genes for the other species were determined by polymerase chain reaction (PCR) and sequencing after obtaining cDNA sequences by reverse transcription PCR (RT-PCR). Total RNA was isolated using TRIzol Reagent (Life Technology, Carlsbad, CA) and was reverse transcribed using SuperScript III Reverse Transcriptase (Life Technology, Carlsbad, CA) with oligo(dT)₂₀ primer (Life Technology, Carlsbad, CA). Degenerate primers were used to amplify products from organisms within different clades (i.e., Drosophila and Lepidoptera). PCR was performed with initial denaturing at 95°C for 2 min, then 35 cycles of denaturing at 95°C for 30 s, annealing at 50–55°C for 30 s, and extension at 72°C for 1 min, followed by a final extension at 72°C for 10 min. PCR and RT-PCR products were cloned into the pGEM-T Easy vector (Promega, Madison, WI) for sequencing. Specific primers were designed based on cDNA sequences, and their sequences are listed in Table S2. Whole genomic DNA was extracted using the MiniBEST Universal Genomic DNA Extraction Kit (TaKaRa, Dalian, China). PCR cycling conditions were as follows: 95°C for 2 min, followed by 35 cycles of 95°C for 30 s, 55–60°C for 30 s, and 72°C for 2–4 min, with a final extension at 72°C for 10 min. Amplification products were cloned into the pGEM-T Easy vector for sequencing.

Sequence alignment and phylogenetic reconstruction

Sequences were aligned using the ClustalW2 program (http://www.ebi.ac.uk/clustalw/index.html).⁴¹ Mutually exclusive exons were translated and the resulting polypeptides were aligned. A maximum-likelihood tree was constructed from the amino acid sequences using MEGA 6.0 software.⁴²

Predictions of RNA pairing, protein 3D structures, and branch points

The RNA pairing interactions between the docking site and selector sequences were predicted using the Mfold program.⁴³ All protein 3D structures were acquired using Swiss-Model in automated mode (www.swissmodel.expacy.org).⁴⁴ The structures were displayed and processed using PyMOL (www.pymol.com), and MRP (isoform M, FlyBase ID: FBgn0032456) modeling was based on template 2hydB (pdb ID). Branch points were predicted using SVM-BP finder (http://regulatorygenomics.upf.edu/Software/SVM_BP/).⁴⁵

Quantification of mRNA splice isoforms

We assayed the ratio of RNA splice isoforms using RT-PCR followed by exon-specific restriction digestion. Total RNA was isolated from various developmental stages, tissue samples, and cells then reverse transcribed using SuperScript III Reverse Transcriptase with gene-specific primers. Specific primers flanking the alternative exons were designed. PCR cycling conditions were as follows: 95°C for 2 min, followed by 35 cycles of 95°C for 30 s, 55–60°C for 30 s, and 72°C for 30–60 s, with a final extension at 72°C for 10 min. PCR products were then digested by exon-specific restriction enzymes. Images were captured using a CCD camera and mRNA isoforms were quantified by comparing the integrated optical density of the detected bands using the GIS 1D Gel Image System (ver. 3.73; Tanon, Shanghai, China). The error bars displayed were calculated from the average of 3 independent experiments. Significance was determined using 2-tailed Student's t-tests with thresholds of *P < 0.05, **P < 0.01, and ***P < 0.001.

Minigene construction, mutagenesis, and transfection

PCR was performed to amplify a DNA segment encompassing the upstream constitutive exon, the downstream constitutive exon, and the intervening sequence for minigene construction. Site mutagenesis was performed on both elements to disrupt the RNA stem, as shown in the schematic diagrams of minigene construction. Compensatory double mutations were performed to restore the RNA stem structure, and all constructs were confirmed by sequencing. Drosophila and A. mellifera MRP1 WT and mutant constructs were inserted into the pMT/V5-HisB vector (Invitrogen, Carlsbad, CA) under control of the metallothionein promoter. Likewise, silkworm and P. xylostella WT and mutant constructs were cloned into the pBluescript-SK+ vector (Stratagene, Santa Clara, CA) under control of the Baculovirus IE1 promoter (obtained from C. Gong, Suzhou University, China). Plasmids were isolated using the Endo-free Plasmid Mini Kit I (OMEGA, Norcross, GA) and transfected into Drosophila S2 cell lines or silkworm cells (BmN) using Lipofectamine 2000 (Life Technology, Carlsbad, CA) according to the manufacturer's protocol. To induce expression of the minigene, copper sulfate was added directly to the S2 cell culture medium to a final concentration of 500 μM. Total RNA was harvested 48 h after transfection using TRIzol reagent. RT-PCR and exon-specific restriction enzyme analysis were performed to detect alternative splicing. The error bars presented in this study were calculated from the average of 3 independent experiments.

Accession numbers

The insect and Crustacea MRP1 gene sequences were deposited into GenBank with accession numbers JF501607-JF501637, KJ776428, KJ776430, KJ776433-KJ776435, KJ776439, KJ776440, KJ776443, KJ776447, KF776448.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Funding

This work was partly supported by research grants from the National Natural Science Foundation of China (31430050, 31630089, 31125011, 31270844).