Transcriptomic comparison of Drosophila snRNP biogenesis mutants reveals mutant-specific changes in pre-mRNA processing: implications for spinal muscular atrophy.

Survival motor neuron (SMN) functions in the assembly of spliceosomal small nuclear ribonucleoproteins (snRNPs) that catalyze pre-mRNA splicing. Here, we used disruptions in Smn and two additional snRNP biogenesis genes, Phax and Ars2, to classify RNA processing differences as snRNP-dependent or gene-specific in Drosophila Phax and Smn mutants exhibited comparable reductions in snRNAs, and comparison of their transcriptomes uncovered shared sets of RNA processing changes. In contrast, Ars2 mutants displayed only small decreases in snRNA levels, and RNA processing changes in these mutants were generally distinct from those identified in Phax and Smn animals. Instead, RNA processing changes in Ars2 mutants support the known interaction of Ars2 protein with the cap-binding complex, as splicing changes showed a clear bias toward the first intron. Bypassing disruptions in snRNP biogenesis, direct knockdown of spliceosomal proteins caused similar changes in the splicing of snRNP-dependent events. However, these snRNP-dependent events were largely unaltered in three Smn mutants expressing missense mutations that were originally identified in human spinal muscular atrophy (SMA) patients. Hence, findings here clarify the contributions of Phax, Smn, and Ars2 to snRNP biogenesis in Drosophila, and loss-of-function mutants for these proteins reveal differences that help disentangle cause and effect in SMA model flies.

INTRODUCTION

Survival motor neuron (SMN) protein functions in the assembly of Sm-class small nuclear ribonucleoproteins (snRNPs), core components of the spliceosome (Fischer et al. 2011; Matera and Wang 2014). Loss of SMN causes the neuromuscular disease spinal muscular atrophy (SMA) (Lefebvre et al. 1995), but the extent to which SMA is caused by disruptions in snRNP levels is not known. Loss of SMN causes transcriptome-wide changes in gene expression and RNA processing that are often assumed to stem from downstream disruptions in snRNP supply (Zhang et al. 2008, 2013; Garcia et al. 2013), but direct transcriptomic comparisons of SMA models with other targeted disruptions in snRNP biogenesis genes have, to date, not been performed.

The nuclear export of Sm-class small nuclear RNAs (snRNAs) is a key step in snRNP biogenesis. Sm-class snRNPs catalyze pre-mRNA splicing in the nucleus, but some steps of snRNP biogenesis occur in the cytoplasm, including assembly by SMN and associated SMN complex proteins (Matera and Wang 2014). The phosphorylated adaptor for RNA export protein (PHAX) mediates the nuclear export of newly transcribed snRNAs by linking them to the Exportin 1 (XPO1/CRM1) nuclear export machinery (Ohno et al. 2000). PHAX interacts with the 5′-terminal cap binding complex (CBC) and adjacent 5′-proximal RNA of RNA polymerase II transcripts, but it is displaced from longer mRNAs by the heterogeneous nuclear ribonucleoprotein C (C1/C2) tetramer (Ohno et al. 2000; McCloskey et al. 2012; Ohno 2012). Hence, PHAX functions to selectively export smaller snRNAs with lengths approximately 100 to 250 nucleotides (nt), making it an attractive target for genetically manipulating snRNP levels.

The nuclear export of Sm-snRNAs follows transcription termination. Prior to nuclear export, snRNAs are cotranscriptionally cleaved from downstream sequences. Unlike mRNAs, snRNAs are cleaved but not polyadenylated (Matera and Wang 2014). The Integrator is a multiprotein complex that cleaves snRNAs to form 3′-ends that lack a poly(A) tail (Baillat et al. 2005; Baillat and Wagner 2015). Proper 3′-end maturation of snRNAs is cotranscriptionally linked to upstream promoter sequences and 5′-ends (de Vegvar et al. 1986; Hernandez and Weiner 1986; Ezzeddine et al. 2011). Recently, the CBC-binding protein arsenite resistance protein 2 (ARS2/SRRT) was implicated in this crosstalk between snRNA 3′- and 5′-ends (Hallais et al. 2013). Knockdown of ARS2 or CBC proteins in mammalian cells increased readthrough of endogenous snRNAs and snRNA 3′-end reporters (Andersen et al. 2013; Hallais et al. 2013). In addition to snRNA 3′-end maturation, ARS2 contributes to the 3′-end maturation of replication-dependent histone mRNAs, at least in mammals, and it also has a conserved role in the regulation of microRNA biogenesis (Grigg et al. 2005; Lobbes et al. 2006; Yang et al. 2006; Gruber et al. 2009, 2012; Sabin et al. 2009). In this context, ARS2 is likely a less selective target than PHAX for disrupting the supply of snRNPs.

Observations in diverse eukaryotes suggest that perturbations in snRNP supply produce transcriptome-wide effects, particularly in alternative splicing and polyadenylation. Mutations in yeast genes encoding core spliceosomal factors disrupted the splicing of a subset of pre-mRNAs, including splicing events in highly expressed pre-mRNAs for ribosomal protein genes (RPGs) (Clark et al. 2002; Pleiss et al. 2007). Inversely, repressing the transcription of these highly expressed RPGs changed the alternative splicing patterns of other yeast pre-mRNAs, suggesting that spliceosomal accessibility controls pre-mRNA splicing patterns (Munding et al. 2013). In Drosophila, core splicing factors were uncovered by an RNAi screen for regulators of a predetermined set of alternative-splicing patterns (Park et al. 2004). In cultured human cells, knockdown of Sm-class snRNP-associated proteins led to numerous changes in alternative splicing (Saltzman et al. 2011). Furthermore, blocking U1 snRNP specifically with an antisense morpholino led to premature cleavage and polyadenylation of mRNAs in cells from diverse lineages (Kaida et al. 2010; Berg et al. 2012). These studies highlight that multiple strategies can be used to effectively disrupt snRNP supply and function, which creates an opportunity to use comparative transcriptomics to comprehensively categorize the snRNP dependence and/or spliceosomal accessibility requirements of disparate RNA processing events.

