Long-read subcellular fractionation and sequencing reveals the translational fate of full-length mRNA isoforms during neuronal differentiation

Characterization of a transcriptome supplemented with long-read-derived novel transcripts

To investigate the relationship between alternative pre-mRNA splicing and isoform-specific mRNA translation, we capitalized on the capability of long-read sequencing to capture complete transcript structures of polyribosome-associated mRNA, without sacrificing throughput by generating both long-read and short-read Frac-seq libraries (Sterne-Weiler et al. 2013). We used human embryonic stem cells (ESCs) and neural progenitor cells (NPCs) as a model system to characterize the translated transcriptome during early neuronal differentiation. The resulting samples were the cytosol, monosome (Mono), light polyribosome (LPR; two to four ribosomes), and heavy polyribosome (HPR; five or more ribosomes) fractions (Fig. 1A). By utilizing the R2C2 method (Volden et al. 2018), our long-read libraries, with a mean read length of 2 kb and a mean library size of 500,000, were reinforced with improved base-calling accuracy (93%) and with high-confidence transcript starts and ends. Fractionation of the long reads was employed to enhance the likelihood of detecting transcripts that may be preferentially associated with distinct subcellular fractions. All long-read libraries were pooled for de novo transcriptome assembly using Mandalorion (Volden et al. 2023). Mandalorion assembled a transcriptome consisting of transcripts with long-read support of five or more. We then used SQANTI3 to quality control and filter the initial transcriptome, resulting in a long-read-derived transcriptome containing only those transcripts with 20 or more total mean normalized short-read counts across all fractions for each cell type. To this, we applied functional annotation using all three modules of the Functional IsoTranscriptomics analysis suite (de la Fuente et al. 2020). The resulting long-read-derived transcriptome was then merged with GENCODE's GRCh38.p13 release 41 primary assembly annotation (Harrow et al. 2012) to account for transcripts that were not captured by long-read sequencing. The following analyses were done in the context of this “comprehensive transcriptome” containing both annotated and long-read-derived novel transcripts.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

Experimental overview and characterization of the comprehensive transcriptome and the cytosol. (A) Schematic of the experiment and subsequent bioinformatic analysis workflow of the resulting cytosolic extract and fractionated, ribosome-associated short and long reads from ESCs and NPCs. (B) The transcriptome classified by SQANTI3-defined structural categories of spliced transcripts, including: full splice match (FSM), incomplete splice match (ISM), novel in catalog (NIC), novel not in catalog (NNC), and intergenic or fusion transcripts (other). FSM and ISM transcripts match annotated splice sites and junctions (in GRCh38.p13 release 41), whereas NIC transcripts comprise novel combinations of annotated splice sites and junctions, and NNC transcripts contain at least one unannotated splice site. (C) The transcriptome classified by productivity based on the detection of complete or incomplete open reading frames (productive or noncoding respectively), premature stop codons (NMD), and retained introns (RIs). (D) Stratification of the transcriptome by the number of isoforms and unique coding sequences per gene. (E,F) Gene-level (E) and transcript-level (F) differential expression between NPC and ESC cytosolic fractions. (G) Top 12 enriched Metascape pathways in differentially expressed genes between NPC and ESC cytosolic fractions.

The comprehensive transcriptome had short-read coverage meeting a one count per million reads (CPM) cutoff for 37,755 transcripts with 33,561 unique CDSs, arising from 13,161 genes (Fig. 1D). Of these, 5875 and 4590 transcripts from 4095 and 3176 genes were uniquely expressed in ESCs or NPCs, respectively. Transcripts were organized into SQANTI3-defined structural categories based on their fidelity to transcript structures in the GRCh38.p13 release 41 primary assembly annotation (Harrow et al. 2012): 91.2% of transcripts matched the annotation; 8.7% were considered novel (containing either novel combinations of known splice sites and junctions or at least one novel splice site); and 0.1% were categorized as either genic or fusions (Fig. 1B). Additionally, transcripts were categorized based on their productivity. We define productive transcripts as those encoding a full-length, canonical protein. Unproductive classes include noncoding (lacking a complete open reading frame), NMD, and retained intron (RI): 66.1% of transcripts were considered productive; 0.7% were predicted to be noncoding; 18.1% were classed as NMD based on the presence of a PTC; 11.5% had a RI; and the remaining 3.6% met both NMD and RI conditions (Fig. 1C).

