Metabolic Labeling of RNAs Uncovers Hidden Features and Dynamics of the Arabidopsis Transcriptome[CC-BY]

Metabolic RNA labeling coupled to RNA sequencing detects low-abundance RNAs and RNA processing intermediates, allowing for inhibitor-free RNA stability measurements in Arabidopsis seedlings. Transcriptome analysis by RNA sequencing (RNA-seq) has become an indispensable research tool in modern plant biology. Virtually all RNA-seq studies provide a snapshot of the steady state transcriptome, which contains valuable information about RNA populations at a given time but lacks information about the dynamics of RNA synthesis and degradation. Only a few specialized sequencing techniques, such as global run-on sequencing, have been used to provide information about RNA synthesis rates in plants. Here, we demonstrate that RNA labeling with the modified, nontoxic uridine analog 5-ethynyl uridine (5-EU) in Arabidopsis (Arabidopsis thaliana) seedlings provides insight into plant transcriptome dynamics. Pulse labeling with 5-EU revealed nascent and unstable RNAs, RNA processing intermediates generated by splicing, and chloroplast RNAs. Pulse-chase experiments with 5-EU allowed us to determine RNA stabilities without the need for chemical transcription inhibitors such as actinomycin and cordycepin. Inhibitor-free, genome-wide analysis of polyadenylated RNA stability via 5-EU pulse-chase experiments revealed RNAs with shorter half-lives than those reported after chemical inhibition of transcription. In summary, our results indicate that the Arabidopsis nascent transcriptome contains unstable RNAs and RNA processing intermediates and suggest that polyadenylated RNAs have low stability in plants. Our technique lays the foundation for easy, affordable, nascent transcriptome analysis and inhibitor-free analysis of RNA stability in plants.


INTRODUCTION
Transcriptomes are highly dynamic and vary greatly among different cell types, environmental conditions, and developmental stages. Plant transcriptome analysis by RNA sequencing (RNAseq) has become an affordable tool for understanding plant development and adaptive growth, for assessing the effects of specific mutations, and for improving genome annotation. Various RNA-seq technologies have not only been used for transcript abundance measurements in plants but also for detecting and annotating transcript sequence variation generated by alternative splicing, differential promoter usage, or the selection of different transcription termination and polyadenylation sites (Sherstnev et al., 2012;Morton et al., 2014;Filichkin et al., 2015;Kappel et al., 2015;Abdel-Ghany et al., 2016;Cartolano et al., 2016;Sun et al., 2017;Tokizawa et al., 2017;Ushijima et al., 2017;Vaneechoutte et al., 2017;Zhang et al., 2017). RNA-seq analysis of mutants impaired in RNA turnover or the enrichment of specific RNA fractions prior to library preparation has also provided important information about the nature, stability, and processing of various RNA molecules as well as their modifications (Drechsel et al., 2013;Schwarzl et al., 2015;Vandivier et al., 2015;Li et al., 2016;Shen et al., 2016;Cui et al., 2017;David et al., 2017;Sorenson et al., 2018).
Most transcriptomic studies in plants focus on detecting steady state transcriptomes, which are defined by the rate of RNA synthesis and degradation. Although steady state transcriptomes provide essential information about RNA composition, information about RNA synthesis and degradation cannot be deduced from steady state transcriptomes. Furthermore, the presence of long-lived and abundant RNAs in steady state transcriptomes limits the detection of less abundant and unstable RNAs including RNA processing intermediates. Certain techniques take such challenges into account. Among these, global run-on sequencing (GRO-seq) has been applied in plants (Erhard et al., 2015;Hetzel et al., 2016;Liu et al., 2018). For GRO-seq analyses, labeled nucleotides such as bromouridine (BrU) are added to isolated nuclei. Transcriptionally engaged RNA polymerase II then completes the transcription of its bound genes and incorporates BrU into the nascent RNAs (Lopes et al., 2017). BrU allows labeled RNAs to be affinity purified with antibodies, followed by sequencing analysis. GRO-seq facilitates the detection of transcriptionally active genes and allows the relative contributions of RNA synthesis and degradation to be unraveled compared to the steady state transcriptome. GRO-seq is, however, an in vitro technique, which limits its applications, including expression analysis under different environmental conditions. Furthermore, RNA processing and degradation cannot be studied in standard GRO-seq experiments. Another technique that specifically enriches nascent RNA molecules is native elongating transcript sequencing (Net-seq; Churchman and Weissman, 2011). For Net-seq analysis, engaged RNA polymerase complexes are immunopurified, followed by the isolation of nascent RNA molecules from the RNA-DNA-polymerase ternary complex. This method can be used to estimate polymerase dynamics and RNA processing (Mayer et al., 2015;Nojima et al., 2015;Harlen et al., 2016;Mayer and Churchman, 2016;Zhu et al., 2018). Net-seq can also be used to detect unstable transcripts that escape detection using conventional RNA-seq (Marquardt et al., 2014;Wery et al., 2018).
Another set of techniques for studying transcriptome dynamics and homeostasis uses metabolic RNA labeling. This approach provides information about neosynthesized RNAs during the labeling period, including nascent RNAs, unstable RNA species, and RNA processing intermediates (Rabani et al., 2011). For metabolic labeling of RNAs, cells are incubated with nucleotide analogs such as BrU, 4-thiouridine (4-SU), or 5-ethynyl uridine (5-EU). These analogs are incorporated during transcription instead of unlabeled uridine. Labeled RNAs can be purified with specific antibodies or by adding a biotin tag to the RNA molecule by click chemistry, followed by streptavidin affinity purification (Jao and Salic, 2008;Duffy et al., 2015;Yildirim, 2015;Sawant et al., 2016;Herzog et al., 2017). Such pulse-labeling experiments using 4-SU combined with DNA array hybridization have revealed RNA synthesis rates in Arabidopsis (Arabidopsis thaliana; Sidaway-Lee et al., 2014). The labeling of nucleolar RNAs with 5-EU has been applied to visualize plant nucleoli (Dvo ráčková and Fajkus, 2018). Although RNA labeling with base analogs occurs in living cells, its current application is mostly limited to cell cultures. In addition, some base analogs have toxic effects (Burger et al., 2013). DNAcontaining organelles pose an additional problem for metabolic labeling of RNA given their additional membrane barriers. Still, successful labeling of RNAs has been reported for mitochondria in human cell lines using 4-SU or 5-EU (Borowski and Szczesny, 2014;Nguyen et al., 2018). By contrast, RNA labeling using modified bases has not been reported for plant organelles.
Metabolic labeling of RNAs also allows RNA stability to be determined. Traditionally, RNA half-lives are determined by following RNA degradation over time after blocking transcription using specific molecules. These experiments use cordycepin, a-amanitin, or actinomycin D to inhibit RNA synthesis. The disadvantages of these methods are that the specific chemicals are often toxic and that they do not exert their full inhibitory capacity immediately after treatment, causing a lag phase before transcription is blocked efficiently (Wada and Becskei, 2017). Transcription and RNA degradation are also tightly linked (Wada and Becskei, 2017); thus, blocking transcription in turn can positively affect RNA stability and biases RNA half-life calculations (Braun and Young, 2014). Inhibitor-free strategies for RNA stability measurements, such as 59-bromouridine immunoprecipitation chase (BRIC)-seq, start with pulse labeling using a nucleotide analog, followed by chasing periods of various lengths using unlabeled nucleotides (Tani et al., 2012b;Imamachi et al., 2014). Transcriptome analysis after chasing periods of different length provides information about changes in RNA abundance over time. RNA half-lives determined by pulse-chase experiments often result in different stability measurements compared to inhibitorbased approaches (Rabani et al., 2011;Tani et al., 2012b).
Here, we present a technique involving the use of 5-EU for in vivo metabolic labeling of plant RNAs, which allows neosynthesized RNAs to be isolated from intact plants. We demonstrate that 5-EU metabolic labeling of Arabidopsis RNAs can be used for pulsechase experiments, which allowed us to determine stabilities of polyadenylated mRNAs genome-wide without the use of chemicals to inhibit transcription.