Here, we carried out transcriptome profiling of Drosophila Smn mutants, along with two additional snRNP biogenesis mutants in order to identify snRNP-dependent versus gene-specific changes in gene expression, premature cleavage and polyadenylation, and alternative splicing. Specifically, transposon insertion mutations in the fly orthologs of Phax and Ars2 were used for comparison with Smn null mutants (Ohno et al. 2000; Kitao et al. 2008; Laubinger et al. 2008; Gruber et al. 2009; Andersen et al. 2013; Hallais et al. 2013; Pabis et al. 2013). The Phax mutant was particularly useful in identifying an overlapping set of snRNP-dependent RNA changes in SMA model flies, and the analysis of the Ars2 mutants revealed a distinct set of RNA processing events that likely depend on the integrity of the CBC and 5′-cap. Additionally, alternative-splicing differences shared between the Phax and Smn mutants could be rescued at both the RNA and protein levels by expression of Phax or Smn transgenes, respectively. In combination with a modest restoration of steady-state snRNA levels, observations from these rescue lines corroborate the connection between snRNP levels and specific pre-mRNA splicing events. However, the link between these changes and SMA phenotypes remains unclear.

RESULTS

Phax and Smn mutants display similar steady-state snRNA decreases

RNA-sequencing (RNA-seq) of poly(A)-selected RNA was used to identify mRNA changes in snRNP biogenesis mutants, and Northern analysis coupled with RNA-seq of total ribosomal RNA-subtracted (rRNA−) RNA was used to quantify steady-state snRNA levels. Browser shots of Phax, Smn, and Ars2 gene loci confirmed the gene disruptions in the respective mutant lines (Fig. 1A; Garcia et al. 2013). In our previous RNA-seq analysis of Smn null mutants, publicly available modENCODE data were used to determine the developmental stage at which the mutants arrested (Garcia et al. 2013). By the same method, the Phax and Ars2 mutants exhibited a similar timing of developmental delay as did the Smn animals (Supplemental Fig. S1A,B). Using developmentally appropriate wild-type controls (Supplemental Fig. S2A), Northern blots demonstrated comparable, near twofold decreases in steady-state U1, U2, and U4 snRNAs in the Phax and Smn mutants (Fig. 1B,C). When assayed by Northern blotting, U5 snRNA levels were essentially unaffected in the Phax mutants (Fig. 1B,C), as were levels of nearly all of the snRNAs in the Ars2 mutants (Fig. 1B,C). When RNA-seq data were used to quantify snRNA levels, each of the three mutants appeared to exhibit snRNA deficits, but the extent varied by the method of normalization (Supplemental Fig. S2B,C). Only small (<15%) changes were observed in specific snRNA paralogs relative to the total (Supplemental Fig. S3), quantified by using reads covering isoform-specific nucleotides (Lu and Matera 2014). Overall, the Phax and Smn mutants displayed the most consistent decreases in steady-state snRNA levels.

FIGURE 1.

Analysis of three different snRNP biogenesis mutants. (A) Browser shots of larval RNA-seq reads at the three disrupted gene loci. Reads from Oregon-R (WT) versus Phax, Smn, or Ars2 snRNP biogenesis mutants are shown. Mapped read tracks were normalized to the median of the middle two quartiles of the mapped WT sequence read counts. (B) Northern blot of RNA from snRNP biogenesis mutants and WT larvae. RNA was extracted from WT larvae at 60 ± 2 h (60 h) post egg-laying, and mutant RNA was extracted at 74 ± 2 h to account for their delayed development (see text). (C) Quantification of Northern blot in B. Sm-class snRNAs were normalized to the Lsm-class U6 snRNA, and U6-normalized WT snRNA levels were set at 100. Asterisks are P-values from a Student's t-test: (*) P-value ≤0.05, (**) P-value ≤0.01, (***) P-value ≤0.001.

snRNP biogenesis mutants share a small number of gene expression changes

To identify shared and distinct changes in gene expression between snRNP biogenesis mutants, we used the TopHat/Cufflinks pipeline to quantify mRNA differences from our poly(A)-selected RNA-seq data (Trapnell et al. 2012). Raw reads from age-matched wild-type (WT) controls and modENCODE second (L2) and early third instar (L3 12h) developmental controls were integrated with sequencing data from the mutants into a single analytical pipeline. For comparison, we focused on genes whose expression changed uniformly and significantly (false discovery rate-adjusted P-value <0.05) relative to each of the controls (L2, L3 12h, and our WT). A relatively small number of genes and undifferentiated loci with multiple genes (186 out of ∼15,000) met these stringent criteria for differential expression (i.e., 186 loci were significantly different from all of the controls; total number from Supplemental Table S1), and most of these differences were mutant specific (Fig. 2A,C; Supplemental Tables S1–S4). A number of stress-responsive transcripts were among those that were specifically elevated in Smn mutants, though some of these were also up slightly in the Ars2 mutants (Fig. 2B). Phax mutants shared only 28 gene expression changes with Smn mutants (Fig. 2C; Supplemental Table S2). This small list of genes was enriched for genes involved in oxidation–reduction, as measured by DAVID gene ontology term analysis (P-value = 0.003) (Huang et al. 2007, 2009a,b). A notable gene expression change shared by the Phax and Smn mutants was a fivefold increase in the Activity-regulated cytoskeleton associated protein 1 (Arc1) transcript (Supplemental Table S2). Arc1 functions in a well conserved hyperlocomotion response to starvation, and hence elevated levels in Phax and Smn mutants might be an indirect consequence of a disruption in nutrient acquisition in these animals (Mattaliano et al. 2007). In summary, gene expression analysis uncovered a small number of overlapping changes that exceeded normal developmental fluctuations. The relationship of these changes to SMA etiology and/or pathology remains to be determined.

