Widespread Accumulation of Ribosome-Associated Isolated 3′ UTRs in Neuronal Cell Populations of the Aging Brain

Isolated 3′ UTRs Accumulate in Specific Aged Neuronal Cell Types of the Mouse Striatum

The striatum is impacted in both Huntington’s disease and Parkinson’s disease, both age-dependent neurodegenerative conditions. To identify cell-type-specific and age-associated changes in translating mRNAs that might contribute to age-associated degeneration of striatal neurons, we profiled ribosome-associated mRNAs from the two major neuronal cell subtypes of this region: spiny projection neurons (SPNs) of the direct and indirect pathways—referred to hereafter as D1 and D2 SPNs, respectively. Cell-type-specific translating mRNAs were isolated in young and aged mice by using translating ribosome affinity purification (TRAP), which makes use of transgenic mouse lines expressing an EGFP-tagged L10a ribosomal protein gene driven by a cell-type-specific promoter (Heiman et al., 2008). D1 and D2 SPNs express almost exclusively Drd1a or Drd2 (dopamine receptors 1 and 2), respectively, the regulatory elements of which were used to drive expression of the tagged L10a protein. This system thus allows for cell-type-specific isolation of ribosome-associated mRNA from each of these cell types, avoiding contamination from non-neuronal cell types (Figure 1A).

TRAP RNA was isolated in biological quadruplicate from 42-day-old (PN42) and 2-year-old D1- and D2-tagged mice, and ribosome-associated (non-poly(A)-selected, rRNA depleted) RNA was sequenced. This sequencing procedure captures all ribosome-associated RNAs, including those without a poly(A) tail in a highly specific fashion (Figure S1A). In addition to gene expression changes across thousands of genes (see below), a meta-gene analysis of sequencing coverage profiles aligned to the stop codon indicated a dramatic increase in sequence reads at and beyond the stop codon in aged D1 SPNs (Figures 1B and S1B). This pattern was not observed in aged D2 SPNs or in samples from young mice (Figure 1C). An inspection of read coverage profiles identified a subset of genes in aged D1 SPNs that had read coverage exclusively or mostly restricted to the 3′ UTR. This pattern was highly consistent across biological replicates (Figure S1C). We confirmed a several-fold increased abundance of 3′ UTRs relative to coding regions in aged D1 SPNs for several of these genes by qPCR analysis (Figure S1D), suggesting that full-length (or ORF-only) mRNAs are reduced but not absent for this set of genes. The signatures associated with isolated 3′ UTRs observed here in aged D1 SPNs are of much greater magnitude and extend across far more genes than has been previously reported (Kocabas et al., 2015; Malka et al., 2017).

To quantify the extent of 3′ UTR enrichment of a gene in a dataset, we developed a simple metric that we call the “termination codon ratio,” R_tc, defined as the log₂ of the ratio of read density in the 3′ UTR of a gene to read density in the body of the gene (including coding region and 5′ UTR, STAR Methods), restricted to constitutive portions of the gene (Figure 1D). R_tc is analogous to measures used previously to quantify ribosome occupancy in 3′ UTRs from ribosome footprint data (Mills et al., 2016), but when applied to TRAP or RNA-seq data, R_tc measures the positional enrichment of transcript fragments rather than ribosome abundance. The distribution of R_tc values was approximately symmetrical, with a mean slightly above 0 across young D1 SPNs (mean R_tc = 0.44) and aged D2 SPNs (mean R_tc = 0.38), but was substantially skewed toward larger values in aged D1 SPNs (mean R_tc = 1.34, p < 2.2e–16, Wilcoxon rank-sum test; Figure 1E). Among aged D1 SPNs, 404 genes exhibited an R_tc ≥5, indicating at least 2⁵ = 32-fold enrichment of reads in 3′ UTR relative to gene body, versus 66 and 52 genes observed at this cutoff in young D1 SPNs and aged D2 SPNs, respectively (Figure 1F). Genes with increased R_tc represented a variety of functions and pathways (none strongly enriched) and were generally shorter, with a higher C+G content and an increased potential for 3′ UTR RNA secondary than expressed genes overall (Figures S1E and S1F).

