Genome-wide analysis of long noncoding RNA stability

Abstract

Transcriptomic analyses have identified tens of thousands of intergenic, intronic, and cis-antisense long noncoding RNAs (lncRNAs) that are expressed from mammalian genomes. Despite progress in functional characterization, little is known about the post-transcriptional regulation of lncRNAs and their half-lives. Although many are easily detectable by a variety of techniques, it has been assumed that lncRNAs are generally unstable, but this has not been examined genome-wide. Utilizing a custom noncoding RNA array, we determined the half-lives of ∼800 lncRNAs and ∼12,000 mRNAs in the mouse Neuro-2a cell line. We find only a minority of lncRNAs are unstable. LncRNA half-lives vary over a wide range, comparable to, although on average less than, that of mRNAs, suggestive of complex metabolism and widespread functionality. Combining half-lives with comprehensive lncRNA annotations identified hundreds of unstable (half-life < 2 h) intergenic, cis-antisense, and intronic lncRNAs, as well as lncRNAs showing extreme stability (half-life > 16 h). Analysis of lncRNA features revealed that intergenic and cis-antisense RNAs are more stable than those derived from introns, as are spliced lncRNAs compared to unspliced (single exon) transcripts. Subcellular localization of lncRNAs indicated widespread trafficking to different cellular locations, with nuclear-localized lncRNAs more likely to be unstable. Surprisingly, one of the least stable lncRNAs is the well-characterized paraspeckle RNA Neat1, suggesting Neat1 instability contributes to the dynamic nature of this subnuclear domain. We have created an online interactive resource (http://stability.matticklab.com) that allows easy navigation of lncRNA and mRNA stability profiles and provides a comprehensive annotation of ~7200 mouse lncRNAs.

Figure 1.

RNA half-life determination following transcription inhibition. (A,B) Transcript decay curves after blocking transcription in N2A with actinomycin D and measuring transcript remaining relative to a control gene by qPCR. Error bars show standard deviation. (A) Myc decay relative to Gapdh. Gapdh is a suitable control gene for transcripts that are not highly stable. Results are from four biological replicates, which were subsequently used for microarray analysis. The fitted curve was modeled by one-phase decay using nonlinear least squares regression. Myc expression was also tested in mock treated time courses, which did not show evidence of transcript decay. (B) Transcript decay curve for Zfas1 relative to Atp5e. Results from three biological replicates. No degradation is observed, and nonlinear regression supports a horizontal line fit. (C) Zfas1 expression over 32-h time course following transcription inhibition from microarray. Four biological replicates; error bars show standard deviations. Nonlinear regression was used to test model fits and supports a linear fit with a positive slope showing apparent up-regulation of expression by 32 h. (D–F) Decay curves and half-lives determined for two random mRNAs (D,E) and one lncRNA (F) transcript from microarrays. All were modeled using one-phase exponential decay. Error bars represent standard deviations.

Figure 2.

Half-lives of lncRNA and protein-coding transcripts. (A) Box-and-whisker plot of coding and lncRNA transcript half-lives. (Whiskers) 1st–99th percentile, with individual transcripts outside this shown as dots. (Box) 25th–75th percentile. Difference calculated using a nonparametric Mann-Whitney t-test. (B) Frequency distribution showing the fraction of protein-coding and lncRNA transcripts in 2-h bins. Plotted points are at the center of the 2-h bin. Only time points with 1% or more of transcripts are plotted. (C) Percentage of unstable (half-life under 2 h) lncRNA and protein-coding transcripts. Significant difference calculated using a χ² test. (D) Percentage of highly stable (half-life over 12 h) lncRNA and protein-coding transcripts. Significant difference calculated using a χ² test.

Figure 3.

Hierarchical clustering of transcript decay rates. Unsupervised hierarchical clustering of all lncRNA transcripts above the expression cut-off. Clustering was performed using cluster3 (de Hoon et al. 2004) and visualized in Java Treeview (Saldanha 2004). All transcripts were set to an expression level of 1 at 0 h, so clustering is determined only by decay rate. Transcripts that decay quickly turn black during the early time points; transcripts that show no degradation remain bright yellow. As clustering was performed on all transcripts above the expression cut-off, transcripts whose half-lives could not be determined are also included.

Figure 4.

Distance-based clustering of transcript decay rates. Transcripts are automatically clustered; with those showing indistinguishable decay profiles over the time course present in the same cluster, while transcripts with similar profiles are found in nearby clusters. The physical distance between individual clusters and between super-clusters represents the degree of difference in the decay profile. Clusters of clusters or super-clusters are created by applying the method used to form the clusters to the clusters themselves. Title gives the cluster number and the number of probes in the cluster. (X-axis) Cumulative expression of all probes in the cluster. Decay profiles are stacked bar graphs with every stack representing a separate transcript; when there are many transcripts in a cluster, the expression level of some transcripts cannot be individually visualized and are seen as areas of black (representing many transcripts). An interactive version of this figure can be found at http://stability.matticklab.com/.

Figure 5.

Effect of lncRNA features on stability. (A) Comparison of the stability of intronic versus intergenic and cis-antisense lncRNAs. Box-and-whisker plot. (Whiskers) 1st–99th percentile, with individual transcripts outside this shown as dots. (Box) 25th–75th percentile. Difference calculated using one-way ANOVA with Kruskal-Wallis nonparametric test and Dunn's post-test to compare individual annotations. (B) Frequency distribution showing the fraction of lncRNA transcripts in 2-h bins. Plotted points are at the center of the 2-h bin. Only time points with 1% or more of transcripts are plotted. (C) Percentage of unstable (half-life under 2 h) lncRNA intergenic, cis-antisense, and intronic transcripts. Significant difference calculated using χ² test. (D) Comparison of the stability of all lncRNA genomic classifications. Box-and-whisker plot. To focus on the center of the distribution, whiskers show 10th–90th percentile, with individual transcripts outside this shown as dots, and only half-lives between 0.2 h and 20 h are shown. Significant differences found between stability of classes using one-way ANOVA with Kruskal-Wallis nonparametric test and Dunn's post-test to compare individual annotations. (*) Level of significance common to all comparisons. (E) Comparison of the stability of localized transcripts. Box-and-whisker plot and statistical testing as per A. (F) Frequency distribution showing the fraction of nuclear, cytoplasmic, and nonenriched transcripts in 2-h bins. Plotted points as per B. (*) P < 0.05, (**) P < 0.01, (***) P < 0.001.

Figure 6.

LncRNAs and decay elements. (A) Comparison of the stability of single exon versus spliced lncRNAs. Box-and-whisker plot. (Whiskers) 1st–99th percentile, with individual transcripts outside this shown as dots. (Box) 25th–75th percentile. Difference calculated using a nonparametric Mann-Whitney t-test. (B) Frequency distribution showing the fraction of single exon versus spliced lncRNAs in 2-h bins. Plotted points are at the center of the 2-h bin. Only time points with 1% or more of transcripts are plotted. (C) Comparison of the stability of lncRNAs containing a major or minor polyA signal versus those with evidence of internal priming. Box-and-whisker plot and statistical testing as per A. (D) Correlation between GC% and lncRNA half-life. Spearman correlation = 0.0852 (P = 0.0145) indicates a small positive relationship between increased half-life and lncRNAs with higher GC%. Spearman correlation utilized because data is non-Gaussian. Trend line shows a semilog fit from nonlinear regression. Axes are log10-linear.

Figure 7.

Neat1 stability. (A) Comparison of lncRNA stabilities in mouse and human from Friedel et al. (2009). Error bars are standard deviations. Significant differences determined by two-way ANOVA using Bonferroni multiple comparisons. (***) P < 0.001, (****) P < 0.0001. (B) Transcript decay curve for Neat1 (both v1 and v2 isoforms) in N2A cells after blocking transcription with actinomycin D and measuring transcript remaining relative to Gapdh by qPCR. Four biological replicates. Error bars are standard deviations. Fit modeled by one-phase decay using nonlinear least squares regression. (C) Neat1 genomic locus showing v1 and v2 isoforms plus positions of PCR amplicons, FISH, and microarray probes. (D) Neat1_v2/ long isoform half-life in N2A cells, qPCR as per B. (E) Comparison of Neat1 (both v1 and v2 isoforms) (95% CI = 10–32 min) and Neat1_v2/ long isoform (95% CI = 40 min–1 h, 57 min) in N2A cells. Error bars show 95% confidence intervals. Unpaired t-test. (F) Comparison of stability of Neat1 isoforms in 3T3 cells. qPCR from three biological replicates. Error bars show 95% confidence intervals. Unpaired t-test. (G,H) Combined RNA protein FISH on N2A cells (G) and 3T3 cells (H). (Left panel) Neat1 RNA; (second panel) localization of paraspeckle protein NONO; (third panel) DAPI nuclear stain; (final panel) overlay.