FIGURE 2.

Gene expression differences in snRNP biogenesis mutants. (A) Heatmap comparison of Cuffdiff FPKM levels of differentially expressed transcripts. Heatmap colors were rescaled for each row and the rows were clustered based on pattern of gene expression between WT and mutants. (B) A heatmap of FPKMs from a set of stress-responsive transcripts in A. (C) Venn diagram of the overlap in gene expression differences in the snRNP biogenesis mutants. Phax = Phax, Smn = Smn, and Ars2 = Ars2. Numbers in parentheses are totals.

snRNP biogenesis mutants exhibit a trend toward shorter mRNAs

In addition to its role in splicing, U1 snRNP has been reported to control mRNA length by preventing premature cleavage and polyadenylation (Kaida et al. 2010; Berg et al. 2012). To determine whether reduced U1 snRNA levels in our snRNP biogenesis mutants correlated with expression of shorter mRNAs, mapped RNA-seq reads were analyzed using the DaPars linear regression algorithm (Masamha et al. 2014; Xia et al. 2014). DaPars output is displayed as a percentage of distal poly(A) site usage index (PDUI) for each transcript, or as the difference in the PDUI value between wild-type and mutant (ΔPDUI). Pairwise comparisons of PDUI values for wild-type and mutants can be used to visualize changes in RNA length. Whereas the Phax and Smn mutants expressed transcripts that were both longer and shorter than their wild-type counterparts, there were greater numbers of shortened mRNAs (Fig. 3A–C; Supplemental Tables S5–S7). Focusing on differences that are greater than twofold, Phax and Smn mutants had a relatively small number of mRNA length changes in common with each other (34) (Supplemental Table S5). However, the changes were largely distinct from those observed in Ars2 animals, who had more than twice as many mRNA length changes as did the Phax and Smn mutants combined (Fig. 3D; Supplemental Table S8). The individual direction of length change, either longer or shorter, was also factored into the consideration of overlapping differences in Figure 3D. The overall trend of shorter transcripts in the mutants does not appear to be a consequence of developmental delay, as parallel DaPars analyses of modENCODE data sets of second to third instar larvae showed a trend toward shorter transcripts in the more mature animals (Supplemental Fig. S4). Among the few overlapping mRNA length changes between Phax and Smn, one of the largest and most apparent changes was a shortened mRNA for CG9662, a putative oligosaccharyl transferase complex subunit (ΔPDUI = −0.96, Supplemental Table S5). These results are generally consistent with the notion that reductions in U1 snRNP levels predispose mRNAs to premature cleavage and polyadenylation. Because the Ars2 mutants exhibited only a small decrease in steady state U1 levels (<10%; Fig. 1C), the mRNA length changes in these animals are more likely due to disruptions of the CBC, which are known to affect a subset of pre-mRNA termination events (Wong et al. 2007; Andersen et al. 2013; Hallais et al. 2013).

FIGURE 3.

mRNA length changes in snRNP biogenesis mutants. Pairwise comparison of significant (false discovery rate adjusted P-value <0.05) differences in mRNA length between WT and Phax (A); WT and Smn (B); and WT and Ars2 (C). RNA length differences were extracted from mapped RNA-seq read files by the DaPars linear regression algorithm (see text). Percentage of distal poly(A) site usage index = PDUI. (D) Venn diagram of RNA length differences in snRNP biogenesis mutants labeled as in Figure 2.

Phax and Smn mutants share a set of alternative splicing pattern changes

Our previous RNA-seq analysis of Smn null animals revealed many small-amplitude alternative splicing differences (compared to wild-type controls), but it was unclear whether any of these changes might be due to decreases in snRNP levels (Garcia et al. 2013). To ascertain whether any of these splicing changes are snRNP-dependent, we used the mixture of isoforms (MISO) probabilistic framework to compare splicing differences in Smn mutants with those identified in Phax and Ars2 mutants (Katz et al. 2010). Simply put, the MISO algorithm estimates the expression of alternatively spliced isoforms in RNA-seq data, and reports it as a fraction of alternatively spliced events relative to the total (Katz et al. 2010). MISO denotes this fraction as the “percent spliced in” or “Psi” (Katz et al. 2010). Differences are reported here as an absolute change relative to the wild-type control, |deltaPsi|. Similar to our previous results with the Smn null animals, we found mostly small changes (less than twofold, |deltaPsi| <0.5) in alternative-splicing events in the Phax and Ars2 mutants (Fig. 4A–C; Supplemental Tables S9–S11). More than half of the splicing changes identified in the Smn null animals overlapped with changes in the Phax mutants, but neither the Phax nor Smn mutants had considerable overlap with those observed in the Ars2 mutants (Fig. 4D; Supplemental Table S11). Of the top 150 differences in the Ars2 background, close to 90% of these changes were located in the 5′-proximal intron (Supplemental Table S11). These changes likely reflect the association of Ars2 with the CBC and the role of this complex in promoter proximal splicing events (Laubinger et al. 2008) rather than Ars2-related changes in snRNPs. In contrast, the changes that were common to Phax and Smn mutants are more likely to be due to decreases in snRNP levels.

FIGURE 4.

Alternative-splicing differences among snRNP biogenesis mutants. Distribution of significant alternative-splicing changes from small |delta Psi| near zero to large |delta Psi| near one: (A) Phax mutants relative to WT; (B) Smn mutants relative to WT; and (C) Ars2 mutants relative to WT. (D) Venn diagram of overlapping alternative-splicing changes in snRNP biogenesis mutants, labeled as in Figure 2.

Transgenic rescue of alternative splicing changes induced by loss of Phax or Smn