Although the presence of 3′ UTR RNAs lacking the coding region has been observed previously in some instances (Carninci et al., 2006; Kocabas et al., 2015; Malka et al., 2017; Mercer et al., 2011), the mechanism of their biogenesis remains elusive. Analyzing total RNA sequencing data from several mouse tissues and the mouse NIH 3T3 cell line, we identified several genes, including Frat2, Hnrnpa0, and Sox12, that had read densities at least several-fold higher in the 3′ UTR than in the coding region in most or all of these samples (Figures S1G and S1H). Thus, the phenomenon of 3′ UTR-enriched genes observed so widely in aged D1 SPNs extends to some genes in other mouse cells and tissues, as observed previously (Kocabas et al., 2015; Malka et al., 2017).

Aged D1 SPNs Exhibit Increased Oxidative Damage Compared to D2 and Young SPNs

Because gene expression signatures can provide information about cellular environmental conditions, we examined gene expression profiles of young and aged mouse SPNs for clues to the conditions that give rise to 3′ UTR enrichment. Young D1 and D2 SPNs exhibited cell-type-specific expression signatures typical of those that have been observed previously by microarray and single-cell sequencing studies (Gokce et al., 2016; Heiman et al., 2008). However, 2-year-old mice exhibited three times as many differentially expressed genes as PN42 mice (2526 genes compared to 835; Figure 2A). Changes in expression levels between young and aged D1 SPNs, rather than changes in D2 SPNs, explained most of these differences, and a pathway analysis of genes differentially expressed between aged SPN subtypes showed an enrichment in nuclear-encoded oxidative phosphorylation (OXPHOS) components (Figures 2B and S2A–S2C). We also observed increased expression of the primary mitochondrial superoxide dismutase Sod2, which plays crucial roles in eliminating reactive oxygen species (ROS) and sensing oxidative phosphorylation-associated distress (Zou et al., 2017) and changes in other reactive oxygen response genes in aging D1 SPNs (Figure S2D). These observations raised the possibility that D1 SPNs might be exposed to higher levels of oxygen free radicals.

To ask whether D1 SPNs are exposed to higher levels of oxidative stress, we quantified lipofuscin aggregates, a byproduct of oxidatively damaged proteins, lipids, and mitochondria, in sectioned striata of both young and aged mouse lines expressing GFP in either a D1 -specific (CP73) or D2-specific (CP101) manner (Figure 2C). Some 95% of the neuronal cells of the striatum are SPNs, allowing us to directly compare the accumulation of lipofuscin in GFP-positive versus GFP-negative neurons in each of these mouse lines. Both aged D1 and aged D2 SPNs exhibited a significantly higher average per cell lipofuscin fraction than young SPNs (p < 1e–10, chi-square test), consistent with previous observations of lipofuscin aggregate accumulation in aging neurons (Nandy, 1971). However, at 9 weeks, 47.5% of D1 SPNs were positive for lipofuscin speckles compared to 37.8% of D2 SPNs (Figure 2D; p < 1e–10, chi-square test). Furthermore, the mean per cell lipofuscin fraction was also significantly higher in aged D1 SPNs compared to D2 SPNs (Figure 2E; p = 0.002, Wilcoxon rank-sum test). These observations support the idea that wild-type D1 SPNs are subject to greater exposure to oxygen free radical species during aging than D2 SPNs.

A Model of Isolated 3′ UTR Formation from Impaired ABCE1 Activity and Cleavage by the No-Go Decay Pathway