We used Salmon (Patro et al. 2017) to pseudoalign the fractionated short reads, with an average library size of 71.5 million reads, to the comprehensive transcriptome, producing transcript-level quantification across the gradient. Using the cytosolic fraction, which represents the raw output of the nucleus, we next tested the baseline transcriptomic differences in NPCs relative to ESCs at the gene level (Fig. 1E) and at the transcript level (Fig. 1F) to reveal upregulation of NPCs and neuronal differentiation markers and downregulation of pluripotency markers. Twenty-two and 24 transcripts from genes that are PSC or NPC markers or that are associated with neuronal differentiation had distinct ribosome association profiles in ESCs and NPCs, respectively. Only three of those transcripts had distinct ribosome association profiles in both cell types, and they all associate with the Mono regardless of expression differences. In these few examples, we did not observe a compensatory relationship between ribosome association and expression levels. Metascape (Zhou et al. 2019) pathways further encapsulated these observations (Fig. 1G). Taken together, these results present the framework for an approach to integrate fractionated long and short reads to study translational control at isoform-level resolution.

Thousands of transcripts exhibit distinct association with particular subcellular fractions

To discover if mRNA transcripts have distinct ribosome association profiles, we clustered transcript-level expression trajectories across the gradient using tappAS (de la Fuente et al. 2020), revealing subpopulations of transcripts with clear enrichment in one subcellular fraction over the others (Fig. 2A; Supplemental Tables 4, 5). A subpopulation of transcripts enriched in both the Mono and LPR fractions stood out as one of the most populous subsets, making up ∼30% of transcripts with enrichment in subcellular fractions in both cell types. Thus, subsequent analyses class Mono-associated transcripts and LPR-associated transcripts as those that are exclusively enriched in those fractions, leaving the set of Mono and LPR-associated (M + L) transcripts as a standalone subpopulation. Thousands of transcripts were considered significantly enriched (log₂FC ≥ 1.0, P-value ≤ 0.05) in subcellular fractions relative to the cytosol (Fig. 2B). Overall, 7.5% and 6.8% of transcripts were significantly associated with a subcellular fraction in ESCs and NPCs, respectively. In the context of nonmutually exclusive enrichment in subcellular fractions, subpopulations of transcripts were generally dissimilar across fractions within cell type, with the exception of the Mono and LPR fractions with Jaccard similarity of 0.41 and 0.33 in ESCs and NPCs, respectively (Supplemental Fig. 1), owing to the substantial M + L transcript subpopulations. When stratified by productivity, Mono- and LPR-associated transcript subpopulations exhibited pronounced incorporation of unproductive classes relative to the cytosol. HPR-associated transcripts displayed a relative reduction of unproductive classes relative to the cytosol in ESCs, whereas an increase is observed in NPCs (Fig. 2C). These findings support the hypothesis that levels of ribosome association may correlate with levels of translatability.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Establishing transcript ribosome association profiles. (A) Clustering of transcripts by their expression trajectories across the gradient relative to the cytosol. (B) Extraction of fraction-associated transcripts based on significant enrichment (log₂FC ≥ 1.0, P-value ≤ 0.05) in the Mono, LPR, or HPR fractions relative to their abundance in the cytosol. Fractions: (C) cytosol, (M) Mono, (L) LPR, and (H) HPR. (C) Categorization of fraction-associated transcripts by productivity. (D) Differential sedimentation of three isoforms in TMEM59. Above spliced isoform models, histograms of short-read support at exons are colored by fraction. The stacked barplot summarizes the proportion of total gene expression each isoform contributes in each subcellular fraction.

AS confers functional consequences to the stability and translation of mRNAs