Labeling Arabidopsis RNAs with 5-EU
First, we tested whether small molecules commonly used for in vivo RNA labeling have toxic effects in Arabidopsis. We germinated the Arabidopsis wild-type Columbia-0 (Col-0) seeds on two different concentrations of the labeling substances 5-EU, BrU, or 4-SU, as well as cordycepin, an adenosine derivative commonly used to induce transcriptional arrest. Seedlings germinated on Murashige and Skoog (MS) medium containing 4-SU or cordycepin died after germination. 5-EU or BrU had no obvious negative effect on germination or seedling development (Supplemental Figure 1). We also grew seedlings on MS plates and transferred 5-d-old seedlings to MS plates containing 5-EU, BrU, 4-SU, or cordycepin. Seedlings transferred to 4-SU or higher concentrations of cordycepin exhibited developmental arrest, did not develop true leaves, and eventually died. By contrast, 5-EU and BrU had no obvious effect on Arabidopsis development. The higher toxicity of 4-SU compared to other labeling agents has also been observed in mammalian cell cultures (Tani et al., 2012b).
Next, we conducted a set of preliminary experiments to investigate whether feeding Arabidopsis seedlings nontoxic 5-EU would lead to in vivo RNA labeling. Arabidopsis seeds were grown in liquid MS medium for 5 d in constant light. We then supplemented the growth medium with different concentrations of 5-EU. We isolated 5-EU-containing RNAs from total RNA by covalently attaching biotin to 5-EU residues, followed by purification using streptavidin-coated magnetic beads. As a negative control, we isolated RNAs from seedlings that were not treated with 5-EU but were subjected to the purification strategy for 5-EU-labeled RNAs. RNAs attached to streptavidin beads were used for on-bead cDNA synthesis and subsequent RT-PCR analysis ( Figure 1A).
RNA from plants treated with 200 mM 5-EU was much more highly enriched compared to the negative control, suggesting that we achieved labeling of endogenous RNAs in Arabidopsis seedlings with 5-EU ( Figure 1B). In a separate set of replicate experiments, we also tested for enrichment of chloroplast RNAs by RT-qPCR. We found up to 70-fold enrichment of three selected mRNAs after 5-EU treatment relative to the untreated controls ( Figure 1C; Supplemental Table 1). Short-term pulse experiments with 500 mM 5-EU showed that 5-EU-containing RNAs had already accumulated 30 to 60 min after the addition of 5-EU to the growth medium ( Figure 1D). The use of oligonucleotide pairs that flank introns revealed that spliced, mature mRNAs as well as unspliced precursor mRNA (pre-mRNAs) accumulated after 5-EU labeling ( Figure 1D). These results suggest that 5-EU not only permits the labeling of RNAs in vivo but also allows RNA processing intermediates to be enriched. Furthermore, these results indicate that the plant splicing machinery processes 5-EU-labeled RNAs.