The association of ribosomes with isolated 3′ UTRs, implied by their isolation via TRAP, and the location of the increase in 3′ UTR enrichment close to the stop codon are consistent with a cytoplasmic, translation-linked mechanism for biogenesis of isolated 3′ UTRs. We considered a number of translation-associated mechanisms involving the 3′ UTR as possible sources of isolated 3′ UTRs. For example, nonsense-mediated mRNA decay (NMD) targets mRNAs with premature termination codons and/or long 3′ UTRs for degradation (Garneau et al., 2007). However, we observed no strong relationship between 3′ UTR length and 3′ UTR enrichment (Figure S1F), and NMD is associated with prematurely terminating ribosomes rather than 3′ UTR-associated ribosomes. Alternatively, translational read-through (i.e., translation extending through the stop codon) would yield 3′ UTR-associated ribosomes. However, read-through by itself does not produce mRNA fragments. Recent studies in yeast and mammalian hematopoietic cells have observed the accumulation of ribosomes in 3′ UTRs following depletion of the yeast ribosome recycling factor Rli1 (Young et al., 2015) and of its human homolog, ABCE1 (Mills et al., 2016). Depletion of ABCE1 leads to the increased accumulation of ribosomes on the stop codon and subsequent reinitiation of translation in the 3′ UTRs of a subset of genes (Guydosh and Green, 2014; Young et al., 2015); the extent of this phenomenon varies widely between genes for reasons that are not understood. The failure of ABCE1 to efficiently recycle ribosomes on particular mRNAs has been associated with increased rates of both read-through and reinitiation, although mechanistic details are lacking (Guydosh and Green, 2014; Mills et al., 2016; Young et al., 2015).

A model in which the formation of 3′ UTR RNA fragments results from impaired ABCE1 activity is also consistent with the observed association between oxidative stress and 3′ UTR enrichment. ABCE1 contains two iron-sulfur clusters that are essential for its activity and are acutely sensitive to oxygen free radicals (Alhebshi et al., 2012). Treatment with ROS-inducing compounds, such as paraquat (PQ), impairs ABCE1 iron-sulfur cofactor uptake, reducing its activity (Alhebshi et al., 2012). In yeast, age-induced mitochondrial dysfunction leads to impaired assembly of iron-sulfur clusters, which are synthesized in the mitochondrial matrix, and targeted mRNA decay of numerous iron-dependent genes, including Rli1 (Puig et al., 2005; Veatch et al., 2009). The induction of oxidative stress with hydrogen peroxide in yeast also increases the proportion of 3′ UTR ribosomes, consistent with defects in ribosome recycling (Gerashchenko and Lobanov, 2012).

Transcripts that fail to undergo translation termination and have stalled ribosomes are normally degraded through the PELO/HBS1-dependent No-Go decay (NGD) pathway, which induces endonucleolytic cleavage of the message just upstream of the stalled ribosome (Doma and Parker, 2006; Tsuboi et al., 2012). Cleavage of the mRNA yields two fragments, with the 5′ fragment degraded 3′-to-5′ by the exosome and the 3′ fragment normally degraded 5’-to-3′ by XRN1 (Tsuboi et al., 2012). This process has also been implicated in mRNA degradation during NMD (Arribereand Fire, 2018). A similar process of ribosome-associated endonucleolytic cleavage just upstream of ribosomes termed “ribothrypsis” has recently been reported and could be considered a generalized form of NGD (Ibrahim et al., 2018).

We hypothesized that the consistent failure of ABCE1 to recycle a ribosome from the stop codon of an mRNA and from its 3′ NGD cleavage product could produce a relatively stable 3′ UTR RNA species. Such an RNA would be protected from exonucleases at its 5′ end by a persistently stalled ribosome (and perhaps secondarily stalled ribosomes upstream) and at its 3′ end by the poly(A) tail (Figure 3A). Inhibition of ABCE1 is known to increase the rates of read-through and reinitiation (Guydosh and Green, 2014; Mills et al., 2016; Young et al., 2015), potentially giving rise to secondary RNA species. Read-through or reinitiation involves ribosome movement into the 3′ UTR, potentially yielding partial 3′ UTR fragments protected at their 5’ ends by secondarily stalled ribosomes (Figure 3A).

Depletion of ABCE1 Induces Cleavage-Associated Fragments and 3′ UTR Enrichment

To directly test our hypothesis that reduced ABCE1 activity can trigger 3′ UTR enrichment, we used small interfering RNAs (siRNAs) to knockdown Abce1 in cultured NIH 3T3 cells (Figures 3B and S3A) and performed qRT-PCR with primers to the Nanos1 gene, which showed robustly increased 3′ UTR enrichment in aged D1 SPNs. This knockdown resulted in a ~50% increase in 3′ UTR enrichment of Nanos1 after 3 days, assessed using primer pairs in the 3′ UTR versus pairs in the ORF. To test whether oxidative stress could inhibit ABCE1 and subsequently induce isolated 3′ UTRs, we treated 3T3 cells with the drug PQ, a potent inducer of ROS. Over a 5-day treatment time course with 100 μM PQ, ABCE1 protein levels decreased by ~4-fold (Figure S3B), possibly by limiting or chelating iron and inhibiting the assembly or stability of ABCE1 (Brumaghim et al., 2003), demonstrating impacts of oxidative stress on the ABCE1 protein. Furthermore, PQ-induced ABCE1 depletion was accompanied by the increased 3′ UTR enrichment of Nanos1 (measured after 72 hr; Figure 3B). Thus, both direct depletion and indirect oxidative stress-induced inhibition of ABCE1 activity can induce 3′ UTR accumulation.