As described previously, Smn is essential for viability in the fly (Chan et al. 2003), and transgenic expression of Flag- or GFP-tagged SMN protein in an otherwise null mutant background fully rescues the larval lethality and other SMA-like phenotypes (Rajendra et al. 2007; Shpargel et al. 2009). Here, we found that Phax is also an essential gene, and that ectopic expression of Phax-GFP using the GAL4-UAS system rescued the lethality of Phax mutants (Supplemental Table S12). We next assayed (and confirmed by semi-quantitative RT-PCR) a number of the most apparent alternative splicing changes identified by RNA-seq (Supplemental Fig. S5). We then tested the ability of the Phax and Smn transgenes to rescue four of these splicing events by quantitative qRT-PCR (Fig. 5B,C). Ubiquitous expression of Phax using an Armadillo-GAL4 driver significantly rescued all of the splicing changes we tested (Fig. 5B). Transgenic expression of Smn, driven by its native promoter, also rescued these alternative-splicing events in Smn null animals (Fig. 5C). Together, the results confirm that the observed changes we identified by RNA-seq are indeed caused by loss of Phax or Smn expression.

FIGURE 5.

Transgenic rescue of alternative-splicing changes in Phax and Smn mutants. (A) UCSC browser shots of mapped read tracks for four of the top 15 identified alternative splicing changes. (B) qRT-PCR analysis of intron retention in Phax mutants versus Phax mutants expressing a UAS-Phax transgene from the armadillo promoter driven GAL4 (Arm > Phax). (C) qRT-PCR of intron retention in Smn mutants versus Smn mutants expressing a wild-type Smn transgene from its native promoter. For B and C, levels of WT intron retention were set at one. (D) Western blot for dUTPase. Alternate protein isoforms reflect alternative-splicing change in dUTPase mRNA. Wild-type Smn transgene in Smn background = Tg:Smn. Long dUTPase protein isoform = PA and short = PB.

One of the splicing changes we detected is expected to produce an alternative (larger) dUTPase protein isoform. Therefore, we assayed the impact of this alternative splicing change by Western blotting, using an antibody that targets both dUTPase isoforms. As shown in Figure 5D, the expected increase in the ratio of the long isoform relative to the short one was observed in the mutants, and this change was largely rescued in both transgenic lines. Thus, the transcriptomic changes we identified in SMA model animals can (and do) translate into detectable differences at the protein level.

An increase in steady-state snRNA levels correlates with rescue of alternative splicing

If the alternative splicing changes, observed in both Phax and Smn mutants, are truly caused by the decreased availability of snRNPs, then rescue of steady-state snRNA levels should mirror the rescue of these alternative splicing changes in the Phax and Smn transgenic lines. Northern blotting was used to determine whether the expression of Phax and Smn transgenes increased steady-state snRNA levels in the mutants (Supplemental Fig. S6A,B). Expression of the Phax transgene by the Armadillo driver significantly rescued U1 and U4 snRNA levels in the Phax mutant background (Supplemental Fig S6A,B). All of the snRNAs tested here were significantly rescued by expression of the wild-type Smn transgene in the Smn null background (Supplemental Fig. S6C,D). For the wild-type Smn rescue line, these observations differ from our previously published results (Praveen et al. 2012), wherein we detected a modest increase in snRNA levels in the rescue line. This difference likely reflects minor improvements to our Northern blotting protocol, specifically: using less total RNA per gel lane and more labeled probe in the hybridization reaction to ensure probe excess. The increases in steady-state snRNA levels that we observed in the Phax and Smn rescue lines correlated well with the phenotypic rescue at both the organismal (viability) and molecular (alternative splicing) levels.

Splicing changes shared between Phax and Smn mutants are also caused by knockdown of splicing factors Prp6 and Prp8

Although the correlation of snRNA levels with particular splicing events provides evidence that these events are acutely sensitive to disruption of the core spliceosomal machinery, it does not rule out other possibilities. For example, the covariance could be due to a direct interaction of splicing factors with Phax and SMN, or a result of indirect factors, such as altered developmental progression. In addition, snRNPs are well known to be stable RNPs, and we do not know at what point reductions in ongoing snRNP supply begin to impair spliceosomal function. To test whether the splicing changes in the Phax and Smn mutants are indeed linked to a reduction in snRNP supply and a corresponding disruption in spliceosome function, we depleted the tri-snRNP and U5-associated proteins Prp6 and Prp8 using RNA interference. As shown in Figure 6A–C, SMN depletion caused a slight (less than twofold) perturbation in the splicing events tested, but the Prp6 and Prp8 knockdowns disrupted the splicing of these transcripts to a much greater extent. Knockdown of Prp8 disrupted the splicing of the events identified above (Abi, αCOP, dUTPase, and Rpn2 intron 1 [int. 1]), but it did not significantly disrupt the splicing of a downstream intron in Rpn2, serving here as a control for broader spliceosomal dysfunction. The Prp6 and Prp8 knockdowns also altered dUTPase protein expression patterns to a greater extent than did knockdown of Smn (Fig. 6B). Steady-state snRNA levels were not significantly affected by any of these knockdowns, with the possible exception of U5 snRNA, which was slightly reduced (∼40%) in the Prp8 depletion (Fig. 6D,E). The lack of snRNA decreases in the SMN depleted samples is consistent with the well-known perdurance of snRNPs in cultured cells, and it likely mirrors the observed low-amplitude disruptions to splicing. Prp6 and Prp8 are core spliceosomal proteins, and their depletion is expected to severely compromise the spliceosome (Grainger and Beggs 2005; Will and Luhrmann 2011). Thus, the splicing changes seen in these knockdowns are perhaps not unexpected, but they nevertheless support the conclusion that the overlapping changes identified in the Phax and Smn null mutants are snRNP dependent.

FIGURE 6.