Comparison of the Nascent and Steady State Transcriptomes
Because we observed the accumulation of labeled RNAs 60 min after the addition of 5-EU to the plant growth medium (Figure 1), Figure 1. In Vivo RNA Labeling with 5-EU.
(A) Schematic representation of the workflow for RNA labeling with 5-EU in 5-d-old Arabidopsis seedlings grown in liquid MS medium. (B) RT-PCR analysis using RNAs isolated from mock and 5-EU (200 mM)-treated Arabidopsis seedlings according to the workflow shown in (A). A conventional cDNA synthesis served as a positive control. Oligonucleotides for RT-PCR were specific for the ACTIN2 mRNA.
(C) RNAs isolated as in (B) were analyzed by RT-qPCR using oligonucleotides specific for chloroplast genes. Mean cycle threshold (Ct) values are shown for the amplification of the chloroplast psbH, psbA, and ndhF mRNA together with the corresponding SDs calculated based on three biological replicates. (D) RT-PCR analysis using RNAs from Arabidopsis seedlings treated for 0, 0.5, 1, and 2 h with 500 mM 5-EU. Three genes were amplified with oligonucleotides flanking two or three exons. Asterisks indicate unspliced pre-mRNAs. we chose this time point for strand-specific RNA-seq library preparation. Such short-term pulses might allow nascent pre-mRNAs and unstable RNAs to be isolated, which are difficult to detect in the presence of other stable RNAs in steady state transcriptomes. To prepare nascent 5-EU-labeled RNA sequencing (Neu-seq) libraries, we purified 5-EU-labeled RNAs from total RNAs as described above. The same total RNA sample was also used to prepare regular RNA-seq libraries, which we directly compared to a transcriptome from the 60-min time point. Because nascent pre-RNAs might not have been polyadenylated, we selectively removed rRNAs and randomly primed all RNAs for RNA-seq library preparation (for details, see Methods). The experiments were performed in three biological replicates. Sequencing reads from the Neu-seq libraries aligned less well than regular RNA-seq libraries (Supplemental Table 2; Supplemental  Table 3), perhaps because 5-EU-labeled RNAs are prone to accumulate more errors during library preparation and sequencing.
To compare the steady state transcriptomes and Neu-seq transcriptomes, we first calculated the transcript abundances for all genes (Supplemental Data Set 1). Transcript per million (TPM) values from the steady state and Neu-seq transcriptomes exhibited high correlation ( Figure 2A). To test whether the Neu-seq library was enriched for rare RNAs, we considered all genes with an expression value < 1 TPM as rare. When we compared the Neuseq and regular RNA-seq libraries, we identified 16,406 genes with TPM $ 1 in both library types ( Figure 2B; Supplemental Data Set 1). The steady state transcriptome contained only 69 RNAs (TPM > 1), which were not detected by Neu-seq (TPM < 1; Figure 2B; Supplemental Data Set 1). Among these, we found several small nuclear RNAs, which are components of the spliceosome (e.g., At2g04375, At1g05057, At5g61455, At2g04455, At4g06375, At5g04085). Such RNAs are probably not produced at high rates but exhibit high stability and are hence more detectable in steady state transcriptomes. On the contrary, the Neu-seq library contained 4058 RNAs (TPM > 1) that were not detectable in the regular RNA-seq library (TPM < 1; Figure 2B; Supplemental Data Set 1). Such RNAs might be nascent or unstable RNAs, which are difficult to detect in steady state transcriptomes. Gene classification analysis using MAPMAN revealed that genes detected as expressed specifically in the Neu-seq library were enriched in categories such as noncoding RNAs including microRNAs (miRNAs; Figure 2C; Thimm et al., 2004). miRNAs are derived from longer primary RNAs (pri-miRNAs). In plants, the RNase III-like enzyme DICER-LIKE1 (DCL1) rapidly processes pri-miRNAs into mature miRNAs (Achkar et al., 2016). pri-miRNAs are present at low abundance and are therefore difficult to detect. In total, the expression of 30 miRNAs could be detected in the Neu-seq transcriptomes, but not in the steady state transcriptome (example provided in Figure 2D; Supplemental Data Set 1), emphasizing the ability of Neu-seq to detect lowabundance RNAs such as pri-miRNAs.
Neu-seq libraries also contained more antisense RNAs than steady state transcriptomes (Figures 2E to 2G; Supplemental Data Set 1). The functions of antisense transcription at individual loci have been described for some Arabidopsis genes. One of the most prominent examples is FLOWERING LOCUS C (FLC). For FLC, transcription in the antisense direction is important for establishing repressive chromatin marks (Whittaker and Dean, 2017). Interestingly, we observed high accumulation of antisense RNAs at the FLC locus in the nascent RNA transcriptome, but not in the steady state transcriptome ( Figure 2E). These results suggest that RNA production at the FLC locus under our experimental conditions mainly occurs in the antisense direction. Similar results were obtained for other annotated antisense genes ( Figures 2F and 2G). To test whether antisense RNAs globally accumulate in the Neu-seq libraries compared to the steady state transcriptome, we analyzed the amounts of RNA-seq reads mapping in antisense orientation to annotated genes and calculated antisense TPM values for all genes. The ratio between sense/antisense TPMs was clearly shifted to sense TPMs, indicating that the majority of Arabidopsis genes is mainly transcribed in the sense direction ( Figure 2H). However, the ratio of sense/antisense TPMs was significantly lower in the Neu-seq library compared to the steady state transcriptome, which indicates that Neu-seq detects antisense RNAs better than standard RNAseq ( Figure 2H). Previous Gro-seq experiments also detected substantial amounts of antisense transcripts in the Arabidopsis transcriptome (Hetzel et al., 2016). These results imply that antisense RNAs are components of the Arabidopsis transcriptome and that their detection requires more sophisticated RNA labeling techniques such as Neu-seq or GRO-seq.

RNA Processing Intermediates Accumulate in the Nascent Transcriptome
Apart from rare RNAs, we also analyzed whether Neu-seq is suitable for detecting RNA processing intermediates such as unspliced pre-mRNAs. To investigate whether Neu-seq libraries are enriched in pre-mRNAs, we first classified reads into exonic, intronic, exon/intron junction, and intergenic reads using Fea-tureCounts (for details, see Methods). The Neu-seq libraries contained 2.5 times more reads mapping to intron-exon boundaries (6.8% of all reads) than steady state RNA-seq libraries (2.8% of all reads; Figures 3A and 3B). These results suggest that Neu-seq specifically covers unspliced pre-RNAs. The number of reads mapping exclusively to exonic or intronic regions did not drastically differ between the two libraries ( Figures 3A and 3B). It is likely that purely intronic reads were not more abundant in Neuseq libraries than in RNA-seq libraries because Arabidopsis introns are relatively short (on average, 165 bases; The Arabidopsis Information Resource [TAIR]), while the read length of our Illumina libraries was 150 bases. Therefore, we expect that retained introns of unspliced pre-RNAs were mainly represented by exon/intron junction reads. Few reads mapped to intergenic regions in the Neu-seq and steady state RNA-seq libraries ( Figure 3B). These findings suggest that the vast majority of transcribed regions in Arabidopsis are annotated and that only a few additional transcribed regions can be detected in Arabidopsis.
Because we observed more reads mapping to intron-exon borders in the Neu-seq libraries, we applied the replicate multivariate analysis of transcript splicing (rMATS) algorithm to detect splicing isoforms that specifically accumulated in the nascent or steady state transcriptome (Supplemental Data Set 2). The majority of differentially spliced RNA isoforms contained retained introns, which were predominately detected in the Neu-seq data set ( Figure 3C). Only a small number of RNAs containing retained (A) Comparison of TPM of all genes detected by RNA-seq and Neu-seq. Arabidopsis seedlings were treated for 1 h with 5-EU. The same rRNA-depleted RNAs served as material for the preparation of RNA-seq and Neu-seq libraries.
(B) Venn diagram depicting the number of genes with TPM > 1 detected by Neu-seq and RNA-seq. CHO, carbohydrate; misc., miscellaneous; mit., mitochondrial; OPP, oxidative pentose phosphate; org., organic acidtricarboxylic acid. (C) Characterization of 4058 genes for which expression was specifically detected by Neu-seq. Gene functional characterization was performed with MAPMAN. Error bars denote SD. ****P < 0.0001, ***P < 0.001, **P < 0.01, *P < 0.05. n.s., not significant. introns were more abundant in the steady state transcriptome ( Figure 3C). Examples of introns found by the rMATS pipeline to be significantly retained in the nascent transcriptome data set are shown in Figures 3D to 3F. The rMATS pipeline allowed us to detect splicing isoforms including alternative 59 or 39 splice sites, mutually exclusive exons, or skipped exons. We detected more RNA containing alternative 59 or 39 splice sites or skipped exons in the nascent versus steady state transcriptome ( Figure 3C). These results suggest that Neu-seq can be utilized to detect splicing intermediates and variants. (G) Overlap of Neu-seq-specific and NMD-specific alternative (Alt.) splicing events. Alternative splicing (AS) events detected in NMD-deficient mutants (Drechsel et al., 2013) were compared with alternative splicing events detected by Neu-seq. Mut. Excl., mtually exclusive; SS, splice site.
We envisioned two distinct scenarios that could explain the accumulation of RNAs with retained introns in the Neu-seq data set. First, such RNAs might contain introns that are spliced less efficiently compared to other introns. In this case, the splicing variants detected in the Neu-seq libraries would be splicing intermediates that would be further processed into mature mRNAs. The second alternative is that alternative RNA isoforms in the Neu-seq data set might be nonfunctional products of RNA splicing that accumulate only transiently and are subject to rapid degradation. The nonsense-mediated mRNA decay (NMD) pathway removes some incorrectly spliced mRNAs and thereby regulates the expression and degradation of thousands of mRNAs in Arabidopsis (Drechsel et al., 2013;Karousis et al., 2016). To distinguish between the two scenarios, we compared Neu-seq data with NMD targets that were identified by transcriptome analysis of NMD mutants (Drechsel et al., 2013). We found that 7% of the genes producing RNAs with retained introns also accumulated in NMD mutants and that 93% of all splicing variants containing retained introns identified by Neu-seq did not overlap with known NMD targets ( Figure 3G). These observations suggest that a small subset of unspliced RNAs identified by Neu-seq are targeted for elimination by the NMD pathway but that the majority of RNA splicing variants detected in Neu-seq libraries might be further processed into mature mRNAs.

