Genome-wide monitoring of alternative splicing

Using our array, we monitored splicing profile changes in Drosophila SL2 cells following RNAi knockdown of four splicing regulators: the SR proteins dASF/SF2 and B52/SRp55, and the hnRNP proteins PSI and hrp48 (Supplementary Fig. 1). Each of these four well-characterized splicing regulators is highly expressed in SL2 cells, and several of them have known pre-mRNA targets. Following treatment with double-stranded RNA (dsRNA) against each splicing factor, efficient protein reduction was confirmed by immunoblot analysis using antibodies specific for each protein (Fig. 1B). RNAi knockdown of each of these splicing factors generated no obvious morphological or growth phenotype in SL2 cells, despite the fact that in Drosophila PSI, hrp48, and B52/SRp55 are essential and dASF/SF2 is likely to be essential (Wang et al. 1996, 1998; Longman et al. 2000). From each RNAi-treated sample and from control cells treated with nonspecific dsRNA, total RNA was extracted, cDNA prepared, and labeled using a protocol developed to give good coverage over the entire length of all mRNAs (Castle et al. 2003; Johnson et al. 2003). Following standard hybridization, scanning, and data extraction, each experiment and each probe signal were filtered for consistency and RNAi target specificity (for a description of the filters, see Supplemental Material). Expression ratios (red/green ratios) of RNAi knockdown of each splicing factor versus no knockdown control were computed for each probe.

Biochemical experiments demonstrate that B52/SRp55 and PSI associate with, and presumably modulate splicing of, at least dozens and perhaps hundreds of distinct pre-mRNAs (Labourier et al. 2002; Kim et al. 2003). This is also likely the case for dASF/SF2 and hrp48. Therefore, reduction of any of these factors may impede or deregulate pre-mRNA processing severely and inconsistently, rendering the array data undecipherable or irreproducible. To address this possibility, we carried out multiple RNAi knockdown experiments for each splicing factor and compared the effect of experiments using simple hierarchical clustering (Eisen et al. 1998). Clustering experiments using data aggregated for each locus, each isoform, or each splice junction (see Materials and Methods) generated nodes specific for the dASF/SF2, B52/SRp55, and hrp48 experiments (Fig. 1C; data not shown). This indicates that the global expression patterns and splicing responses assayed on this array are largely distinct for the splicing factor that has been knocked down. The results of this high-level analysis support the notion that knockdown of each of these splicing factors results in a characteristic, interpretable, and reproducible splicing response. Interestingly, the PSI replicates were more variable in this analysis. This is likely due to the relatively small number of PSI-affected splicing events (see below). Despite this variability, PSI experimental replicates did identify several mRNAs that were previously found to be associated with PSI in an embryonic nuclear RNP fraction and whose expression was deregulated in PSI-mutant flies (Supplementary Table 1; Labourier et al. 2002).

RNAi knockdown of each splicing factor causes a general decrease in processed mRNA

Five distinct, exogenous and heterologous, in vitro-transcribed, positive control RNAs were amplified and labeled along with the SL2 total cellular RNA to provide an independent assessment of any global changes in gene expression that would otherwise be masked by the normalization procedure. Following a standard assumption, we normalized the expression data from each array so that the average log expression ratio is 0; i.e., no net change in expression (see Materials and Methods). Interestingly, after normalization, the average log expression ratio of RNAi to control of nearly all positive control RNAs in each sample, across the range of expression intensities, was positive (Fig. 2A). Since identical amounts of the positive control RNAs were introduced in both the reference and experimental samples, these observations indicate an average reduction of the detectable, expressed genes on the array following knockdown of any of these four splicing factors. It is unclear to what extent this represents a biologically direct effect or a secondary consequence of target gene transcript reduction. Anecdotally, a previous screen identifying four B52/SRp55 target pre-mRNAs showed that the predominant effect of B52/SRp55 deletion was a reduction in target gene mRNA levels (Kim et al. 2003).

Open in a separate window

Figure 2.