snRNP-specific protein knockdown in S2 cells. (A) qRT-PCR of untreated cells versus those treated with dsRNA for Smn, Prp6, or Prp8 mRNAs. (B) Western blot of dUTPase levels in RNAi-treated cells. Anti-SMN and anti-α-tubulin verify SMN knockdown and load, respectively. (C) qRT-PCR verification of Prp6 and Prp8 mRNA knockdowns. Levels in untreated cells were set at 100. (D) Northern blot of snRNA levels from S2 cell knockdowns, quantified in E. U6-normalized snRNA levels of WT were set at 100.

Hypomorphic Smn missense mutants rescue snRNP levels and do not display defects in alternative splicing

As mentioned, the snRNP biogenesis mutants described above are severe loss-of-function alleles. SMA is, however, a hypomorphic condition, as complete loss of SMN activity is lethal in all organisms tested. To better model the disease in the fly, we previously engineered an allelic series of transgenic animals that express Smn missense mutations known to cause SMA in humans (Praveen et al. 2012, 2014). In addition to the wild-type rescue transgene (Smn), we analyzed three hypomorphic lines that display intermediate (SMA type II) phenotypes: Smn, Smn, and Smn (Praveen et al. 2014). These lines cover a range of phenotypic outcomes: Smn is semi-viable (∼30% eclose as adults), Smn is pharate-lethal (very few animals eclose), and Smn mutants all die as pupae (Praveen et al. 2014). Importantly, in this system, the wild-type and mutant transgenes are integrated at the identical chromosomal locus and are expressed from the native promoter in an otherwise Smn null background (Praveen et al. 2012, 2014). Using qRT-PCR, we analyzed the alternative splicing patterns of the four transcripts described above. As shown in Figure 7A, none of the Smn missense lines exhibited a greater than twofold change (in either direction) relative to Smn. Consistent with this finding, the missense mutants also displayed little-to-no difference in steady-state snRNA levels, compared to the wild-type rescue line (Fig. 7B,C). To ensure that differences in developmental timing did not confound the analysis, animals were harvested after they began to display the wandering behavior that is characteristic of the late third instar stage of Drosophila development. We also measured snRNA levels of animals at slightly earlier and later developmental time points, and we found little difference between the mutant and the control lines (Supplemental Fig. S7). Steady-state snRNA differences in earlier animals may reflect subtle differences in developmental progression from wild-type controls, but these decreased levels were largely absent from the later pupal stage, which only exhibited significantly lower U4 levels in the single Smn line.

FIGURE 7.

Analysis of steady-state snRNA levels and alternative splicing of target genes in SMN missense mutants. (A) qRT-PCR analysis of intron retention of flies expressing SMA patient-derived missense mutations. Transgenic animals of the following generalized genotype were used: Smn, Flag-Smn, where Tg represents a WT, V72G, Y107C, or T205I transgene. Intron retention in the Smn transgenic rescue line was set at one. (B) Northern blots of snRNA levels in the missense mutant lines. (C) Quantification of snRNA levels from panel B. RNA levels in the Smn (WT) transgenic line were set at 100.

Smn missense mutants display limited overlap in splicing pattern changes

The large-amplitude changes identified in the severe loss-of-function mutants were not evident in the Smn hypomorphs; however, these lines could exhibit other splicing pattern or gene expression changes that might account for differences in their phenotypic outcomes. To identify RNA changes in the hypomorphic lines, we sequenced total (rRNA-subtracted) RNA from whole animals isolated ∼12 h after puparium formation, a developmental stage that Smn null animals never reach. The hypomorphs displayed numerous changes in their splicing patterns relative to the wild-type rescue line (1662 combined; Supplemental Tables S13–16). However, the missense mutants displayed few large amplitude (greater than twofold) changes (107/1662), and hardly any of these events were shared between them (42/1662) (Fig. 8A; Supplemental Table S13). Of the forty-two overlapping splicing changes, only six were previously identified as being snRNP-dependent (four) or Smn gene-specific (two) by comparison to data from the severe Phax and Smn loss-of-function lines (Fig. 8B). Browser shots of representative examples are shown in Figure 7C. We used qRT-PCR to validate one of the two Smn gene-specific changes. The scarface transcript (Supplemental Fig. S8; Supplemental Table S13) exhibited an increase in intron retention between the third and fourth exons in the Smn transgenic line (deltaPsi = −0.61 and a >100-fold change by qRT-PCR). In the fly, Scarface functions as a negative regulator of JNK (c-Jun N-terminal kinase) stress signaling (Rousset et al. 2010). As measured by RNA-seq, the Smn (deltaPsi = −0.23) and Smn (deltaPsi = −0.07) mutants exhibited smaller changes; these events were below the limit of detection by qRT-PCR (Supplemental Fig. S8A; Supplemental Table S13). Consistent with the notion that this is not a snRNP-dependent change, an increase in scarface intron retention was not observed following knockdown of Prp8 (Supplemental Fig. S8B). Thus, RNA-seq profiling of hypomorphic, SMA-causing missense mutations is in agreement with the qRT-PCR data described above, and together these findings demonstrate that the splicing changes identified in the Smn null animals are not conserved among the intermediate SMA models.

FIGURE 8.

Alternative splicing changes in Smn missense mutants. (A) Venn diagram of overlapping alternative-splicing changes in Smn missense mutants. (B) Venn diagram of the overlap in shared splicing differences in the Smn missense mutants with either snRNP-dependent changes or Smn gene-specific changes categorized by comparison of the Phax and Smn total rRNA(−) RNA-seq. (C) Browser shot examples of an Smn gene-specific change in alternative-splicing in scarface (top) or a snRNP-dependent change in CG43394 within a PAPLA1 intron on the opposite strand (bottom). Genotypes are labeled by Smn, which was expressed in the Smn null background.

Differential expression analysis reveals an activation of stress signaling in Smn missense mutant animals