Pulse-Chase Labeling of Arabidopsis RNAs and Global Decay Analysis of Polyadenylated RNAs
Metabolic labeling of RNA also allows stability measurements of RNA molecules to be performed via pulse-chase experiments. We performed 5-EU immunoprecipitation chase (ERIC)-seq, a technique similar to BRIC-seq. To ensure a high degree of labeling with 5-EU, we incubated Arabidopsis seedlings for 24 h with 200 mM 5-EU and then replaced the 5-EU-containing growth medium with growth medium containing 20 mM unlabeled uridine. Samples were collected 0, 1, 2, 6, 12, and 24 h after the start of the chase ( Figure 4A). We isolated poly(A)-RNA from total RNA because we were mainly interested in the stability of mature, polyadenylated RNAs and not in RNA degradation intermediates. Thus, this analysis mainly addressed RNA deadenylation kinetics. A recent study in yeast revealed a strong correlation between mRNA deadenlyation and mRNA stability, suggesting that our data can serve as a proxy for mRNA stability (Chan et al., 2018). We spiked poly(A)-RNA with 5-EU-labeled LUCIFERASE (LUC) RNA (for details, see Methods). The LUC spike-in control allowed us to compare the abundance of RNAs between different samples and to account for differences in the efficiency of the click-chemistry reaction, RNA purification procedures, and library preparation. 5-EU-labeled RNAs were purified from the poly(A)-RNA/LUC mixture as described above and used for strand-specific RNA-seq library preparation. The experiment was performed in biological triplicates (for details, see Methods). Gene expression values were calculated and normalized to the amount of the LUC spike-in in the respective sample (for details, see Methods). We calculated the mean TPM values of three replicates for each time point. Genes that exhibited an expression value of 0 TPM in any of the time points were removed. We also filtered out genes that accumulated to more than 125% after the beginning of the chase. In total, 19,507 genes remained for which decay rates could be determined. To compare genes with different expression strengths, we normalized all expression values relative to the expression at time point 0 (Supplemental Data Set 3). The low amount of sequencing reads from samples 24 h after the chase might have had a negative influence on ERIC-seq data quality (Supplemental Table 3). Therefore, we analyzed the ERIC-seq data set without the 24-h data point in all subsequent analyses.
We defined five clusters of mRNAs by performing K-Means clustering of gene expression values. The mRNAs in the five clusters exhibited distinct degradation kinetics ( Figure 4B). To determine the RNA half-lives and decay patterns of each cluster, we applied an exponential decay model and locally weighted scatterplot smoothing (LOWESS) to calculate the half-life time (t 1/2 ) of each cluster. Genes in cluster A exhibited the fastest decay, with an average half-life less than 1 h (t 1/2 5 0.73 h; Supplemental Figure 2). Cluster B to E contained more stable RNAs, featuring average halflives of 1.27, 2.13, 3.38, and 5.29 h, respectively (Supplemental Figure 2). The exponential decay model and LOWESS in cluster A to D yielded similar half-lives, suggesting that RNAs in these clusters decayed exponentially (Supplemental Figure 2). The exponential decay model and LOWESS resulted in different average half-lives in cluster E (5.29 versus 4.2 h). This might be caused by the unexpected abundance pattern of RNAs in this cluster: RNAs in cluster E decayed after 1 h but re-accumulated 2 h after the chase, resulting in an internal peak (Supplemental Figure 2). This unexpected behavior might be explained by reincorporation of 5-EU nucleotides, which are generated after the degradation of 5-EU-labeled RNAs. Although we chased with a 100 times higher concentration of the unlabeled nucleoside uridine, we cannot rule out the possibility that these uridine nucleosides are quickly converted into unlabeled uridine nucleotides, which are then used for transcription. High levels of reincorporation of 5-EU nucleotides might occur, especially in RNAs produced from genes with a high transcriptional rate. Importantly, decay rates for genes that fall into cluster E (total number 1114) should be interpreted carefully (Supplemental Data Set 4).
We also calculated the RNA half-lives from all genes individually (Supplemental Data Set 5). We applied two different models to calculate half-lives: an exponential decay model and a compound model. The exponential decay model is often used to describe RNA decay patterns and to predict RNA half-life values. To address potential problems with internal peaks of RNA accumulation during RNA decay, especially for genes in cluster E, we also applied compound model fitting by a convolution of Voigt curve and exponential decay curve (Shippony and Read, 1993). A Voigt model has successfully been applied to detect peaks in mass spectrometry data (Wijetunge et al., 2015). We compared the performance of the exponential and the compound model by comparing how well the fitted curve represents the actual data (Hu and Bentler, 1999). The chi 2 goodness-of-fit test revealed that chi 2 values for the exponential fitting are larger than the chi 2 values of the compound model, suggesting that the compound model described the actual data more precisely (Supplemental Figure 3). Nevertheless, half-life values determined by exponential and compound model fitting showed high correlation (Pearson correlation, 0.7; Spearman correlation, 0.74; Figure 4C). The RNA half-lives determined by the compound model were shifted to shorter half-lives ( Figure 4C). These results suggest that both models describe the decay rates of RNAs in a similar manner but that the compound model predicts shorter RNA half-life values for some genes. We provide all estimated average RNA halflives determined by compound and exponential models in Supplemental Data Set 5. In general, RNAs exhibited a wide range of predicted half-lives from minutes to a few hours (Supplemental Data Set 5). Because we found many RNAs with half-lives less than 1 h, we suggest that a higher timescale resolution within the first hours of the chase should be obtained in future experiments.
Among the 20 most highly cited genes in Arabidopsis, 12 genes are present in our data set. The other genes are probably missing because they are expressed only in specific tissues or after specific treatments. We chose these 12 genes to show mRNA decay patterns besides their half-lives and confidence intervals provided in Supplemental Data Set 6. Most of the mRNAs exhibited an almost perfectly exponential decay ( Figure 4D). Similar plots can be generated with the information contained in Supplemental Data Sets 7 and 8, which list the normalized expression values and the fitted expression values, respectively, including statistical measures (SD, SE) for each time point.
Because we isolated poly(A)-RNAs to generate ERIC-seq libraries, we cannot reliably estimate the stabilities of nonpolyadenylated RNA, such as organellar RNAs. Given that chloroplast RNAs were efficiently labeled by 5-EU, we tested whether we could detect decay by performing a separate set of experiments on samples chased with nonlabeled uridine for 0, 2, and 6 h. Total RNAs were randomly primed for cDNA synthesis, which also allows nonpolyadenylated RNAs to be detected. A synthetic, EUlabeled transcript named aadA served as a spike-in after RNA extraction for all samples, and the data were normalized based on the aadA signal. We chose three mRNAs for our RNA decay analysis by RT-qPCR, which are all stabilized in a proteindependent (i.e., potentially regulated) manner (Lezhneva and Meurer, 2004;Johnson et al., 2010;Kupsch et al., 2012). All three mRNAs showed decreasing signals over time (Supplemental Figure 4). This analysis demonstrates that it is possible to perform RNA turnover measurements in chloroplasts based on metabolic labeling.