Global microarray analysis. (A) Plot of total expression values versus log expression ratio of each positive control probe reveals a strong positive skew. Equal amounts of each positive control RNA were added to control and experimental samples. (B) The distribution of the net expression for each splice junction probe across all experiments is shown in the histogram. Cut-offs for up- or down-regulation were set at 2 standard deviations from unchanged. (C) Size of the sets of alternative and constitutive junctions that are strongly and consistently up- and down-regulated following knockdown of each splicing factor. dASF/SF2 affects the largest number of splicing events; PSI affects the smallest number. (D) Number of splicing events strongly and consistently affected by RNAi knockdown of each of two splicing factors. In parentheses are the expected sizes of each set assuming that each splicing factor affects in-dependent sets of splicing events. For each combination except B52 and PSI, the number of affected events is larger than expected by chance. (E) Antagonistically affected splice junctions. The numbers of splice events that were increased following knockdown of one splicing factor and decreased following knockdown of another are shown. Also shown is the number expected under the random model in which each splicing factor's affects are independent of the affects of the other splicing factors.

Each splicing factor affects a distinct set of splice junctions

Following data extraction and normalization, we computed log expression ratios for each probe in experimental versus control experiments. In order to determine which splicing events are affected by knockdown of which splicing factors, we calculated a value we call the net expression for each splicing junction (see Materials and Methods). The net expression for each splicing junction is the log expression ratio for that splice junction minus the average log expression ratio for all other probes on the isoforms containing that junction. This value indicates the increase or decrease in abundance of the specific splicing junction in question above and beyond any increase or decrease in abundance of the isoforms that contain the splice junction, and is designed to highlight individual splice junction changes by removing any differences in isoform expression or overall RNA levels. Because each splice junction is considered separately, this strategy should be minimally affected by incomplete data about which isoforms exist. The net expression value was calculated for all alternative and all constitutive junctions using the same formula. The distribution of the net expression values in all experiments (Fig. 2B) was used to generate cut-offs for classifying the affect of knockdown of each splicing factor on each splice junction. In order to find pre-mRNA splicing events that were strongly and consistently affected following knockdown of each splicing factor, net expression differences for each junction were compared in the biological replicates (see Materials and Methods; Supplementary Table 1).

The splicing events that were strongly (≥2 standard deviations) and consistently (in both replicates) affected by knockdown of each splicing factor were identified (Supplementary Table 1). The knockdown of dASF/SF2 affected the largest number of splicing events (319 events), while PSI affected the smallest number of splicing events (43 events) (Fig. 2C). This result is consistent with the notions that dASF/SF2 is a general regulator of alternative splicing, affecting a large number of targets, and that PSI is a more specialized regulator of alternative splicing. B52/SRp55 and hrp48 affected intermediate numbers of splice events (107 and 90 events, respectively). Interestingly, several of the splicing events detected for B52/SRp55 and for PSI were on genes previously found to interact with these splicing regulators (Supplementary Table 2).

The lists of splicing events consistently affected by each splicing factor were examined across experiments to determine which are significantly affected by knockdown of more than one of these splicing factors (Fig. 2D). The number of such splicing events was then compared with the number that would be expected by chance under a random model; i.e., assuming the splice junctions affected by one splicing factor are chosen independently of the splice junctions affected by any other splicing factor. The combination of dASF/SF2 and B52/SRp55 produced the most striking result: 22 splicing events were similarly affected by knockdown of either SR protein compared with 1.87 expected by chance (P-value ≤ 0.00001). This result is particularly interesting considering that, in vitro, SR proteins can complement one another for activity on several RNA targets (Fu et al. 1992; Zahler et al. 1992). Since dASF/SF2 and B52/SRp55 are the closest Drosophila SR protein paralogs (data not shown), this functional overlap is perhaps not surprising. However, the unique character of each SR protein is demonstrated by the fact that the majority of the splicing events strongly and consistently affected by RNAi of either dASF/SF2 or B52/SRp55 alone are not strongly or consistently affected by the other alone (22 shared targets out of 319 dASF/SF2 and 127 B52/SRp55 affected targets, respectively) (Fig. 2C,D).