In addition to changes in splicing patterns, our genome-wide analysis of RNA from the severe loss-of-function mutants revealed potential snRNP-dependent and Smn gene-specific changes in steady-state mRNA levels. Differential expression analysis of the Smn missense mutant lines revealed numerous changes in mRNA abundance in the most severe line, Smn (599), but far fewer in the other less severe lines, Smn (88) and Smn (59) (Fig. 9A,B; Supplemental Table S17–S21). Despite the low number of changes in the least severe Smn line, more than half of these overlapped with the other hypomorphic mutants (Fig. 9B). Surprisingly, many of the stress-responsive transcripts identified in the severe Smn null mutants were also among the mRNAs that were up-regulated in all three Smn transgenic lines, relative to the wild-type rescue line (Fig. 9C; Supplemental Table S18). As corresponding increases in mRNA levels were not found in the Phax mutants (Fig. 2B), these events were classified as Smn gene-specific. Thus, unlike the snRNP-dependent splicing changes (which are not conserved in the Smn hypomorphs), the intermediate models of SMA exhibited clear RNA signatures of stress that correlate with disease severity (Fig. 9C) and are independent of decreases in snRNP supply (Fig. 2B).

FIGURE 9.

Differential gene expression in Smn hypomorphic animals. (A) Heatmap comparison of Cuffdiff FPKM levels of differentially expressed transcripts. Heatmap colors were rescaled for each row and the rows were clustered based on pattern of gene expression between Smn rescue line and the Smn missense mutants. (B) Venn diagram of overlapping differences in mRNA levels. (C) A heatmap of FPKMs from a set of stress-responsive transcripts. Genotypes are labeled by Smn, which was expressed in the Smn null background.

DISCUSSION

Connecting SMA pathology to SMN's function in snRNP assembly and downstream expression of specific mRNA isoforms remains a necessary benchmark for defining how splicing changes might contribute to disease etiology. Therefore, a first step toward uncovering the molecular mechanisms that elicit SMA outcomes would be to separate snRNP-dependent phenotypes from those that are SMN specific. Transcriptomic profiling of Smn and Phax mutants enabled us to identify a robust set of snRNP-dependent alternative splicing changes in SMA model flies. We show that reductions in steady-state snRNA levels correlate with the appearance of these alternative-splicing events, and that rescue of snRNA levels in transgenic animals correlates with the restoration of normal splicing. Furthermore, depletion of core splicing factors (Prp8 or Prp6) essentially bypasses the snRNP depletion step, causing many of the same alternative-splicing changes. Together, these observations support a close connection between snRNP levels and proper gene expression in vivo. Whether or not disruptions to this ubiquitously required process contribute to SMA phenotypes is less clear.

Ars2, 5′ exons, and U1 telescripting

Changes in splicing and mRNA length in Phax and Smn animals were largely distinct from the more abundant changes in the Ars2 mutants (Figs. 3D, 4D). Ars2 is a component of the 5′ cap binding complex (CBC), and this association may underlie the transcriptomic changes that we identified in this mutant (Gruber et al. 2009; Gonatopoulos Pournatzis and Cowling 2014). The Ars2 splicing changes were mostly in 5′-proximal introns, and this finding is consistent with CBC mutants in yeast, Xenopus, and plants (Gonatopoulos Pournatzis and Cowling 2014). The association of Ars2 with the CBC likely also underlies the numerous mRNA shortening events we observed. These findings are consistent with previous studies in yeast and human cells that suggest the CBC protects a subset of weak poly(A) signals from premature cleavage and termination (Wong et al. 2007; Andersen et al. 2013; Hallais et al. 2013). Importantly, both types of Ars2 dependent mRNA processing events were largely distinct from the snRNP-dependent changes observed in the Phax and Smn mutants.

The distinct mRNA processing changes in Ars2 mutants may also reflect species-specific differences in the role of the CBC and Ars2. Unlike ARS2 knockdown experiments in mammalian cells (Gruber et al. 2012; Hallais et al. 2013), we did not observe defects in snRNA or histone mRNA processing. This finding is consistent with the observation that Ars2 protein is dispensable for proper histone mRNA maturation in Drosophila cells (Sabath et al. 2013), and it suggests that aspects of Ars2 function are not conserved between vertebrates and invertebrates.

Our identification of a large number of mRNA shortening events that are shared between Phax and Smn mutants supports a recently discovered role for the U1 snRNP in the protection of pre-mRNAs from premature cleavage and polyadenylation (PCPA). PCPA was previously observed by subtractive hybridization PCR, coupled with high-throughput sequencing of RNA from cells treated with antisense oligomers targeting the 5′-end of U1 snRNA (Kaida et al. 2010; Berg et al. 2012). This protective role for U1 snRNPs was termed “telescripting” (Kaida et al. 2010; Berg et al. 2012) and our study is the first to validate this process outside of cell culture models. Although we observed an overall trend toward shorter transcripts in our data sets (Fig. 3), we found only a few large-magnitude (greater than twofold) changes (Supplemental Table S2). This is perhaps unsurprising, given that treatment of cultured cells with relatively high concentrations of anti-U1 oligomers provides an acute and nearly complete inhibition of U1 snRNA function, whereas the mutations in Phax and Smn do not perturb U1 to the same extent in the mutant animals.

RNA signatures of disease

Alternative splicing depends on the regulated accessibility of pre-mRNAs to the spliceosome. Studies in a variety of eukaryotes have shown that changing this accessibility, either directly or indirectly, alters the splicing patterns of specific pre-mRNAs (Clark et al. 2002; Park et al. 2004; Pleiss et al. 2007; Saltzman et al. 2011; Munding et al. 2013). The extent to which snRNP levels contribute to normal splicing regulation and/or disease pathology is still being worked out; however, changes in the pattern of alternative-splicing have been reported in various SMA model systems (Zhang et al. 2008, 2013; Bäumer et al. 2009; Lotti et al. 2012). The splicing changes reported in these SMA models are typically assumed to stem directly from loss of SMN and snRNPs, but such changes are often difficult to distinguish from broad, indirect effects of disease pathology or developmental delay (Bäumer et al. 2009; Garcia et al. 2013). By identifying overlapping splicing changes in Phax and Smn mutants, we are now able to distinguish between snRNP-dependent changes and gene-specific effects.