Quality Assessment of the ERIC-Seq Data Set
We determined how half-life measurements obtained by ERICseq vary among the biological replicates. To address this issue, we analyzed each replicate separately and compared the resulting RNA half-lives from three biological replicates to obtain statistical measures (Supplemental Data Set 6). For this analysis, we further filtered out all genes exhibiting an expression value of 0 in any of the analyzed samples at any time point. In total, we calculated halflives from three biological replicates for 14,355 genes. From these values, we calculated mean RNA half-lives, their SDs, and confidence intervals (Supplemental Data Set 6). To determine the variability of the half-life measurements among all genes, we calculated the coefficient of variation (CV), which represents the ratio of the SD and the mean. The CV allows us to directly compare the variability of RNA half-life values among the three replicates. Small CVs indicate low variation, while high CVs suggest high variance among the biological replicates. To visualize this, we depicted RNAs half-life measurements from the ERIC-seq data set that exhibited a CV of 0.2, 0.3, 0.4, 0.5, 0.6, and 1.0, respectively (Supplemental Figure 5A). In general, we found that 86% of triplicate RNA half-life measurements exhibited a CV < 0.5 (Supplemental Figure 5B). However, RNA half-life data should be interpreted carefully, especially for genes exhibiting a very high CV (Supplemental Data Set 6).

