Long-read sequencing of nascent RNA reveals coupling among RNA processing events

Abstract

Pre-mRNA splicing is accomplished by the spliceosome, a megadalton complex that assembles de novo on each intron. Because spliceosome assembly and catalysis occur cotranscriptionally, we hypothesized that introns are removed in the order of their transcription in genomes dominated by constitutive splicing. Remarkably little is known about splicing order and the regulatory potential of nascent transcript remodeling by splicing, due to the limitations of existing methods that focus on analysis of mature splicing products (mRNAs) rather than substrates and intermediates. Here, we overcome this obstacle through long-read RNA sequencing of nascent, multi-intron transcripts in the fission yeast Schizosaccharomyces pombe Most multi-intron transcripts were fully spliced, consistent with rapid cotranscriptional splicing. However, an unexpectedly high proportion of transcripts were either fully spliced or fully unspliced, suggesting that splicing of any given intron is dependent on the splicing status of other introns in the transcript. Supporting this, mild inhibition of splicing by a temperature-sensitive mutation in prp2, the homolog of vertebrate U2AF65, increased the frequency of fully unspliced transcripts. Importantly, fully unspliced transcripts displayed transcriptional read-through at the polyA site and were degraded cotranscriptionally by the nuclear exosome. Finally, we show that cellular mRNA levels were reduced in genes with a high number of unspliced nascent transcripts during caffeine treatment, showing regulatory significance of cotranscriptional splicing. Therefore, overall splicing of individual nascent transcripts, 3' end formation, and mRNA half-life depend on the splicing status of neighboring introns, suggesting crosstalk among spliceosomes and the polyA cleavage machinery during transcription elongation.

Figure 1.

Transcriptome analysis of S. pombe chromatin reveals cotranscriptional splicing activity. (A) Schematic of nascent RNA (nRNA) purification from chromatin for short- (RNA-seq) and long-read RNA sequencing (LRS). (B) Enrichment of genomic DNA (DNA) and nRNA in the chromatin fraction and depletion of rRNA (18S and 28S) and tRNA revealed by gel electrophoresis and staining with GelStar (Lonza). Western blot analysis with antibodies specific for chromatin-associated proteins Pol II and Histone 3 (H3) and cytoplasmic marker proteins GAPDH and RPL5. (C) Nascent and mRNA-seq read coverage (RPM) over a three-intron gene and a convergent intronless gene. The pooled coverage from three biological replicates for each cellular fraction is shown. To assess splicing levels in nRNA and mRNA, splicing per intron (SPI) was calculated from intron junction reads (Herzel and Neugebauer 2015). SPI values are shown for these representative introns underneath the RNA-seq coverage track. (D) Cumulative SPI distribution for S. pombe introns (nRNA, n = 4481 introns; mRNA, n = 2181). Mean values (line) with standard deviation (shading) are shown for three biological replicates. Gray dashed lines indicate the median splicing levels in the two populations (nRNA 0.59, mRNA 0.95). The inset shows how the SPI is calculated from the number of “spliced” and “unspliced” junction reads spanning a particular intron.

Figure 2.

Pre-mRNA features that correlate with the extent of cotranscriptional splicing. (A) Cotranscriptional splicing levels differ based on gene position. The box plot shows the distribution of nRNA SPIs for the group of single intron genes and first, internal (second, third, and other), or last introns in multi-intron genes. The box width corresponds to the respective group size. (B) mRNA splicing levels differ between introns in different gene positions. Analogous data representation as in A. (C) One quarter of introns are significantly less or more spliced than the next downstream (3′) intron in nRNA, as depicted in a volcano plot (three biological replicates). (D) 15/42 analyzed gene architecture features correlate significantly with differentially spliced intron pairs (sequence-based in black font and RNA-seq-derived in gray font). The smaller intron position for “low” spliced introns in a pair (first intron – 1, second intron – 2, etc.) is consistent with enrichment of first introns in the “5′ less spliced” group. The median modified Z-score is shown for each feature with significant difference between the “low” and “high” groups, and the respective negative log₁₀ of the Bonferroni-corrected P-value is given. For two features, no change in the median modified Z-score is visible. The respective feature distribution is presented as a box plot below. Asterisks indicate significance of direct neighbors according to the Wilcoxon rank-sum test: (*) P < 0.05, (**) P < 0.01, (***) P < 0.001, (****) P < 0.0001 after Bonferroni-correction, in A, B, and D.

Figure 3.

Single molecule long-read sequencing reveals predominantly “all or none” splicing of multi-intron transcripts. (A) Full-length transcripts mapping to nine genes with more than one intron are aligned underneath each gene diagram and color-coded according to their splicing profile. Thin lines in transcripts indicate that the intron sequence is absent (due to splicing). Transcripts overlapping with <2 introns (“others,” gray) cannot be used to analyze the order of intron splicing. The inset shows five transcripts of gene SPBC428.01c, highlighting two subclasses of partially spliced transcripts and that transcript 3′ ends refer to Pol II position. (B) All unspliced, all spliced, and partially spliced fraction in the transcriptome. The large pie chart shows the proportion of the three splicing categories from nRNA LRS. The small pie chart depicts the fraction of partially spliced transcripts, which show intron removal “in order” (all introns are spliced upstream of an unspliced intron in a particular transcript; compare zoom-in of SPBC428.01c in A), “not in order” (at least one intron is unspliced upstream of a spliced intron in a particular transcript), or a mixed pattern (compare inset of SPBC428.01c in A).