To test the above model, we performed ribosome profiling in Abce1 knockdown cells (“siA”) and under induced oxidative stress by using PQ as described above. Results were compared to control siRNA (“siC”). Both PQ and siA treatments reduced ABCE1 levels by >90% (Figures S3A and S3B). Ribosome-protected RNAs were isolated from these samples and sequenced using standard techniques targeting ribosome-protected RNA fragments 15–32 nt long (STAR Methods).

To identify bulk changes in ribosome positioning in response to our treatment conditions, we generated metagene plots of the density of ribosome footprint 5′ ends from fragments 27–31 nt long, corresponding to the width of a ribosome. These plots showed the expected 3-nt periodicity over transcripts, and metagene plots centered on the stop codon showed a higher density of stop codon stalled ribosomes after knockdown of Abce1 (Figure S3C), as reported previously (Mills et al., 2016; Young et al., 2015). We also observed a substantial increase in the density of 3′ UTR ribosomes after Abce1 knockdown and PQ treatment (Figure 3C).

Our model of isolated 3′ UTR formation predicts that siABCE1 (small interfering RNA to ABCE1) or PQ treatment should induce 3′ UTR ribosomes specifically in genes susceptible to 3′ UTR enrichment compared to other genes. To test this hypothesis, we compared the density of 3′ UTR ribosomes in genes with high R_tc values in aged D1 SPNs to genes with low R_tc values in aged D1 SPNs (top and bottom 15%). The average density of 3′ UTR ribosomes induced by PQ or siABCE1 was ~3- to 4-fold higher in high Rtc genes compared to low Rtc genes and was shifted significantly higher relative to control treatments (Figures 3D and 3E). These observations support our model, indicating that 3′ UTR ribosomes are induced under conditions that trigger 3′ UTR accumulation specifically in those genes, which exhibit aging-associated increases in R_tc values.

The bulk of ribosome-protected fragments recovered from our ribosome profiling assay were 27–31 nt long, corresponding to the region protected by an 80S ribosome. However, a small proportion (~0.5%–2.5%) of recovered fragments corresponded to a second peak at ~15 nt (Figure S3D). These short fragments have been shown to result from cleavage upstream of a stalled ribosome and subsequent protection by an upstream translating ribosome, which stalls when the A site of the ribosome reaches the cleavage site (Figure 3F; Guydosh and Green, 2014; Young etal., 2015). This stalled ribosome can in turn induce mRNA cleavage upstream, resulting in a chain reaction for hundreds of nucleotides upstream of the initial cleavage event and the enrichment of short ribosome-protected fragments throughout this region.

We quantified short cleavage fragments as a proportion of the total number of recovered ribosome-protected fragments (Figure 3F). ABCE1 knockdown and PQ treatments increased the proportion of cleavage fragments by ~2- to 3-fold. A metagene plot of these short cleavage fragments centered at the stop codon exhibited a peak just upstream of the stop codon and extending ~200-nt upstream (Figure 3G). This pattern was also observed in a reanalysis of the shABCE1 (short hairpin RNA to ABCE1) knockdown in K562 cells (Figures S3E and S3F; Mills et al., 2016). Short cleavage fragments were enriched to the largest extent just upstream of the stop codon, whereas PQ and ABCE1 knockdown increased the density of short cleavage fragments both upstream and downstream of the stop codon (Figure 3G, inset), consistent with the variability observed in the precise boundaries of 3′ UTR enrichment. Similar to full-length ribosome-protected fragments, short cleavage fragments were more enriched in response to PQ or ABCE1 knockdown in the 3′ UTRs of high R_tc genes (Figure S3G). These results suggest that endonucleolytic transcript cleavage occurs frequently near stop codons and is increased under conditions of oxidative stress or direct depletion of ABCE1 for a subset of mRNAs.