Comparison of mRNA Half-Lives Determined by ERIC-Seq and Transcriptional Arrest Experiments
We compared RNA half-lives determined by ERIC-seq with previously published half-lives for Arabidopsis mRNAs determined by transcriptional arrest experiments. For this comparison, we obtained mRNA half-lives calculated in experiments using actinomycin D (Narsai et al., 2007) and cordycepin (Sorenson et al., 2018). The overall Pearson correlation of the half-lives determined by ERIC-seq and transcriptional arrest experiments were 0.33 and 0.25 (Spearman correlation, 0.29 and 0.32, respectively), which indicates low correlation between the different methodologies ( Figures 5A and 5B). Also, the half-lives determined by Narsai et al. (2007) and Sorenson et al. (2018) exhibited only a low correlation with each other (Pearson correlation, 0.36; Spearman correlation, 0.45; Figure 5C). These results indicate that half-lives determined by all three approaches can be quite different.
Interestingly, the RNA half-lives determined by ERIC-seq are much shorter than those determined after treatment with cordycepin or actinomycin D (Figures 5A and 5B). However, general trends deduced from previous studies were confirmed using ERIC-seq: Genes with no or few introns produced RNAs with shorter half-lives ( Figure 5D). Previous studies also found that RNA stability is correlated with distinct biological functions (Narsai et al., 2007, Sorenson et al., 2018. To test this notion, we performed gene classification analysis using MAPMAN on the 10% least and 10% most stable RNAs, respectively (Thimm et al., 2004). Genes that produce unstable RNAs were enriched for stress-related genes, genes involved in hormonal pathways, and genes involved in signaling ( Figure 5E). On the contrary, genes that produce very stable RNAs were specifically enriched for genes responsible for mitochondrial and plastid functions as well as primary metabolism (such as the tricarboxylic acid cycle; Figure 5F).

DISCUSSION
In vivo RNA labeling is commonly used in the yeast and metazoan RNA research community, but it is not well established in the plant field. Here, we used 5-EU, a uridine homolog, for in vivo RNA labeling in Arabidopsis seedlings. Unlike 4-SU, we did not observe toxic effects of 5-EU treatment during germination or seedling development. In addition, our results indicate that 5-EU is taken up by the plant, is incorporated into RNAs generated from nuclear and chloroplast genes, and does not interfere with RNA processing steps such as splicing. 5-EU is effectively coupled to biotin using click chemistry, allowing for the efficient purification of labeled RNA molecules. 5-EU labeling of RNAs in vivo combined with RNA-seq is a powerful tool, providing important insights into the dynamics of the Arabidopsis transcriptome.  (D) Box plots depicting the half-lives (determined by ERIC-seq) of RNAs derived from genes with zero, one, two, three, four, or more than four introns. Box plots show the median and outliers were removed. Whiskers indicate the upper and lower quartiles. A Wilcoxon-Mann-Whitney test was performed. n.s., not significant. (E) and (F) Functional gene characterization by MAPMAN. Genes that produce the 10% least stable mRNA (E) and the 10% most stable mRNAs (F), respectively, were analyzed for an enrichment of functional gene categories. All significant categories (P < 0.01) are depicted. CHO, carbohydrates; misc., miscellaneous; mito., mitochondrial; TCA, tricarboxylic acid.

Neu-Seq: From Identifying Rare and Unstable RNAs to Measuring Transcriptional Activity
We analyzed an Arabidopsis transcriptome with Neu-seq and compared it with a regular steady state transcriptome detected by standard RNA-seq. Neu-seq detected more low-abundance RNAs, including pri-miRNAs and antisense transcripts, than the standard RNA-seq experiment. Our Neu-seq libraries did not contain many more sequencing reads emerging from intergenic regions, most likely due to the high-quality annotation of the Arabidopsis genome. However, Neu-seq could still be useful for detecting unknown transcribed regions that might only be visible in certain Arabidopsis mutants defective in RNA turnover and in other plants with much less well-annotated genomes. In general, Neu-seq provides a simple alternative to GRO-seq and related techniques for identifying nascent transcripts and thus can serve as a proxy for determining transcriptional activity. Future experiments could also include much shorter 5-EU-labeling periods, which might provide an even more comprehensive view of the nascent transcriptome.
Another promising application of Neu-seq is the potential to assess chloroplast transcription. Chloroplast transcription depends on nuclear and chloroplast-encoded RNA polymerases and is regulated in response to various cellular and external signals, often via changes in the chloroplast redox state (Pfannschmidt et al., 1999;Liere et al., 2011). Chloroplast transcriptional activity has thus far been assayed using run-on techniques that depend on broken chloroplasts. This disruptive treatment precludes the analysis of the impact of changing external conditions on transcription (Legen et al., 2002). The finding that chloroplast RNAs can be labeled in vivo with externally supplied 5-EU opens up options to determine transcriptional activity under diverse growth conditions. How 5-EU traverses the double-membraned chloroplast envelope is currently unknown. It is clear that pyrimidines must be imported into chloroplasts, since the final steps of pyrimidine de novo synthesis occur in the cytosol (Witz et al., 2012). Salvage enzymes such as uracil phosphoribosyl transferase as well as plastidic uridine kinases and pyrimidine degradation enzymes are found in chloroplasts and are essential for photosynthesis and plant development (Mainguet et al., 2009;Zrenner et al., 2009;Chen and Thelen, 2011;Cornelius et al., 2011;Ohler et al., 2019). The proton antiporter PLUTO is currently the only known transport protein in Arabidopsis that mediates the import of uracil into plastids for salvage reactions and degradation (Witz et al., 2012). 5-EU is a likely substrate for this transporter as well, which would explain the success of 5-EU-based labeling of chloroplast transcripts.

Neu-Seq for RNA Processing Studies
The Neu-seq transcriptome contained 2.5 times more reads mapping to intron-exon borders compared to steady state transcriptomes, indicating that unspliced pre-mRNAs are enriched in the Neu-seq transcriptome. mRNA isoforms specific for the Neu-seq transcriptome showed a small overlap with known NMD targets. This suggests that the majority of mRNA isoforms specific for the Neu-seq transcriptome will be further processed into mature mRNAs. Most of the splicing isoforms detected by Neu-seq still contained introns, which is in line with the idea that these introns will be further spliced. We also detected more alternative 59 and 39 splicing events in the Neu-seq transcriptome than in the steady state transcriptome. An attractive hypothesis for the emergence of alternative 59 and 39 isoforms is recursive splicing (Georgomanolis et al., 2016). This noncanonical splicing reaction removes long introns in a stepwise manner. The stepwise processing steps generate low-abundance splicing intermediates, which contain alternative 59 or 39 splice sites. The extent of recursive splicing in plants is not known, but the analysis of nascent 5-EU-labeled transcriptomes in conjunction with specialized bioinformatic approaches for detecting recursive splicing will likely answer this question in the future (Pulyakhina et al., 2015). In general, Neu-seq represents a major improvement that should facilitate future studies on plant RNA processing. A more detailed analysis of RNA processing kinetics will require shorter and multiple 5-EU-labeling periods (Schwarzl et al., 2015). 5-EU-labeling experiments could be extended to the analysis of RNA splicing in RNA processing mutants, which might elucidate the exact functions of specific proteins in regulating RNA processing (Tani et al., 2012a;Maekawa et al., 2015).

ERIC-Seq for Genome-Wide Measurements of RNA Stability
Discrepancies between RNA stability measurements determined by different methods are very common, perhaps due to biases introduced by different technique for determining RNA half-lives. For all methods, certain molecules must enter the cell to execute specific functions. The rate with which a small molecule (such as cordycepin or 5-EU) is taken up, and the lag phase until the small molecule reaches sufficient levels to perform its function (e.g., block transcription), differ depending on the molecule. Thus, depending on the method and cell line used, a wide range of different half-lives are observed. In mammals, the same mRNA exhibited half-lives of more than 10 h in early experiments compared to only 10 min in later studies (Yang et al., 2003;Rabani et al., 2011). We also observed that RNA half-lives determined by ERIC-seq dramatically differ from those obtained in previous studies, which used cordycepin or actinomycin to block transcription (Narsai et al., 2007;Sorenson et al., 2018). Our results suggest that plant RNA half-lives determined by metabolic labeling experiments are generally shorter than previously reported RNA half-lives.
One possible explanation for the discrepancy between the results of transcriptional arrest and metabolic labeling experiments (such as ERIC-seq) is the crosstalk between transcription and RNA degradation, which stabilizes RNAs after blocking transcription, as observed in yeast (Sun et al., 2013;Braun and Young, 2014). Another possibility is that cordycepin and 4-SU have strong toxic effects on Arabidopsis. The application of 5-EU for RNA stability studies might also introduce certain bias. For instance, we cannot exclude the possibility that 5-EU-containing RNAs are less stable than native RNAs in vivo. However, the observation that 5-EU, unlike 4-SU and cordycepin, does not have toxic effects on Arabidopsis makes this explanation less likely. Also, the analysis of polyadenylated mRNAs versus total mRNAs might influence mRNA half-life measurements in plants. A recent study in yeast revealed that polyadenylated transcripts are slightly less stable than poly(A)-less transcripts (median mRNA half-lives 4.0 min without poly(A) selection compared to 3.6 min with poly(A) selection; Chan et al., 2018). Whether this holds true for plants will be an important topic for future research.
Apart from methodology, differences in plant growth and treatment conditions among the different studies might also influence RNA stability. The degradation of some mRNAs might also occur at random and might not be accurately assessed with any of the methods described above. Measurements of such stochastic effects would require more sophisticated methods (Baudrimont et al., 2019). Nevertheless, general conclusions based on experiments with transcriptional inhibitors still hold true for the results we obtained by ERIC-seq. For instance, mRNAs with long half-lives often encode proteins involved in photosynthesis, electron transport, or other basic metabolism (Narsai et al., 2007;Sorenson et al., 2018). mRNAs with short half-lives contain fewer introns than long-lived RNAs and often encode proteins involved in signal transduction, stress defense, or hormone signaling (Narsai et al., 2007;Sorenson et al., 2018). Given these similar patterns, we would argue that RNA half-lives measured by transcriptional arrest-based experiments yield valuable information but may greatly overestimate RNA stability. It might be possible to measure the true half-life of an Arabidopsis mRNA in the future using a combination of genome-wide techniques and single-molecule analyses such as the recently established single RNA imaging and turnover assays (Horvathova et al., 2017).
In this study, we demonstrated the versatility of 5-EU for in vivo RNA labeling in Arabidopsis seedlings and its application for mRNA detection, processing, and stability analyses. This technique allows RNA production and stability to be studied under different conditions such as temperature stress or hormone treatments. Moreover, 5-EU labeling could be used to analyze RNA splicing and processing mutants and to explore the exact function of a specific protein in regulating RNA metabolism (Tani et al., 2012a;Maekawa et al., 2015). In the future, Neu-seq and ERIC-seq could be adapted for use with RNAs that require specialized library preparation protocols, such as miRNAs or circular RNAs Marzi et al., 2016). ERIC-seq combined with methods for accurate mRNA isoform detection, including 39 end sequencing, will allow the RNA stability of specific isoforms to be determined. Finally, 5-EU could be used for labeling experiments in which unlabeled and labeled RNAs are sequenced from the same sample to calculate mRNA half-life. We are looking forward to seeing the future possibilities of 5-EU in vivo labeling of RNAs to address yet-unanswered scientific questions in the field of plant RNA metabolism.