Many of the alternative splicing changes shared by Phax and Smn occur in genes of known importance to neuromuscular biology. For example, the mammalian ortholog of the αCOP protein has a demonstrated interaction with SMN, and it colocalizes with SMN in neuronal projections (Peter et al. 2011; Custer et al. 2013). This interaction was recently shown to be required for neurite outgrowth in mammalian cell and zebrafish models of SMA (Li et al. 2015). Orthologs of the Abi protein appear to have important roles in neurogenesis via the regulation of actin dynamics (Proepper et al. 2007; Lin et al. 2009; Liebau et al. 2011). Dominant mutations in GARS cause the distal neuropathy Charcot–Marie–Tooth disease type 2D, and recent observations suggest that GARS may colocalize with SMN in human cells (Motley et al. 2010; Drew et al. 2011; Grice et al. 2015). Mouse models of SMA exhibit disruptions to ubiquitin homeostasis, and evidence of SMN interactions with proteasomal components, including Psmd1, the mouse ortholog of Drosophila Rpn2 (Aghamaleky Sarvestany et al. 2014; Wishart et al. 2014).

Despite our finding of splicing disruptions in genes of known importance to SMA and neuromuscular biology in severe loss-of-function mutants, these splicing changes were absent from less severe models of SMA. Splicing patterns are not conserved from flies to humans, and we would not expect to find similar splicing disruptions in mammalian models of SMA. Finding these genes at the top of our data set may merely reflect the sensitive nature of the splicing of genes that are important for neuromuscular biology, and the fact that transcriptome diversity is greatest in neuronal tissues, including the brain (Mohr and Hartmann 2014).

Although we found evidence for both snRNP-dependent and SMN-specific RNA changes in SMA models, it remains to be determined which of these, if any, contributes to disease pathology. Our analysis of severe and intermediate SMA models suggests that snRNP-independent, Smn gene-specific RNA changes may play a bigger role in SMA pathology than previously envisioned. The SMN-specific activation of stress signaling pathways is conserved in SMA models (Wu et al. 2011; Genabai et al. 2015) and correlates well with disease severity in Smn hypomorphs (Fig. 8C). This activation of stress signaling appears to be a common signature of neuromuscular disease (Deshpande and Rodal 2016).

Based on our finding of a potential activation in stress signaling in fly models of SMA, an important question moving forward will be to understand how loss of SMN causes this stress response. By text mining with the stress signatures identified in the SMA models described here, we found that mutations in the Activating transcription factor 3 (Atf3) gene cause a similar dysregulation in immune and metabolic homeostasis (Rynes et al. 2012). Among other things, this signaling pathway appears to regulate the cytosolic formation of SMN-containing U-bodies in human cells in response to bacterial pathogens (Tsalikis et al. 2015). If SMN does indeed intersect with a normal host cell response to bacteria, this would provide a parsimonious model for the activation of stress signaling in SMA models, and another avenue to explore potential connections to neuromuscular disease.

In conclusion, a direct comparison of snRNP biogenesis mutants revealed that defects in snRNP supply cause complex changes to the transcriptome. However, an analysis of three different Smn hypomorphic lines showed that RNA processing changes are unlikely to be the primary drivers of SMA pathophysiology. Our studies support a model of SMA etiology wherein general disruptions at the mRNA level and specific disruptions at the SMN protein–protein level are responsible for SMA phenotypes. In future experiments, the key will be to uncover the tissue specific molecular activities of SMN that are important for neuromuscular function. Ongoing work with these and other intermediate models of SMA should help clarify the relative contributions of the aforementioned factors to SMA pathology and etiology.

MATERIALS AND METHODS

Fly mutants and husbandry

Stocks were cultured on molasses and agar at room temperature (25 ± 1°C) in half-pint bottles. Oregon-R was used as the wild-type allele. The Smn microdeletion allele was a gift from S. Artavanis-Tsakonis (Harvard University, Cambridge, USA) (Chang et al. 2008). Transposon insertion alleles Smn (f01109), Ars2 (c00735), and Ars2 (e01128) were obtained from the Exelixis collection at Harvard Medical School (Rajendra et al. 2007; Sabin et al. 2009). The transposon insertion allele Phax (P{lacW}PhaxSH0641) was a gift from S. Hou (Oh et al. 2003).

The wild-type and mutant Smn transgenes are expressed using the native Smn promoter. Transgenic constructs were integrated into the PhiC31 landing site located at 86F8, which was previously recombined into an Smn background (Praveen et al. 2012, 2014). For RNA-seq analysis, Smn transgenes in the Smn null background were crossed to Smn animals, and trans-heterozygous animals were put into TRIzol <12 h after pupation, as brown prepupa. For qRT-PCR and developmental northerns, Smn transgenes in the Smn background were self-crossed, and homozygous animals were put into TRIzol at 72–77 h post egg laying or <12 h after pupation before air bubble migration.

For exogenous expression of a wild-type Phax transgene, an intron-less Phax cDNA was subcloned into pBID-UASC-GV (Wang et al. 2012), and subsequently injected and integrated into the 51C1 PhiC31 landing site by BestGene Inc. The UAS:Phax-mVenus transgene at 51C1 was subsequently recombined into the Phax background. The armadillo promoter-GAL4 driver P{GAL4-arm.S}11 (Sanson et al. 1996) was also recombined into the Phax mutant background.

RNA-seq and bioinformatics analyses