Because we observed transcript subpopulations with distinct ribosome association profiles, we postulated that alternative mRNA isoforms may likewise sediment discretely. To test this hypothesis, we calculated the expression of individual isoforms relative to all isoforms from the same gene and measured the difference between their gene fractions in subcellular fractions relative to those in the cytosol (Supplemental Tables 6, 7). TMEM59 is an example of a gene with three isoforms, two of which exhibit differential sedimentation in ESCs (Fig. 2D). TMEM59 expression in ESCs is composed of a M + L-associated isoform, an HPR-associated isoform, and a cytosol-associated (not differentially sedimenting) isoform. Endogenous post-transcriptional silencing of TMEM59 by miR-351 in murine neural stem cells has been implicated to promote neuronal differentiation (Li et al. 2012), although the two differentially sedimenting isoforms share all but the last base of their 3′ UTRs. But along similar lines, it may be the case that cis-regulatory differences in their 5′ UTRs and CDSs influence the isoform-specific nature of their sedimentation.

We found 3321 (26.5%) and 2254 (19.2%) genes in ESCs and NPCs, respectively, exhibiting differential isoform sedimentation (|Δ gene fraction| ≥ 0.1 and Q-value ≤ 0.05) in a subcellular fraction relative to the cytosol (Fig. 3A). Within those genes, 4906 and 3229 transcript isoforms preferentially sedimented (Δ gene fraction ≥ 0.1) with a subcellular fraction in ESCs and NPCs, respectively. These instances substantiate AS as an architect of isoform-specific translational control. We observed decreasing concordance of gene fraction changes across the gradient between cell types, with a Pearson correlation coefficient of 0.59, 0.54, and 0.09 in Mono, LPR, and HPR, respectively. In fact, 3085 transcripts in 1506 genes exhibit divergent patterns of isoform sedimentation (Supplemental Fig. 2). Together, these findings suggest that isoform-specific sedimentation is likely cell type specific, possibly owing to differences in the composition and environment of trans-acting factors.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Differential isoform sedimentation across the gradient and functional outcomes of alternative splicing (AS). (A) Volcano plots representing differential isoform sedimentation by changes in isoform gene fraction relative to the cytosol; 3321 and 2254 genes in ESCs and NPCs, respectively, exhibit significant differential isoform sedimentation (|Δ gene fraction| ≥ 0.1, Q-value ≤ 0.05). The central plot shows changes in isoform gene fraction, with Q-value ≤ 0.05, of isoforms expressed in both ESCs and NPCs. (B) The first two bar plots show the top six enriched Metascape pathways in genes containing significant instances of differential isoform sedimentation exclusively in ESCs or NPCs. The third bar plot, labeled “divergent,” depicts enriched Metascape pathways in genes displaying contrasting patterns of isoform sedimentation between ESCs and NPCs. (C) The first stacked bar plot categorizes significant AS events (|ΔΨ| ≥ 0.1, adjusted P-value or Q-value ≤ 0.05) as AS; AS coupled with translational control (ASTC), meaning splicing events that are differentially included across the gradient; NMD; and AS coupled with both translational control and nonsense-mediated decay (ASTC + NMD). The following two bar plots show the breakdown of event types comprising each category in ESCs and NPCs. (D,E) UCSC Genome Browser snapshot of long-read and short-read coverage at SRSF7 (D), exhibiting subcellular fraction–associated inclusion of a conserved retained intron, and at ATRAID (E), exhibiting subcellular fraction–associated alternative first exon usage. Fractions: (C) cytosol, (M) Mono, (L) LPR, and (H) HPR. (F) Luciferase assay measuring the translational impact of using either the distal or the proximal ATRAID 5′ UTR in HEK293 cells.

Additionally, Figure 3A illustrates that isoforms associated with the Mono and LPR fractions have a greater magnitude of gene fraction differences than the HPR fraction relative to the cytosol. This finding suggests that the Mono, LPR, and HPR fractions gradate toward increasing similarity in transcript abundance and isoform ratios with the cytosol, which is consistent with findings from other groups (Floor and Doudna 2016). Considering these results, we posit that the average number of ribosomes per mRNA in the cytosolic fractions of our ESC and NPC samples may be similar to that of the HPR fraction. Among genes displaying differential isoform sedimentation, pathways involved in chromatin organization, organophosphate biosynthesis, and phosphorylation were enriched in ESCs specifically, whereas DNA damage, stress response, and cell cycle pathways were enriched in NPCs (Fig. 3B). Genes demonstrating divergent isoform sedimentation across cell types were enriched in similar pathways, with the addition of RNA metabolism. When comparing the Gene Ontology of cognate subpopulations of subcellular fraction–associated transcripts across cell types, we observed that cell cycle, translation, and mRNA processing–related pathways were consistently represented, whereas cell type–specific pathway enrichment was much more apparent in the cytosolic fraction (Supplemental Fig. 3).