Plant Material, Growth Conditions, and 5-EU Treatment
The Arabidopsis (Arabidopsis thaliana) wild-type seeds (Col-0) were surface sterilized using 70% ethanol. For toxicity assays, we grew Arabidopsis on half-strength MS medium containing the indicated amount of 5-EU, BrU, 4-SU, and cordycepin. Alternatively, we grew Arabidopsis on halfstrength MS medium for 5 d and transferred the seedlings to half-strength MS medium containing 5-EU, BrU, 4-SU, and cordycepin. Plants were grown under long-day conditions (16-h-light/8-h-dark cycle) under fluorescent tubes (120 to 160 mmol/m 2 /s) at 20°C/18°C day/night for 14 d.
For 5-EU-labeling experiments, 25 mg of seeds was added to 50 mL of half-strength MS in a 125-mL Erlenmeyer flask and stratified for 2 d at 4°C in the dark. Plants were grown in constant light with constant shaking on a platform shaker (130 rpm). The seedlings were grown for 5 d. Seedlings from two (for Neu-seq) or three flasks (for ERIC-seq) were pooled in a single flask for 5-EU treatment. All experiments were performed in three biological replicates. For each biological replicate, separate batches of seedlings were grown in separate flasks and were treated separately. All experiments started on the same day. For Neu-seq experiments, seedlings were treated with 500 mM 5-EU (dissolved in DMSO) for 1 h and snap-frozen in liquid nitrogen. For the ERIC-seq and RT-qPCR experiments, seedlings were treated with 200 mM 5-EU for 24 h. After the 5-EU treatment, seedlings were washed three times with half-strength MS medium and placed in halfstrength MS medium containing 20 mM uridine. Aliquots of seedlings were harvested after 0, 1, 2, 6, 12, and 24 h and snap-frozen in liquid nitrogen. See Supplemental Protocol for a detailed hands-on protocol.
The protocol was adjusted for tests of chloroplast RNA labeling and turnover. First, 30 to 35 mg of Arabidopsis Col-0 seeds was surface sterilized using 70% ethanol. The seeds were stratified in 40 mL of halfstrength MS medium in a culture flask as described above. For the 5-EU-labeled (1EU) experiments, two flasks of 5-d-old seedlings were pooled for one biological replicate experiment. For the nonlabeled (-EU) control experiments, only one flask of seedlings was used per biological replicate. Aliquots of 5-EU seedlings were harvested after 0, 2, and 6 h. Control reactions without label were performed for the 0-h time point.

RNA Extraction
RNA was extracted using an RNeasy Plant Mini Kit (Qiagen) for total RNA analysis and with a Direct-zol RNA Miniprep Kit (Zymo Research) for chloroplast RNA analysis according to the manufacturers' protocols. Remaining genomic DNA in the RNA samples was removed by DNase treatment (Thermo Fisher Scientific) followed by cleanup with an RNA Clean and Concentrator Kit (Zymo Research) according to the manufacturer's protocol. To analyze chloroplast RNA, TURBO DNase (Thermo Fisher Scientific) treatment was performed prior to cDNA synthesis, and the DNase was heat inactivated (see RT-qPCR Analysis).

Isolation of 5-EU-Labeled Nascent RNA and Illumina Library Preparation
For the Neu-seq libraries, rRNA was removed from 10 mg of DNase-treated RNA using RiboZero according to the manufacturer's protocol (Illumina). Thirty nanograms of rRNA-depleted RNA samples was used for steady state transcriptome library preparation, and the remaining RNA was subjected to 5-EU-labeled RNA isolation. We used a Click-iT Nascent RNA Capture Kit (Life Sciences) to attach biotin to 5-EU-containing RNA according to the manufacturer's recommendations. In brief, biotin molecules were covalently attached to 5-EU. Streptavidin-coated magnetic beads specifically bound to biotinylated RNAs, and unlabeled RNAs were washed away. 5-EU-labeled RNAs bound to beads and rRNA-depleted total RNAs were used as starting materials for analysis using a Total RNA Illumina ScriptSeq Plant Leaf Kit according to the manufacturer's recommendations.
For the pulse-chase experiments, poly(A)-RNAs were isolated from 5 mg of DNase I-treated RNAs using the Next Poly(A) mRNA Magnetic Isolation Module (New England Biolabs) according to the manufacturer's protocol except that RNAs were eluted from the beads in RNase-free water. Each sample was spiked with 10 pg of 5-EU-labeled LUC RNA (see below) and subjected to library preparation as described above. We observed a consistent recovery of LUC spike-in among the biological replicates, except for the last time point (Supplemental Figure 6).
Chloroplast RNAs were isolated by closely following the method described above, with the following deviations: For the RT-qPCR experiments, 2.5 mg of total RNA was used for the attachment of biotin. Each sample was spiked with 0.1 ng of 5-EU-labeled aadA RNA (see below). The immobilized beads were incubated for 5 min on a rotator in Wash Buffer 1 and subsequently in Wash Buffer 2 before performing the standard washing steps of the Click-it protocol. Furthermore, low nucleic acid binding tubes were used for all incubations. For two of the three biological replicates, instead of the Dynabeads MyOne Streptavidin T1 magnetic beads of the Click-iT kit, Dynabeads MyOne Streptavidin C1 magnetic beads (Thermo Fisher Scientific) were used.

RT-qPCR Analysis
Five-EU-labeled (1EU) and non-labeled RNAs (-EU) bound to the beads were resuspended in RNase-free water and treated with TURBO DNase (Thermo Fisher Scientific). The DNase was heat inactivated, and reverse transcription was performed with random hexamer primers using a Re-vertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific) according to the manufacturer's protocol. The cDNA was diluted 1:10 for control and 1:100 for time-course experiments. RT-qPCR was performed with an oligonucleotide concentration of 300 nM using Luna Universal qPCR Master Mix (New England Biolabs;de Longevialle et al., 2008). The results were analyzed using Applied Biosystems 7500 Fast Real-Time PCR Software v2.3 (Thermo Fisher Scientific). All oligonucleotides used in this study are listed in Supplemental Table 4 and were previously described byde Longevialle et al., 2008;de Francisco Amorim et al., 2018).