RNA was isolated from 74 ± 2-h-old wild-type and snRNP biogenesis mutant larvae by homogenization in TRIzol (Invitrogen), according to the manufacturer's protocol. Staged animals were isolated from either the wandering third instar stage of development or just before bubble migration in early tan pre-pupae (<12 h after puparium formation). RNA was DNased with Amplification Grade DNase I (Invitrogen) and TURBO DNase (Ambion). A TruSeq RNA Sample Preparation Kit v2 (Illumina) was used for: poly(A)-enrichment or substituted here with Ribo-Zero (Epicentre) rRNA subtraction, barcoding for multiplexing, and cDNA library preparation. Paired end (2 × 50) sequencing was performed on an Illumina HiSeq 2000 platform. TopHat and Cufflinks were used for gene expression analyses, according to the bioinformatic pipeline from Trapnell et al. (2012). DaPars was used to identify differences in RNA length and alternative polyadenylation (Masamha et al. 2014; Xia et al. 2014). MISO (Katz et al. 2010) was used to measure changes in annotated alternative-splicing events. For comparison, we used a “strong” Bayes factor cutoff for our poly(A)-selected RNA-seq analysis (Jeffreys 1961). We used a “positive” Bayes factor as our cutoff for the total rRNA(−) RNA-seq analysis (Kass and Raftery 1995). For developmental comparison, all publicly available RNA-seq data from L2 to L3 puff stage 7–9 were downloaded from the NCBI Sequence Read Archive (SRA) and converted to fastq format with the SRA toolkit (Graveley et al. 2011) (http://www.modencode.org/). The heatmap was generated with the “pheatmap” package from bioconductor (http://www.bioconductor.org).

Northern blotting

Northern analyses were performed using standard protocols (Garcia et al. 2009). Briefly, RNA was isolated as indicated above and ∼300 ng of total RNA was separated on 10% Novex TBE Urea gels (Invitrogen), transferred to GeneScreen Plus charged nylon membranes (PerkinElmer), dried, and UV cross-linked. Northern blots were probed with 5′-end-labeled oligonucleotides for the respective snRNAs, labeled with T4 polynucleotide kinase (NEB) and [γ-32P]ATP (PerkinElmer). The snRNA-specific oligonucleotide probe sequences are listed in Supplemental Table S22. Hybridizations were performed at 65°C, and sequential washes in 2× SSC and 0.33× SSC were performed at 60°C. Damp blots were exposed to storage phosphor screens (GE Healthcare) and analyzed with a Typhoon TRIO+ (GE Healthcare). Bands were quantified using ImageQuant TL (GE Healthcare), and P-values for three biological replicates were determined with a one-tailed Student's t-test for unequal variance of two samples.

RT-PCR analysis

Semi-quantitative RT-PCR was performed with intron flanking oligos, using the 2× Apex Taq Master Mix (Genesee Scientific). Quantitative RT-PCR experiments were carried out as in Praveen et al. (2012). Briefly, random hexamers were used to prime the reverse transcription of isolated RNA according to the manufacturers protocol SupersScript III (Invitrogen). Real-time PCR reactions of cDNA were carried out on a 7900HT Fast Real-time PCR machine (Applied Biosystems), according to the Maxima SYBR Green/Rox qPCR Master Mix (2×) (Thermo Scientific) two-step protocol. Oligo pairs amplified cDNAs from intron inclusive versus spliced mature mRNAs. Three biological replicates were tested for each genotype. The ΔΔCt method was used to quantify differences, and P-values were determined with a one-tailed Student's t-test for unequal variance of two-samples. Gene-specific primer sequences are listed in Supplemental Table S22.

Western blotting

Western analyses were performed using standard protocols and modifications to the methods of Praveen et al. (2012). Briefly, larval protein lysates were prepared by crushing the animals in lysis buffer (50 mM Tris–HCl [pH 7.5], 150 mM NaCl, 1 mM EDTA, 1% NP-40, 0.1% SDS, 1% deoxycholate) with 5× protease inhibitor cocktail (Invitrogen) and clearing the lysate by centrifugation. The anti-dUTPase antibody was a generous gift from B.G. Vértessy, and it was used at a dilution of 1 to 50,000 (Békési et al. 2004). Affinity purified anti-dSMN was used at a 1 to 2500 dilution (Praveen et al. 2012). Anti-α-tubulin antibody (Sigma) was used at a 1 to 50,000 dilution.

Cell culture and RNAi

The S2 cell line was obtained from the Drosophila Genome Resource Center and cultivated according to published protocols (Rogers and Rogers 2008). Gene-specific RNAi knockdown was also performed as described in Rogers and Rogers (2008). Briefly, S2 cells were maintained in Sf-900 II SFM media (Gibco) supplemented with 1× Pen Strep (Gibco). For RNAi, 1 mL of (1 × 106 cells/mL) was aliquoted into each well of a six-well tissue culture plate. Double-stranded RNA (dsRNA) was in vitro transcribed overnight (∼16 h) from PCR templates with the T7 MEGAscript kit (Invitrogen), and dsRNA products were subsequently treated with TURBO DNase (Ambion), according to the respective manufacturers’ protocols. Oligonucleotide sequences for amplifying the PCR templates are listed in Supplemental Table S22. PCR templates were amplified from the following plasmids: Smn, pAFW-SMN; Prp6/CG6841, cDNA clone OFa18909 (GenScript); or Prp8, cDNA clone OFa09689 (GenScript). Purified dsRNA was phenol–chloroform extracted, ethanol precipitated, resuspended in RNase-free water (Ambion), and quantified by spectrophotometer. Approximately 15 μg of dsRNA was added to 1 × 106 cells per each well of a six-well plate. The dsRNA was added once each day for three successive days. On the fourth day, cells were pelleted and resuspended in TRIzol (Invitrogen). RNA and proteins were extracted according to the manufacturer's protocol.

DATA DEPOSITION

The poly(A)-RNA-seq of wild-type Oregon-R and Smn null animals were deposited in the NCBI Gene Expression Omnibus (Edgar et al. 2002) and are accessible through GEO Series accession number GSE49587 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE49587). The additional data discussed in this publication are accessible through GEO Series accession number GSE81121 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81121).

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.