PSI and hrp48 are known to coregulate alternative splicing of the P-element transposase pre-mRNA by binding to an exonic splicing silencer (Siebel et al. 1994, 1995). This analysis suggests their partnership may extend beyond P-element splicing. Of the 43 consistent and strong PSI targets, seven were found to also be strongly and consistently under the control of hrp48, whereas only 0.257 would be expected under the random model of independent effect (P-value ≤ 0.00017). Furthermore, for the 25 splice events that were consistently and strongly decreased following PSI knockdown, 21 were reduced following hrp48 knockdown by an average of -1.93 standard deviations (Supplementary Fig. 3). Similarly, for the 18 splice events that were consistently and strongly increased following PSI knockdown, 14 were increased following hrp48 knockdown by an average of 1.30 standard deviations (Supplementary Fig. 3). Therefore, nearly all of the splicing events under the control of PSI are similarly controlled by hrp48, suggesting that hrp48 may be an obligate partner for PSI. Interestingly, PSI does not appear to be an obligate partner for hrp48 as there are many hrp48 splicing events not similarly affected by PSI (Supplementary Fig. 3).

Traditionally, SR proteins and hnRNP proteins have been viewed as antagonistic partners, regulating in opposite directions, many of the same alternative splicing units (Smith and Valcarcel 2000; Black 2003). One of the best-studied models of antagonistic regulation is the HIV pre-mRNA in which binding of hnRNP proteins to cis-acting splicing silencer elements can be counter-acted by binding of SR proteins to nearby enhancer elements to regulate utilization of adjacent splice sites (Caputi et al. 1999; Jacquenet et al. 2001; Zhu et al. 2001; Caputi and Zahler 2002; Zahler et al. 2004). We analyzed antagonistic regulations by looking at splicing events that were increased following knockdown of one splicing factor and decreased following knockdown of a different splicing factor (Fig. 2E). Surprisingly, very few splicing events were found to be consistently and strongly regulated in an antagonistic relationship by any combination of these splicing factors. At a cut-off of 2.0 standard deviations, only dASF/SF2 and hrp48, which are the Drosophila homologs of the two canonical antagonistic splicing factors ASF/SF2 and hnRNP A1 (Mayeda and Krainer 1992; Cáceres et al. 1994; for review, see Black 2003), are found in more than a single antagonistic splicing event (Fig. 2E). At a more permissive cut-off, 1.5 standard deviations, more such antagonistic affects can be seen (Supplementary Fig. 4). These data indicate that antagonism between SR proteins and hnRNPs appears to be highly specific for both interacting partners (Supplementary Fig. 3) and somewhat uncommon.

RT-PCR validates the microarray results

A validation of the microarray results was performed on several individual genes by RT-PCR designed to amplify the different affected mRNA isoforms. From the analysis described above, six genes were chosen on the basis that they were previously unknown targets of any of the four factors tested and their structures were amenable to RT-PCR analysis using a single pair of primers for each target (Fig. 3; Supplementary Fig. 2). All the selected targets confirmed the microarray results with differences in expressed isoforms ranging from less than twofold (PSI-specific target CG4912) (Fig. 3C) to 25-fold (B52/SRp55-specific target CG6084) (Fig. 3A).

Open in a separate window

Figure 3.

RT-PCR validation of selected targets. Shifts in alternatively spliced isoforms predicted from the microarray analysis were monitored by RT-PCR for three different targets using oligonucleotides flanking the affected alternative splice sites. (Bottom panel) Together with the RT-PCR gel analysis, the alternative spliced junction expression computed from the microarray data are shown. The densitometry of the gels are shown on the left expressed as a log₂ ratio of the two measured isoforms. (A) CG6084 is a predicted target of the SR protein B52/SRp55. B52/SRp55 knockdown promoted skipping of the alternative cassette exon. (B) CG6143 (PEP) is a predicted B52/SRp55 target whose cassette exon is included upon knockdown of B52/SRp55. (C) CG4912 is a predicted target of the hnRNP protein PSI. The knockdown of PSI promotes skipping of the alternative cassette exon.

RNAi, RNA extraction, RNA labeling, and array hybridization

Production of dsRNA and RNAi was as described (Clemens et al. 2000). RNAi were done for 4 d by incubating 10 μg of the different dsRNA with 0.5 × 10⁶ serum-free adapted Schneider SL2 cells (Invitrogen) with addition of 10 μg of dsRNA after 48 h. At day 4, 10% of the cells were recovered and lysed in protein gel loading buffer, while the remaining cells were used for total RNA extraction using the Qiagen mRNA Easy purification kit with on column DNAse digestion following the manufacturer's protocol. In vitro transcribed RNAs of the human U17, U19, and 7Sk small RNA, as well as the human and tetrahymena telomerase RNA (300, 100, 10, 2.5, and 25 fmol each, respectively), were added as internal quality and sensitivity control to 10 μg of total RNA and were amplified, labeled, and hybridized as described (Castle et al. 2003) on a 44-k custom Agilent oligonucleotide array.