To examine the types of AS that give rise to the diversity of the transcriptome, we categorized AS events as AS (0.1 ≤ Ψ ≤ 0.9, adjusted P-value ≤ 0.05 within condition), ASTC (|ΔΨ| ≥ 0.1, Q-value ≤ 0.05 across subcellular fractions), NMD (events that introduce a PTC), and ASTC + NMD (NMD events that adhere to the mentioned cutoffs for significant ASTC events) (Fig. 3C; Supplemental Tables 8, 9); 11.8% and 7.0% of significant AS events were classified as either ASTC or ASTC + NMD in ESCs and NPCs, respectively. We found that alternative first exon, RI, and skipped exon (SE) events feature most prominently among ASTC events, whereas SEs and RIs comprise the majority of ASTC + NMD events; 2456 and 1257 CDS-altering events (A3, A5, MS, MX, and SE event types) and 526 and 359 terminal events (AF and AL event types) were linked to translational control (either ASTC or ASTC + NMD) in ESCs and NPCs, respectively. Because of the mentioned similarity between the HPR fraction and the cytosol, the majority of ASTC/ASTC + NMD events were Mono (79.1% in ESCs, 62.1% in NPCs) and LPR associated (42.5% in ESCs, 54.8% in NPCs).

In our data set, SRSF7 contained one complete and one partial RI event associated with NMD via induction of a PTC in the highly conserved SRSF7 intron 3 locus, which has been previously described to contain a conserved poison exon (Fig. 3D; Lareau et al. 2007; Königs et al. 2020). Preferential association of PTC-containing isoforms with the Mono, as well as modestly with the LPR fraction, is consistent with our understanding of NMD's effect on translation (Maquat et al. 1981; Celik et al. 2017; Nickless et al. 2017). ATRAID, a relatively poorly understood gene implicated to play roles in all-trans retinoic acid–induced apoptosis, osteoblast differentiation, and some cancer types (Ding et al. 2015; Zhang et al. 2023), demonstrates marked patterns of alternative first exon usage across the gradient (Fig. 3E). The proximal first exon of ATRAID was preferentially spliced into the Mono-associated isoform, which may indicate its reduced translation. Indeed, a Renilla-firefly luciferase assay comparing Renilla incorporating either the proximal or the distal 5′ UTR of ATRAID in HEK293 cells exhibited significant differences in fluorescence (Fig. 3F). The two isoforms of ATRAID may potentially be functionally different, as the distal first exon contains an uORF, which may encode a signal peptide with potential importance to its localization with lysosomes (Ding et al. 2015). Collectively, these results demonstrate the widespread functional impacts of AS to the cytosolic fate of mRNAs.

Intrinsic features and cis-elements correlate with transcript polyribosome profiles

Given that AS defines the cis-regulatory landscape of mature mRNAs, we sought to identify intrinsic transcript features that may encode the underlying regulatory grammar of ASTC. We define intrinsic transcript features as measurements and functional elements that are native to the sequence of a spliced transcript. To extract the relative predictive weight of features on ASTC, we employed random forest classifiers (RFCs) to perform feature selection. Across the transcriptome, we measured the length and GC-content of the transcript, CDS, and UTRs. Additionally, our feature set included 5′ and 3′ UTR motifs, uORFs, repeat sequences, coding capacity, codon frequency, the presence of PTCs, and RIs (Supplemental Table 3). Given these features, RFCs were assigned binary classification tasks to predict the correct subcellular fraction for transcripts between every combination of subcellular fraction–associated transcript subpopulations at an 80:20 train:test split using 300 estimators. From the results, we extracted permutation feature importance, with 50 repeats, and found that the number of exons, length of the CDS and UTRs, and the GC-content of the UTRs were important features for our models to correctly classify transcript polyribosome profiles. We note, however, that although our RFC models outperformed unskilled models/chance levels, the combination of generally adequate receiver operating characteristic curves with suboptimal Precision-Recall performance highlights class imbalances and the need for more data points in each fraction and further indicates the requirement of additional features beyond those included in this study (Supplemental Fig. 4).