Endonucleolytically cleaved mRNAs are susceptible to degradation by 5’ terminal exonuclease (TEX) treatment because they are not protected at their 5’ ends by a cap. We isolated total RNA from the above cellular conditions and sequenced this material before and after TEX treatment. In untreated samples, we identified a subset of candidate induced 3′ UTR genes with increased R_tc relative to the control in both the Abce1 knockdown and PQ conditions (Figure S3H). In the TEX-treated cells, however, no difference in R_tc was observed among control, PQ, and Abce1 knockdown treatments, consistent with high R_tc values resulting from endonucleolytic cleavage.

Together, these observations provide support for our model in which age-associated oxidative stress and the subsequent impairment of the activity of the ribosome recycling factor ABCE1 triggers 3′ UTR ribosomes, mRNA cleavage in the vicinity of the stop codon, and, ultimately, the accumulation of isolated 3′ UTRs for many genes in the aging brain. However, the observed patterns are complex and other processes may contribute as well.

Isolated 3′ UTRs Accumulate in the Aging Human Brain and Vary Between Brain Regions

To determine whether isolated 3′ UTRs accumulate in the aging human brain, we analyzed RNA-seq data from various human cells and tissues from individuals of varying ages. An analysis of the total RNA sequencing data from immunopanned human astrocytes derived from 16 individuals (Zhang et al., 2016) identified dozens of genes with read coverages several-fold higher than their 3′ UTRs, with a pronounced increase in read density beginning at or just beyond the stop codon, as illustrated for SOX9 (Figure 4A). Furthermore, the extent of 3′ UTR enrichment in SOX9 increased with age (Figure 4A, inset; p = 0.00016, F-test). More broadly, an analysis of 1,380 poly(A)-selected RNA-seq samples representing 30 different tissues collected as part of the Genome-Tissue Expression (GTEx) project (GTEx Consortium, 2015) identified brain region- and age-specific patterns of 3′ UTR enrichment. The R_tc values among brain-derived tissues were positively associated with age, with a significantly increased slope compared to non-brain tissues (Figures 4B, S4A, and S4B; p < 2.2e–16, F-test). Under our model of isolated 3′ UTR biogenesis, reduced ABCE1 activity should be associated with higher levels of 3′ UTR enrichment. Consistent with this prediction, among the 320 GTEx human neuronal tissues, ABCE1 expression was significantly lower in samples with higher mean R_tc (Figure 4C, p = 1.7e–25, F-test). ABCE1 gene expression levels also declined significantly as a function of sample age (Figure S4C; p = 2.9e–6, F-test) but correlated more strongly with sample R_tc than with age.

Among the 12 different regions of the brain analyzed by the GTEx consortium, the slope of R_tc versus age was generally positive but exhibited significant variation (p < 2.2e–16, ANOVA), with only the cerebellum exhibiting a negative correlation with age (Figure 4D). These findings suggest that, although 3′ UTR enrichment is associated with increased age, this phenomenon is much more prevalent in neuronal than non-neuronal tissues, with different neuronal cell types exhibiting this phenomenon to different extents. In the aging brain, signatures of oxidative damage accumulate in more metabolically active brain regions (Corral-Debrinski et al., 1992). In particular, regional levels of aerobic glycolysis correlate spatially with amyloid-β deposition and are associated with increased levels of oxidative stress and susceptibility to neurodegenerative disease (Vlassenko et al., 2010). We noted that the regions of the brain with the highest and lowest basal rates of aerobic glycolysis (Goyal et al.,2017)—the frontal cortex and the cerebellum, respectively—exhibited the highest and lowest levels of age-associated 3′ UTR accumulation, respectively. These observations suggest that metabolic differences may contribute to brain region-specific differences in 3′ UTR accumulation.

Neuronal Isolated 3′ UTR Accumulation Correlates with Mitochondrial Gene Expression