Array analysis

After hybridization, arrays were scanned and images analyzed following the manufacturer's recommendation (Agilent). Linear-LOWESS dye normalization was performed using all probes except for negative controls. The Pearson correlation coefficient was computed for all probe expression values using the two instances of each probe on the array as the x and y values and found to be >0.97 for all experiments. One B52/SRp55 experiment was removed from the analysis on the basis of a nonconsistent effect in the red channel (Supplementary Fig. 1A). Array data from each experiment were then analyzed to determine the extent of specific knockdown of each target. Experiments that failed to yield a strongly negative log expression ratio for probes of each RNAi target (log-ratio ≤ -0.2) or that also yielded a strongly negative log expression ratio for probes of a different RNAi target were discarded (Supplementary Fig. 1B). This removed the first PSI experiment, the second hrp48 experiment, and the first ASF/SF2 experiment. Remaining for further analysis were two experiments of each splice factor target. Since each probe occurs twice on the array, a single expression average was computed. Any probe whose expression at one position in either red or green channels was >150% of its expression in that channel at its other position was removed from the analysis. Hierarchical clustering of experiments using average linkage clustering was performed. A single average expression value for each locus, each isoform, or each junction was generated. Locus and isoform averages were computed by taking the average log expression ratio of all probes from each locus or from each isoform.

To separate gene-level or isoform-level expression changes from splicing changes, the average log expression ratio of all junctions for each isoform was subtracted from the expression ratio of each junction. This value, which we call the net expression of each junction, is given by the formula log expression ratio (x) - average of log expression ratios (y), where x is the splice probe in question and y is the set of all probes that are on all isoforms that contain junction x. This approach deemphasizes changes in splicing events that are correlated with other splicing events of the same transcript. However, it makes minimal assumptions about the set of isoforms present for a given locus. Therefore, it is unharmed by missing transcript isoform models since each splicing event is considered individually. The distribution of the net expression was found for each experiment and for all experiments and was found to be similar (data not shown). Therefore, we used the data compiled for all experiments to generate statistical cut-offs.

For each splicing junction probe, the net expression was compared between biological replicates of each splice factor knockdown experiment. Bona fide targets of each splicing factor are expected to be consistently affected in biological replicates whereas noise is expected to vary among replicates. Therefore, splicing junction net expression values were used to filter for those within 2 standard deviations of the net-expression value of each other. We compiled lists of splicing events that were strongly (net expression value deviated >2 standard deviations from 0) and consistently (in both biological replicates) affected by knockdown of each splicing factor (Supplementary Table 1). Comparisons of each pair of lists were used to determine overlap in splicing factor targets. The statistical significance was determined using the hypergeometric distribution with Bonferonni correction for three observations on each splicing factor. These results are given in Supplementary Table 4.

RT-PCR

For each experiment, 5 μg of RNA was extracted and primed with a dT₁₆ oligonucleotide and reverse-transcribed using SuperScript II (Invitrogen) following the manufacturer's protocol. Amplification was performed following standard conditions for 21 cycles using one-twentieth of the cDNA reaction in the presence of 5 μCi of α³²P-dCTP (3000 Ci/mmol). Oligonucleotide sequences can be obtained upon request. RT-PCR products were fractionated on a 6% acrylamide-bisacrylamide gel run in 1× TBE buffer. All gels were dried, exposed, and scanned on a Typhoon phosphorimager (Amersham-Pharmacia). Densitometry of unsaturated exposure was performed using ImageQuant (Amersham-Pharmacia).

We thank G. Crooks, M. Levine, and M. Eisen for helpful comments, as well as B. Hogg, D. Fu, C. O'Connor, and D. Cunningham for providing constructs for the in vitro transcribed control RNAs. This work was supported by a U.S. National Institutes of Health grant (R01GM61987). M.B. was the recipient of a Human Frontier Science Program long-term fellowship.