In Vitro Transcription of 5-EU-Labeled Spike-In RNAs
To generate 5-EU-labeled LUC RNA, we used a control plasmid provided in the TnT Quick Coupled Transcription/Translation System (Promega). The plasmid was linearized using the restriction enzyme SacI and purified using a DNA Clean and Concentrator Kit (Zymo Research). We used a TranscriptAid T7 High Yield Transcription Kit (Thermo Fisher Scientific) according to the manufacturer's protocol for the Synthesis of Non-Radioactively Labeled RNA Probes. 5-Ethynyl-UTP was purchased from Jena Bioscience. A fraction of the reaction was analyzed using a standard agarose gel. The RNA was purified using an RNA Clean and Concentrator Kit. (Zymo Research).
As a spike-in control for chloroplast RNA turnover analysis, an aadA PCR product was amplified using a reverse oligonucleotide including the T7 promoter and purified using a GeneJET PCR Purification Kit (Thermo Fisher Scientific). The 5-EU-labeled aadA spike-in was transcribed in vitro from the PCR product using T7 RNA Polymerase (Thermo Fisher Scientific). Labeling and purification were performed as described for LUC RNA.

Analysis of Neu-RNA-Seq Libraries
To analyze steady state RNA-seq and Neu-seq libraries, we developed a pipeline consisting of library processing, mapping to the reference genome, and analyzing the read count composition of exonic and intronic features. Mapped reads were collected with featureCounts from the RSubread package version 1.5.2 with the following parameters: -T 6 -F SAF -J -O and gene coordinates or exon and intron coordinates as features (Liao et al., 2013). TPM and log 2 TPM values were calculated from collected feature read counts. We added 0.001 to every value to avoid 0 read counts. Log 2 TPM values of features coming from the Neu-seq and steady state transcriptome libraries were plotted with ggplot version 2.2.1 (Wickham, 2009). Spanning read composition analysis was performed using featureCounts' "detailed read assignment" result files. In these files, reads were collected and grouped into exonic, intronic, and exonic-intronic categories based on their position among features. Average values from triplicate measurements from Neu-seq and steady state transcriptome libraries were compared using a Fisher's exact one-tailed test in R (https://www.R-project.org/).

Analysis of Alternative Splicing
To analyze alternative splicing events, mapped libraries (after removing tRNA and rRNA contamination) were analyzed using rMATS version 3.2.5 with the following parameters: -t single -len 150 (Shen et al., 2014). Events of exon skipping, mutually exclusive exon usage, alternative 59 or 39 splice sites, and intron retentions were extracted based on P-value < 0.05 and false discovery rate of 0.05. Detected alternative splicing events are listed in Supplemental Data Set 2.

Illumina Read Processing of ERIC-Seq Libraries
Sequencing libraries of 5-EU-labeled RNAs 0, 1, 2, 6, 12, and 24 h after the chase (three biological replicates each) were sequenced on an Illumina HiSeq3000 system. In total, 395,194,241 reads were sequenced. Adapter removal, trimming of low-quality bases, and tRNA and rRNA filtering and mapping were performed as described above. In total, 339,411,658 reads remained for further analysis. Supplemental Table 3 lists the number of reads in each library.

Analysis of ERIC-Seq Libraries
Mapped reads were collected with featureCounts from the RSubread package version 1.5.2 with the following parameters: -T 6 -F SAF -J -O (Liao et al., 2013). Differential expression analysis was performed using DESeq2 version 1.16.1, substituting size factors with the estimated size factors from the luciferase gene read counts (Love et al., 2014). All genes that showed a zero-expression value in any of the replicates were removed from downstream analysis. Shrunken log 2 values of the expression values were extracted from the DEseq2 object. Additional normalization of shrunken log 2 values was executed: time point value changes were calculated relative to time point 0 (T0 5 1). If the expression of a gene exceeded 1.25 in any of the later time points, the gene was removed from further analysis. Normalized and filtered values were used for visualization using pheatmap version 1.0.8 and model fitting (Kolde, 2012) and for statistical evaluation (Supplemental Data Set 2).
Based on the five separate clusters obtained by heatmap analysis using K-Means clustering, LOWESS and exponential model decay fitting were performed on the five clusters in R to compare overall average half-life values using the following formulas: lowess(x5X, y5Y, f50.2) and nls(Y ; a*exp(-k * X), time_course_object, start5c(a51, k51)).
The individual half-life values of all genes were calculated using an exponential decay model and a composite model with an exponential background and Voigt peak model. All calculations were performed using the python lmfit package version 0.9.7 (Supplemental Data Sets 5 to 8; Newville et al., 2016). The resulting half-lives and chi 2 values (Hu and Bentler, 1999) were compared and Pearson's/Spearman's correlations were calculated between the two different data sets (Supplemental Data Set 5).
Half-life values determined by ERIC-seq were compared with the published RNA half-life values described by Sorenson et al. (2018) and Narsai et al. (2007). Data from Col-0, referred to as the sov data set in Sorenson et al. (2018), was used for comparisons.

Visualization of RNA-Seq Results
For visualization snapshots of the RNA-seq data, we used the Integrative Genomics Viewer version 2.3.80 (Robinson et al., 2011;Thorvaldsdóttir et al., 2013).
Supplemental Figure 2. RNA decay among the individual clusters.
Supplemental Figure 3. Comparison of exponential and compound models for describing ERIC-seq data.
Supplemental Figure 4. Metabolic labeling for RNA stability measurements of chloroplast RNAs.
Supplemental Figure 5. Quality assessment of the ERIC-seq data set.
Supplemental Figure 6. Recovery of the LUC spike-in in the ERICseq experiment.
Supplemental Table 2. Neu-seq library statistics table.
Supplemental Table 3. ERIC-seq library statistics table.
Supplemental Data Set 1. Transcript abundance (in TPM) in steady state transcriptomes and Neu-seq transcriptomes.
Supplemental Data Set 2. Alternative splicing events detected in Neu-seq libraries.

Supplemental Data Set 3. ERIC-seq analysis.
Supplemental Data Set 4. Genes associated with clusters A-E.
Supplemental Data Set 5. Calculation of half-lives using an exponential and compound decay model.
Supplemental Data Set 6. Independent calculation of half-lives in three replicates using an exponential decay model.

ACKNOWLEDGMENTS
Technical support by Irina Passow is gratefully acknowledged. This work was supported by the Deutsche Forschungsgemeinschaft (TRR175-A02 to C.S.-L. and LA2633-4/1 to S.L.) and by the Max Planck Society Chemical Genomics Centre II through its supporting companies Astra-Zeneca, Bayer CropScience, Bayer Healthcare, Boehringer-Ingelheim, and Merck (to S.L.). We thank all members of the lab for critical reading of the article.