Our findings that CDS and 3′ UTR length positively correlate with association with heavier polyribosome fractions is consistent with previous reports (Fig. 4A; Floor and Doudna 2016). This is not to be confused with ribosome density, which other groups have shown to be inversely correlated with CDS length (Arava et al. 2003; Zhao et al. 2017). Although a longer CDS can theoretically accommodate a greater number of ribosomes, the increased potential for incorporation of nonoptimal or rare codons may trigger codon usage-dependent negative impacts to translation initiation and elongation (Lyu et al. 2021). Additionally, longer CDS and transcript lengths have been observed to be negatively correlated with translation initiation rates in the context of intra-polysomal ribosome reinitiation (Rogers et al. 2017). Also, although CDS length and 3′ UTR length, individually, positively correlate with ribosome association, we did not observe clear instances of linkage between SE and AL event inclusion. En masse, inference of ribosome association based on the CDS alone is likely too simplistic to make accurate predictions.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Analysis of features correlated with ribosome association profiles. (A) The number of exons and summaries of length and GC-content, with 90% confidence interval, of the CDS, 5′ UTR and 3′ UTR of transcripts associated with subcellular fractions. Fractions: (C) cytosol, (M) Mono, (ML) Mono + Light, (L) LPR, and (H) HPR. (B) Measurement of the change in isoform gene fraction relative to the cytosol and differences in CDS, 5′ UTR and 3′ UTR length, and GC-content of differentially sedimenting isoforms relative to the dominant isoform in the cytosol. Kernel densities for all coding isoforms are drawn with a 0.2 threshold. Subplots in the bottom left of each plot summarize the relative abundance of observations in each quadrant of their respective main plot, colored by fraction. (C) HOMER-derived de novo sequence motif and known RBP motif enrichment ratios in skipped exons enriched in subcellular fractions versus skipped exons not enriched in each given fraction. “M” and “L” refer to the Mono and LPR-associated FASE sets, whereas “HvM” and “HvL” refer to the HPR-associated FASE sets relative to the Mono and LPR fractions, respectively. (D) Enrichment of motifs using the same approach as in C, but in 30 nt windows of the CDS, 5′ UTR, and 3′ UTR of isoforms exhibiting divergent sedimentation profiles across cell types. The target sets were made from isoforms that preferentially sediment with each given subcellular fraction in one cell type, and the background sets were made from isoforms exhibiting the same in the other cell type. The standard-scaled enrichment ratio colorbar is shared by C and D.

To more clearly understand changes in feature length that may impact ribosome association, we also measured the change in CDS, 5′ UTR, and 3′ UTR length relative to the dominant cytosolic isoform among isoforms belonging to genes with differentially sedimenting isoforms (termed, gene-linked isoforms). Mono- and LPR-associated isoforms displayed a clear signal of relatively shorter CDS and longer 3′ UTR, whereas HPR-associated isoforms remained largely similar or equivalent to the dominant cytosolic isoform (Fig. 4B). Relatively longer 3′ UTRs in gene-linked isoforms are connected to strong effects on ribosome association, and isoforms with 5′ UTRs ≥ 1000 nt in length have been observed to be relatively poorly ribosome associated relative to their shorter 5′ UTR-containing counterparts within the same gene (Floor and Doudna 2016). These phenomena could be owing, in part, to the potential increased inclusion of cis-regulatory elements in UTRs including miRNA target sites, uORFs, and iron-responsive elements, which can negatively impact mRNA stability and translation. Our summary analyses of GC-content measurements yielded less clear patterns in relation to ribosome association (Fig. 4B). Nonetheless, ∼30% of isoforms preferentially sedimenting in lowly ribosome-associated fractions (Mono and LPR) exhibited decreasing 3′ UTR GC-content relative to the dominant cytosolic isoform, which may support findings that relate lower 3′ UTR GC-content to increased association with P-bodies and enhanced susceptibility to miRNA targeting (Courel et al. 2019). The other 70% of lowly ribosome-associated isoforms showed the opposite characteristic regarding 3′ UTR GC-content, which is concordant with reports that suggest an inverse relationship between 3′ UTR GC-content and mRNA stability (Litterman et al. 2019). Overall, broad measurements like length and GC-content were predictive of ribosome association to some degree but appear to lack the granularity required to definitively elucidate the mechanisms underlying instances of ASTC. Comparisons between cell types regarding feature measurements also reveals cell type–specific differences in trends that indicate further layers of complexity (Supplemental Figs. 5, 6).