Figure 4.

High degree of splicing co-association in multi-intron transcripts. (A) Predicted and observed “all unspliced,” “all spliced,” and “partially spliced” fractions in the transcriptome. Left panel: Prediction of splicing categories from nRNA-seq, assuming each intron is spliced independently of neighbors. Predicted values were calculated for the first two introns of genes with two and more introns for nRNA. Middle panel: Proportion of the three splicing categories from nRNA LRS as in Figure 3B. Right panel: Proportion of the three splicing categories from total mRNA LRS (data from Kuang et al. 2017). (B) Predominant “all or none” splicing in multi-intron genes. The predicted fraction of “all or none” splicing (fraction of all spliced or unspliced transcripts calculated from nRNA-seq, assuming splicing independence) is plotted against the observed fraction of “all or none” splicing for 100 genes with 10 or more LRS transcripts. (C) Co-association of splicing is close to the maximal possible co-association value. A co-association score was calculated as the log₂-fold change of the observed (LRS) to predicted (RNA-seq) fraction of “all or none” splicing. The maximal possible co-association score for 2-, 3-, and 4-intron genes was calculated [log₂(1/gene splicing)] and plotted versus the mean gene splicing value (solid line). Co-association scores for individual genes fall below or on top of this line, suggesting maximal co-association for most genes (median 86%, calculated as the ratio of the observed co-association score over the maximal co-association score at a particular gene splicing). Outliers are indicated as lighter points and can be explained by inconsistencies in gene annotation.

Figure 5.

Inhibition of pre-mRNA splicing increases “all unspliced” nascent transcripts and reduces mRNA levels. (A) Left panel: Schematic of the workflow in the temperature-sensitive prp2-1 mutant and WT strains at nonpermissive temperatures. Right panel: RT-PCR shows increased levels of unspliced nascent RNA in the prp2-1 mutant after 2 h of growth at 37°C for two introns compared to the WT 972h- strain. −RT control was loaded in lanes (empty) adjacent to the +RT samples. (B) Nascent RNA WT and prp2-1 LRS read coverage over two genes. (C) Bar plot indicating splicing frequency for single intron genes in prp2-1 relative to WT at 37°C calculated from LRS data. The standard deviation from three biological replicates is given, and asterisks indicate significance according to the Student's t-test ([***] P < 0.001). (D) Comparison of the “all unspliced,” “all spliced,” and “partially spliced” fractions in the transcriptome in WT and prp2-1. (E) Reduced mRNA levels for intron-containing genes compared to intronless genes in prp2-1 at 37°C. Cumulative distribution of expression changes between prp2-1 mutant and WT for intron-containing and intronless genes (data from three replicates from Lipp et al. 2015; P-value from Kolmogorov–Smirnov test between intronless and intron-containing group).

Figure 6.

Cotranscriptional splicing correlates with higher mRNA levels. (A) Box plot of nRNA and mRNA gene splicing levels after grouping according to cytoplasmic mRNA levels relative to nRNA levels ([***] P < 0.001, Wilcoxon rank-sum test). (B) Experimental outline to induce changes in gene expression upon caffeine treatment in S. pombe cells. (C) Correlation of mRNA expression values between the two conditions identifies 1477 differentially expressed genes (FDR-adjusted P-value < 0.05, FDR ≤0.05; 566 of those are intron-containing). (D) Cumulative distribution of nRNA gene splicing differences between caffeine treatment and control. Only genes without significant changes in nRNA levels but significant differences in mRNA expression were considered ([****] P < 0.0001, Kolmogorov–Smirnov test between “mRNA expression up” [n = 52] and “mRNA expression down” [n = 65] group).

Figure 7.

Coupling between cotranscriptional splicing, polyA site cleavage, and mRNA stability. (A) Examples of intron-containing genes with unspliced transcripts extending over the polyA site. Black triangles mark the polyA cleavage sites. Full representation of all sequenced transcripts in Supplemental Figure S8B. (B) Pie charts reflecting the fraction of spliced and unspliced transcripts in single or multi-intron genes with 3′ ends either within the gene body or downstream from the polyA site cleavage site. (C) 3′ end profiles downstream from annotated polyA cleavage sites for different transcript classes. Data are binned in 20-nt intervals and normalized to the first bin (−20 nt – 0 nt from polyA cleavage site). (D) WT and Δrrp6 total RNA-seq read coverage over the same example gene as in A is shown (log-scale). Counts per nucleotide were normalized to library size. Data reanalyzed from Zhou et al. (2015). The inset zooms into the region downstream from the annotated polyA site. (E) Nascent transcript levels with 3′ ends extending over the polyA site are increased in the exosome mutant Δrrp6. RT-qPCR from Δrrp6 and WT strains confirmed higher levels of nascent RNA uncleaved at the polyA cleavage site, using qPCR primers to generate amplicons (black line above gene diagram) bridging (polyA) or downstream (post-polyA) from the polyA cleavage site (RT with random hexamers). SDs from four biological replicates are given. Asterisks indicate significance according to Student's t-test: (**) P < 0.01, (***) P < 0.001.