The aging human brain is known to experience sustained exposure to ROS (Mattson and Magnus, 2006). We sought to determine whether signatures of increased oxidative stress were associated with increased 3′ UTR enrichment. To identify genes associated with isolated 3′ UTRs, we calculated the mean R_tc value across genes in each GTEx brain sample as a summary statistic. This vector of mean R_tc values was then correlated to the vector of gene expression values in matched samples, for every gene. The resulting correlations and p values were subsequently ranked to identify genes of interest (Figure 4E). Strikingly, 45 of the 100 genes most positively correlated with mean per sample Rtc across GTEx samples were core mitochondrial genes or oxidative stress response genes, including 16 mitochondrial-encoded transcripts. Applying the same approach to RNA-seq data from cultured mouse cortical neurons treated with a panel of drugs (Figure S4D), we again observed that mitochondrial genes were among the most strongly correlated with increased R_tc (Figure S4E). In both of these analyses, although mitochondrial genes encoded on both the nuclear and mitochondrial genomes were robustly enriched, mitochondrial-encoded transcripts ranked most highly. These results confirm an association between 3′ UTR enrichment and upregulation of mitochondrial genes and oxidative stress response genes in both humans and mice.

Evidence that 3′ UTR-Encoded Peptides Originate from Isolated 3′ UTR Genes

In both human and mouse, some 3′ UTR sequences accumulate more than others (Figure 1F and 4B). By examining the 3′ UTR sequences of high R_tc genes, we observed that genes with higher R_tc were most likely to have TGA stop codons and least likely to have TAA (Figures S5A and S5B; p < 1e–10, D1 SPNs; p = 0.0043, GTEx, Cochrane-Armitage Test). A stop codon sequence can influence the fidelity of termination (Dabrowski et al., 2015) and may contribute to gene-specific differences in isolated 3′ UTR accumulation. We also observed that the median distance from the canonical stop codon to the first 3′ UTR stop codon (in any frame) was ~2-fold longer for genes with higher 3′ UTR enrichment in both humans and mice (but still much shorter than the main ORF) (Figure 5A). This observation raised the possibility that 3′ UTR ribosomes may be actively translating in some cases. In yeast, 3′ UTR-encoded peptides are produced when the ABCE1 homolog Rli1 is inhibited (Young et al., 2015).

To explore the possibility that peptides are produced from ribosomes associated with isolated 3′ UTRs in mammals, we analyzed tandem mass-spectrometry (MS/MS) data derived from post-mortem human brain surveyed across 7 different tissues (Carlyle et al., 2017). Spectra were mapped to both the UniProt protein database as well as a custom database of translated 3′ UTR ORFs (STAR Methods). As a positive control for our computational pipeline, we detected an average of 9,648 high-confidence annotated proteins per sample, consistent with previous analyses. We also detected 470 high-confidence 3′ UTR ORF-derived peptides (at a false discovery rate of 0.01), which were on average 52 amino acids long and lacked recognizable protein domains (e.g., see Figures 5B, 5C, and S5C). Genes with associated 3′ UTR peptides had significantly higher R_tc values compared to similarly expressed genes whose annotated proteins were detected by mass spectrometry (Figure 5D, p = 4.9e–7, t test), suggesting that these peptides may be produced from isolated 3′ UTRs or that the production of isolated 3′ UTRs and 3′ UTR-encoded peptides involves related mechanisms. These peptides were derived from all three ORFs relative to the upstream canonical ORF in similar proportions (Figure S5D). Thus, the active translation of isolated 3′ UTRs may be a biological consequence of increased oxidative stress, but whether the resulting peptides have function is not clear.

CONTACT FOR REAGENT AND RESOURCE SHARING

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Christopher Burge ([email protected]).

EXPERIMENTAL MODEL AND SUBJECT DETAILS

All mouse experiments were conducted with the approval and oversight of the MIT Animal Care and Use Committee.

METHOD DETAILS

QUANTIFICATION AND STATISTICAL ANALYSIS

Statistical parameters and specific tests are reported in the text and figures including the statistical significance. Significance was defined as p < 0.01. All statistical test were performed using the R programming language. No statistical methods were employed a-priori for randomization, stratification, sample size estimation, or exclusion of any samples.

DATA AND SOFWARE AVAILABILITY

Sequencing datasets generated in this study have been deposited under accession number GEO: GSE97461. All software used in this manuscript is referenced in the methods. All custom code and pipelines are available upon request.