To look beyond length and GC-content measurements, we identified sequence motifs that are associated with subcellular fractions. To do this, we took fraction-associated skipped exons (FASEs), exons that were determined to be significantly enriched in a subcellular fraction (ΔΨ ≥ 0.1, Q-value ≤ 0.05 across subcellular fractions) relative to the cytosol, and sliced them into 30 nt windows. Each set of windows was complemented with a background set consisting of windows made from SEs that were not significantly enriched in their given subcellular fraction. The HPR fraction was tested for enrichment against the Mono and LPR fraction (HvM and HvL, respectively) to produce HPR FASE sets owing to a dearth of HPR-enriched exons relative to the cytosol. The resulting sets of FASEs included about 415 and 240 exons on average for each subcellular fraction in ESCs and NPCs, respectively.

Using HOMER (Heinz et al. 2010) on each set of windows, we discovered 111 de novo motifs, in total, that were significantly enriched (P-value ≤ 0.05, FDR ≤ 0.2) in the target sequences over their respective background sets. We used Tomtom (Gupta et al. 2007) to identify the best matches (P-value ≤ 0.05) between the de novo motifs and known RBP motifs in the Homo sapiens data set (Ray et al. 2013). We next combined the set of motifs with the CISBP-RNA H. sapiens RBP motif set (Ray et al. 2013) and used SEA (Bailey and Grant 2021) to measure their enrichment (P-value ≤ 0.05, enrichment ratio ≥ 1.1) in each set of FASEs (Fig. 4C; Supplemental Table 10). Several motifs exhibited enrichment in at least one fraction, with motif enrichment bisecting into clusters of Mono and LPR FASE sets and of HPR-associated FASE sets regardless of cell type.

We applied the same approach to the CDS, 3′ UTR, and 5′ UTR of divergently sedimenting isoforms between ESCs and NPCs to identify cell type–specific motifs that may underlie the differences in their sedimentation (Supplemental Table 11). Target sequences were generated from isoforms preferentially sedimenting with each given fraction in one cell type versus those preferentially sedimenting with the same fraction in the other cell type. Motif enrichment in these sets of sequences clustered more distinctly by cell type than by fraction, and most motifs were exclusively enriched in one cell type and not the other. We acknowledge, however, that motif analyses are limited by the fact that RBP binding specificities are often multivalent and difficult to predict. Nonetheless, we report the presence of statistically significant sequence motifs enriched in FASEs, as well as those that are differentially enriched and utilized between divergently sedimenting isoforms. Altogether, 19 of the 43 known RBPs identified as fraction or cell type specific have been previously implicated in translational control or observed to associate with polyribosomes (Supplemental Tables 10, 11; Paronetto et al. 2006; Smart et al. 2007; Markus and Morris 2009; Jin et al. 2011; Svitkin et al. 2013; Aviner et al. 2017; Ueno et al. 2019; Rizzotto et al. 2020; Sévigny et al. 2020; Zhang et al. 2021; Anisimova et al. 2023). Additional studies are necessary to test the role of these potential factors in ASTC. Because we were specifically interested in motifs related to ribosome association, we did not perform motif analysis on introns. As a whole, these results suggest that ribosome association is impacted by the composition of intrinsic transcript features, likely with combinatorial effects.