Global SLAM-seq for accurate mRNA decay determination and identification of NMD targets.

Gene expression analysis requires accurate measurements of global RNA degradation rates, earlier problematic with methods disruptive to cell physiology. Recently, metabolic RNA labeling emerged as an efficient and minimally invasive technique applied in mammalian cells. Here, we have adapted SH-linked alkylation for the metabolic sequencing of RNA (SLAM-seq) for a global mRNA stability study in yeast using 4-thiouracil pulse-chase labeling. We assign high-confidence half-life estimates for 67.5% of expressed ORFs, and measure a median half-life of 9.4 min. For mRNAs where half-life estimates exist in the literature, their ranking order was in good agreement with previous data, indicating that SLAM-seq efficiently classifies stable and unstable transcripts. We then leveraged our yeast protocol to identify targets of the nonsense-mediated decay (NMD) pathway by measuring the change in RNA half-lives, instead of steady-state RNA level changes. With SLAM-seq, we assign 580 transcripts as putative NMD targets, based on their measured half-lives in wild-type and upf3Δ mutants. We find 225 novel targets, and observe a strong agreement with previous reports of NMD targets, 61.2% of our candidates being identified in previous studies. This indicates that SLAM-seq is a simpler and more economic method for global quantification of mRNA half-lives. Our adaptation for yeast yielded global quantitative measures of the NMD effect on transcript half-lives, high correlation with RNA half-lives measured previously with more technically challenging protocols, and identification of novel NMD regulated transcripts that escaped prior detection.

INTRODUCTION

Post-transcriptional controls make major quantitative contributions to regulation of gene expression. These include regulation of the decay and translation rates of an mRNA species, as well as its splicing and intracellular location. The level of an RNA species is determined by its rates of synthesis and degradation. While control of transcriptional initiation has been extensively studied, understanding of RNA decay is lagging behind. Turnover rates of RNA molecules are intimately linked to other post-transcriptional processes, such as translation (Pelechano et al. 2015; Chan et al. 2018; Hanson and Coller 2018), intracellular localization (Bovaird et al. 2018), and sequestering in RNA granules (Sheth and Parker 2003; Huch et al. 2016; Escalante and Gasch 2021). Established methods to quantitate the decay rate of RNA species in vivo have limitations that hamper progress in this regard. The stability of individual mRNAs can be reliably quantitated in low throughput by placing them under control of a regulatable promoter, and measuring half-life after transcriptional shut-off (Baudrimont et al. 2017). For global analyses, RNA synthesis can be arrested by inhibition of RNA polymerases, and the degradation rate of all RNA species monitored alternatively metabolic labeling approaches can be used. The budding yeast Saccharomyces cerevisiae has been at the forefront of developing technologies to study most aspects of post-transcriptional regulation, by virtue of its genetic tractability and ease of performing genome-wide studies. In S. cerevisiae, several global studies of mRNA decay rates have been performed. Beyond the basal rates, the effects on global mRNA degradation by environmental stress (Molin et al. 2009; Romero-Santacreu et al. 2009; Miller et al. 2018) or defects in RNA decay pathways (Guan et al. 2006; Celik et al. 2017) have been studied. In S. cerevisiae, a heat-sensitive allele of the largest subunit of RNA polymerase II, rpb1-1, can be used to block mRNA synthesis at the restrictive temperature (Nonet et al. 1987). The downside of this approach is that temperature upshift activates heat shock responses and changes the physiology of the cell. Chemical RNA polymerase II inhibitors, 1,10-phenanthroline, thiolutin, and 6-azauracil, have likewise been used for this purpose. Again, complications arise since each of them selectively induces transcription of specific gene groups (Grigull et al. 2004; Eshleman et al. 2020), confounding their analysis. Both for chemical and genetic blockade of RNA synthesis, perturbations of the physiological state of the cell are thus problematic. These weaknesses of the methods are enhanced when studying transient phenomena such as changes in mRNA stability under environmental perturbations (Rodríguez-Gabriel et al. 2006; Molin et al. 2009; Romero-Santacreu et al. 2009). First, the impact of the heat shock or the chemical inhibitor themselves will obscure that from the environmental perturbation; second, the resolution in time is poor and comparable to the mRNA half-life itself and to the time scale of its changes on environmental shocks.

Attempting to avoid perturbations from the method of study itself, decay rates for RNA species can be indirectly calculated from data on their synthesis rates and steady-state levels, for example, from genomic run-on experiments using in vivo radioactive labeling (Marín-Navarro et al. 2011; Jordán-Pla et al. 2019). More recently, metabolic labeling of RNA using nucleoside analogs such as 4-thiouridine (4sU) has been used to study the transcriptional process (Russo et al. 2017). Alternatively, other nucleoside analogs including 5-bromouridine (Tani et al. 2012) or 5-ethynyluridine (Abe et al. 2012) can be used, nevertheless studies in mammalian cells have relied mostly on 4sU. This nucleoside analog cannot be imported into unmodified yeast cells, unless a specific transporter is expressed to facilitate the uptake of 4sU as was done in Miller et al. (2011) using human equilibriative nucleoside transporter (hENT1). On the other hand, the enzyme uracil phosphoribosyltransferase will metabolize the base analog 4-thiouracil (4tU) into the nucleotide analog 4-thiouridine monophosphate, allowing its incorporation into nascent RNA in yeast, but not in mammalian cells. Consequently, 4tU is instead the preferred labeling choice for budding or fission yeast since it alleviates the need for a transporter and can be used with unmodified cells.

The thiolation of 4tU enables it to engage in reversible covalent conjugation to thiol-containing molecules such as streptavidin. This can be used for biochemical enrichment using solid phase immobilization of thiolated mRNAs. With pulse-chase labeling protocols, this allows detailed kinetic studies of transcription and RNA degradation. Such protocols have been used for global measurements of transcript stability in yeast (Chan et al. 2018). In that work, the mean half-life of S. cerevisiae mRNA was estimated to be under 5 min, considerably lower than in previous studies. A mean value of 23 min was measured using temperature shift with rpb1-1 (Wang et al. 2002), and a mean of 35.3 min was calculated from experiments where tagged Rbp1 was depleted from the cell nucleus (Geisberg et al. 2014). One possible reason for such large discrepancies between half-life estimates using inhibition of RNA synthesis versus using metabolic labeling is that RNA polymerase inhibition may also affect RNA degradation. It was observed that application of thiolutin to yeast cells increased apparent mRNA half-life in a dose-dependent manner (Pelechano and Pérez-Ortín 2008). It has been proposed that the processes of mRNA synthesis and degradation promote each other, and that mRNA synthesis inhibition leads to reduced mRNA decay and vice versa (Haimovich et al. 2013). In a further development of metabolic labeling (“SLAM-seq”), after extraction of RNA, the incorporated thiolated uracil derivatives are instead alkylated with iodoacetamide (Herzog et al. 2017). In the subsequent reverse transcription reaction, guanine will base pair with the iodoacetamide uracil derivative, resulting in a T-to-C conversion in the position originally occupied by 4tU. This eliminates errors and time consumption from the biochemical purification steps.

Nonsense-mediated decay (NMD) was originally discovered in S. cerevisiae through the accelerated degradation of mRNAs with premature stop codons (PTCs) (Losson and Lacroute 1979). The sources of PTCs can be transcriptional error, defective splicing, or translational shifting. It has subsequently emerged that NMD can be triggered by other structures than PTCs, and affects a large proportion of mRNAs. Thus, long 3′-untranslated regions (3′-UTRs) can also bring about NMD (Hogg and Goff 2010; Hurt et al. 2013). It has been calculated that up to 20% of all eukaryotic RNAs are affected directly or indirectly by NMD (Mendell et al. 2004; Guan et al. 2006). The core proteins required for NMD are Upf1-3, present throughout eukaryotes and first discovered in S. cerevisiae (Leeds et al. 1992). Among these, the RNA helicase Upf1 is the most highly conserved in evolution. The Upf2 and Upf3 proteins are thought to bridge the effector Upf1 and the structure triggering the NMD reaction, stimulating its helicase activity (Chamieh et al. 2008). It has been proposed that the Upf proteins are also required for processes beyond NMD, for example, telomere maintenance, through incompletely elucidated mechanisms (Askree et al. 2004). On the other hand, these effects may also be mediated through telomere maintenance mRNAs being NMD targets. In metazoans, the SMG1 protein kinase activates Upf1 through phosphorylation; yeast Upf1 also becomes phosphorylated but its kinase is as yet unidentified (Lasalde et al. 2014).

The total population of NMD targets is not fully characterized in any organism. Such a characterization is an extensive task, as many RNA species arise from one transcriptional unit, and these can have different propensity to be targeted by NMD. Moreover, the extent to which a particular RNA species is degraded by NMD may vary between conditions, as NMD is also regulated by for example, stress conditions (Goetz and Wilkinson 2017). By measuring steady-state levels of S. cerevisiae mRNA species in wt and upf mutants, direct and indirect NMD targets were found to comprise around 15% of the transcriptome (Celik et al. 2017).

Here, we applied SLAM-seq for S. cerevisiae using 4tU; previously, this protocol was optimized for mammalian cells using 4sU. We have determined individual mRNA half-lives in S. cerevisiae cells in the BY4741 background and compared wild-type (wt) with the upf3Δ mutant, defective in NMD. We provide the first application of SLAM-seq to study mRNA stability in yeast. Further, we assign a direct quantitative measure of the change of mRNA half-life in an NMD defective background, rather than only observing changes in RNA abundance, which also includes indirect effects through for example, mRNAs encoding transcription factors as NMD targets. SLAM-seq obviates the need for biochemical purification of derivatized mRNA by utilizing high efficiency chemical conversion. We paired SLAM-seq with 3′-UTR end sequencing after having scrutinized the annotation of 3′-UTRs genome-wide to allow cost effective assessment of the stability of polyadenylated mRNA transcripts (Herzog et al. 2020), and capture a total of 3033 transcripts for which we can assign a half-life estimate with high precision in both the wt and the NMD mutant in biological triplicate. We find an average mRNA half-life of 11 min (median 9.4 min) in wt and 13 min (median 11.6 min) in the upf3Δ mutant. We identify transcripts from 580 genes as being under the control of NMD under standard conditions. These include 355 transcripts that have been found in earlier studies, and 225 not previously recognized.

RESULTS

SLAM-seq utilizes a highly efficient derivatization of thiol-containing uridines with iodoacetamide that eventually leads to specific mis-incorporation of a cytosine (C) instead of a thymine (T) (T-to-C conversion) upon reverse transcription (Herzog et al. 2017). The chemical conversion greatly reduces the labor associated with multitime point metabolic labeling experiments by eliminating the need for biochemical enrichment techniques and allowing the use of a lower starting amount for construction of cDNA libraries. Moreover, modeling of RNA decay is simplified due to the elimination of accounting for biotin conjugation and separation of biotinylated RNAs efficiency in the modeling parameters.

SLAM-seq has been utilized in mammalian and Drosophila cells lines using 4sU (Herzog et al. 2017; Muhar et al. 2018; Reichholf et al. 2019); however, no reports exist for the use for SLAM-seq in yeast. We reasoned that SLAM-seq would have a distinct advantage over transcriptional inhibition methods due to its minimally invasive nature, and would be less laborious than currently used pull-down methods. Given that this method can be easily utilized to measure RNA half-life, we measured the RNA decay in the BY4741 wt strain and the isogenic upf3Δ mutant. The aim was to identify NMD substrates since half-life measurements are the golden standard for this as opposed to just measuring changes in steady-state levels in the transcriptome (Fig. 1).

FIGURE 1.

SLAM-seq workflow. Log phase cells are labeled with 4-tU followed by chasing the label with a high concentration of unmodified uracil. 4-tU is taken up by cells and incorporated into RNA. Samples are taken during the chase and metabolically inactivated using chilled ethanol. RNA is then extracted, and incorporated 4-tU is alkylated with iodoacetamide (IAA), yielding a stable derivative of the thiol-containing 4-tU. Reverse transcriptase misincorporates a guanine rather than adenine at derivatized positions during library construction, causing a T to C conversion. These are then quantitated using the SLAMDUNK pipeline.

Confidence of half-life estimates and coverage of RNA species using SLAM-seq in yeast

We started by treating cells in logarithmic growth phase with a low concentration of 4tU for 2 h followed by chasing the labeled transcripts with a high concentration of uracil. A low level of 4tU (0.2 mM) was used with a labeling time of 2 h since it has been shown that concentrations up to 1 mM do not affect growth for at least 3 h in the low-uracil media used here (Chan et al. 2018). Hence, we used the longest applicable labeling time with a safety margin of 1 h to ensure consistent conversion rates in the RNA. The samples were incubated for 2 min prior to the onset of sampling to account for the inefficient chase that has been previously observed (Chan et al. 2018), and time points were sampled at 0, 1, 4, 8, 18, and 58 min in order to reasonably capture both short- and long-lived RNAs. Subsequently RNA was extracted, alkylated, reverse transcribed, and sequenced. We observed a T-to-C conversion rate ranging between 2.54% and 2.83% for BY4741 at time point 0 and between 2.41% and 2.55% in upf3Δ, similar to observations in mouse embryonic stem cells (Herzog et al. 2017) with the conversion rate decreasing with increasing chase time as expected (see Materials and Methods for details and Supplemental Files S1–S4 for quality control reports and sample indices). We started by comparing the performance of SLAM-seq to other methods. We limited the analysis to transcripts with a measureable half-life in all three biological replicates with a coefficient of variation less than 0.35. A full unfiltered data set can be found in Supplemental File S5.

Using all six time points, we were able to fit 4049 transcripts in BY4741 that meet the criteria above, representing 67.5% of all annotated ORFs in the wt or 72.8% of the ORFs available in our 3′-UTR annotation, and 3638 transcripts in upf3Δ (60.7% of all annotated ORFs and 65.4% of ORFs available in our 3′-UTR annotation). The linear distribution of the half-life was right-skewed as expected (Fig. 2A). The wt BY4741 strain showed a sharper peak than upf3Δ, with the individual biological replicates showing high correlation coefficients, ranging between 0.88 and 0.93 in BY4741 and between 0.83 and 0.88 in upf3Δ (Fig. 2B). Gene Ontology (GO) term enrichment analysis of the 1% most long-lived transcripts indicated that those encoding proteins involved in housekeeping functions such as glycolysis and translation had the longest half-life (Fig. 2A; a full enrichment table is in Supplemental File S6), while the 1% with the shortest half-life did not show any particular enrichment.

FIGURE 2.

Half-life distribution and comparison to other methods. (A) Half-life distribution in wt and upf3Δ. The dashed line indicates mean and solid line the median; red is for wt and blue for upf3Δ. Partial enrichment analysis for the longest lived 1% of the transcripts is shown as an inset (full GO term enrichment for longest lived 1% is listed in Supplemental Table S6). (B) Inter-replicate correlation using Pearson's correlation coefficient. (C) Comparison of half-lives calculated using SLAM-seq to other methods using Spearman's rho. GRO (Pelechano et al. 2010), cDTA (Sun et al. 2013), DTA (Miller et al. 2011), MGC (Baudrimont et al. 2017), Chan et al. with and without poly(A) selection (Chan et al. 2018), Munchel et al. (2011), RATE-seq (Neymotin et al. 2014), Anchor away (Geisberg et al. 2014), rpb1-1 (Wang et al. 2002), phenanthroline (Grigull et al. 2004). Chan et al. (2018) RNA total and Munchel et al. (2011) were not significantly correlated (P = 0.996), hence the correlation was left blank. (D) Correlation of half-lives to other transcripts features using Spearman's rho. RNA abundance was calculated using the mean CPM across all the time points of wt, protein abundance was from Ho et al. (2018), codon optimality was from Drummond et al. (2006), and ORF length from annotations of S. cerevisiae S288C (assembly R64).

Correlation of SLAM-seq with other methods for mRNA turnover measurements

Next, we compared SLAM-seq and other mRNA decay measurements using half-lives calculated from the BY4741 strain. Figure 2C shows the Spearman correlation between various studies using transcriptional inhibition, transcriptional shutoff or metabolic labeling. As expected, transcriptional inhibition methods correlated very poorly with each other and with other methods, while multiplexed gene control and variants of metabolic labeling correlated well with each other and with SLAM-seq, the rho values ranging between 0.5 and 0.8. It has been previously suggested that a Spearman rho value of 0.7 and above indicates high intermethod reliability (Wada and Becskei 2017). We observe high correlation (0.78) with multiplexed gene control (MGC), the metabolic labeling study with improved RNA separation (Chan et al. 2018); 0.78, and the dynamic transcriptome analysis (DTA) method (Miller et al. 2011) (0.8), indicating that SLAM-seq is able to classify transcripts as stable or unstable in good agreement with other recent methods. However, the absolute magnitude of measured half-life is different between the methods. Chan et al. (2018) and studies using MGC calculate a shorter half-life than DTA and SLAM-seq (see Discussion). We also compared our results to those of Sun et al. (2013), who used comparative dynamic transcriptome analysis (cDTA) to determine degradation rates in wt and mutants affected in RNA degradation, and noted a high correlation (0.72) with their results. The correlation also remained high (>0.7) for half-lives calculated for upf3Δ using SLAM-seq to the half-lives calculated for upf2Δ and upf3Δ using cDTA (Spearman rho correlation values can be found in Supplemental Fig. S12). Finally, we checked for half-life correlation with other ORF features (Fig. 2D). The results mirrored that of Chan et al. (2018), namely a weak inverse correlation of mRNA half-life with ORF length (rho = −0.26) and a positive correlation with protein abundance (rho = 0.6), RNA abundance (rho = 0.58) and codon optimality (rho = 0.52) all of which indicates that highly expressed transcripts and transcripts with optimal codons tend to have longer half-lives.

Assignment of transcripts as NMD targets

We found 3033 transcripts passing the above filtration parameters that were in common between upf3Δ, and could be directly compared. We found 346 transcripts to have a significantly altered half-life in upf3Δ (q value of <0.05 and half-life fold change of at least 1.5). However, the criteria used to filter the transcripts for comparison was too restrictive for reporting NMD targets, since if two replicates reported a similar half-life in wt but the other replicate did not pass the fitting filter (R2 ≥ 0.6) and the replicates in upf3Δ passed all the other criteria, the transcript would still be excluded from the analysis due to the lack of one half-life measurement. Therefore, we used greater than or equal to twofold half-life fold change as criteria for the other transcripts that were not part of the previously compared 3033 transcripts and had at least two measurements in one strain and three measurements in the other. Combining the above analyses, we found a total of 598 transcripts with altered half-life between the two strains (full list can be found in Supplemental File S7. The distribution of half-life ratios between upf3Δ and wt for all mRNAs that met either of the above criteria is shown in Supplemental Fig. S13). Additional transcripts that had less than 1 count per million (CPM) in any of the time points and hence filtered out in the wt even though a number of them had increased read counts in upf3Δ, were not considered in this set since we could not calculate their half-lives. These transcripts, thus representing weakly expressed mRNAs, were only detectable in upf3Δ mutants. The transcripts that failed the CPM filtration in wt can be found in Supplemental Table S8. The majority of the transcripts with an altered half-life relative to wt had an increased half-life in upf3Δ (n = 580) with only 18 transcripts showing a decreased half-life in upf3Δ and having no particular enrichment. Considering that the half-life of NMD substrates increases rather than decreases, we only considered the transcripts with increased half-life ratio for further analysis.

We were able to recover transcripts that represent the common structural classes regulated by NMD (Fig. 3): (1) nonsense codon-containing YML002W. In our strain background, YML002W is continuous with YML003W, and a single thymidine deletion in YML003W causes a frameshift and a premature stop codon, hence rendering the transcript sensitive to NMD; (2) mRNAs utilizing frameshifting in their translation (EST3); (3) uORF-containing transcripts (CPA1); (4) transcripts with long 3′-UTR (COX19). Moreover, we recover eight additional transcripts that were previously identified in low-throughput studies (Fig. 4). We compared our identified NMD targets with previous studies (Table 1) and found an overall overlap of 61.2%, indicating that we were able to successfully identify NMD substrates based on changes in half-life. The group of targets were analyzed for enrichment of GO terms, and were found to be enriched for two subgroups of interest: protein sumoylation by molecular function, and spliceosome by cellular component (Table 2).

FIGURE 3.

Transcripts representing common structural classes known to be regulated by NMD. Background subtracted and chase onset normalized T-to-C conversion rates were fitted to first order exponential decay. YML002W: nonsense codon-containing in this strain background. EST3: translational frame-shifting. CPA1: uORF-containing. COX19: long 3′-UTR. Error bars represent the standard deviation. Error bars for YML002 in upf3Δ were not added since it was calculated from two replicates.

FIGURE 4.

Transcripts found to be NMD regulated in previous low throughput studies. Background-subtracted and chase onset normalized T-to-C conversion rates were fitted to first order exponential decay. The first row is control transcripts not regulated by NMD. Error bars represent standard deviation. Transcripts with no error bars were calculated from two replicates. Transcripts were found to be NMD regulated in the following studies: FZF1, INO4, PDR3, and PDR8 (Guan et al. 2006); RAD54 and RAD57 (Janke et al. 2016); MPA43 (Kebaara and Atkin 2009); SPO16 (Zaborske et al. 2013).

TABLE 1.

Overlap with other studies aiming to identify NMD targets

TABLE 2.

Enrichment analysis for transcripts with increased half-life in upf3Δ

DISCUSSION

Recent understanding and chemical derivatization of RNA species has led to the emergence of new techniques that simplify the measurements of previously difficult to capture parameters such as RNA-half life on a genome-wide scale. SLAM-seq has been heavily used in mammalian cells (Baptista and Dölken 2018; Matsushima et al. 2018; Muhar et al. 2018); however, no reports for simpler eukaryotic organisms are available. In this work, we conducted an assessment of the applicability of SLAM-seq to S. cerevisiae, and found excellent correlation with recent studies utilizing metabolic labeling. Moreover, we applied SLAM-seq to measure changes in RNA half-lives in NMD defective cells, and were able to identify NMD substrates based on actual RNA half-life change as opposed to increased steady-state RNA levels. Using SLAM-seq, we could avoid the problems associated with other techniques such as transcriptional inhibition (e.g., rpb1-1), which essentially measures RNA half-life in dying cells, hence giving rise to physiologically irrelevant half-life measurements that are poorly reproducible, and also bypassed the requirement for ample starting material and spike-in normalization required for pull-down approaches. Using SLAM-seq, we were able to calculate half-lives with high intersample and intermethod correlation as opposed to the poorly correlating values reported previously in transcriptional inhibitions studies, due to the problems highlighted above.

SLAM-seq was paired with 3′-end mRNA sequencing (QuantSeq). This alleviates the cost of full RNA-seq, especially given the multitimepoint experiments that are required for half-life determination. It also facilitates downstream analysis due to the lack of a requirement for transcript length normalization. The QuantSeq protocol captures fully processed RNA pol II transcripts, which can be advantageous if that is the subset of interest. However, the inability to distinguish RNA isoforms may pose a problem for certain applications. In our case, we were unable to detect transcripts belonging to an NMD class consisting of inefficiently spliced transcripts that leak to the cytoplasm (Celik et al. 2017) with the exception of HRB1. This could indicate that the mature form of HRB1 could also be under NMD regulation. Alternatively, since HRB1 has lower expression than other transcripts belonging to the same class (e.g., YBR089C-A,YFL034C-A, and YGR148C), the increase in unspliced isoforms can be detectable with half-life change since it is able to compete with the mature form for labeling and hence can contribute to half-life calculation. For the other more strongly expressed transcripts, this effect is too low to be seen. Nonetheless, this does not detract from the value of the technique as S. cerevisiae is intron-poor. Further, this problem can be remedied with other types of RNA-seq libraries as long as they incorporate a reverse transcription step (Herzog et al. 2017).

We found a high correlation in the classification of transcripts into stable and unstable transcript with other metabolic labeling studies. However, the absolute value for RNA half-life was not consistent between the studies. This has been an ongoing problem for half-life determination in yeast due to earlier reports indicating a higher half-life value, but considering that earlier studies utilized transcriptional inhibition, it is not surprising to see higher half-life values that poorly correlate with our study and even other studies using derivatives of transcriptional inhibition. In this work we report a mean half-life value of 11 min (median 9.4) which agrees more closely with Miller et al. (2011) in which a mean value of 14 min (median 11 min) was reported, while the other high correlating studies (Baudrimont et al. 2017; Chan et al. 2018) reported much lower values than all previous measurements. These differences might be attributable to methodic differences; an earlier study (Munchel et al. 2011), which had lower separation efficiency of labeled transcripts compared to the second-generation assay (Chan et al. 2018), calculated values that did not correlate well with any of the previous studies. This highlights the technical difficulties associated with pull-down approaches and the effect of the pull-down efficiency on half-life calculation. However, it is noteworthy that our study was the only one besides Munchel et al. (2011) to use a pulse-label approach while all the other studies generally used approach to equilibrium. One rationale behind approach to equilibrium is to counter the recycling of labeled nucleotides which may compromise the chase (Nikolov and Dabeva 1985; Wolfe et al. 2019). However, as we show here, the effect seems to be minimal when using SLAM-seq. Baudrimont et al. (2017) similarly reported a lower half-life mean than previous estimates agreeing more closely to that of Chan et al. (2018). It has been argued that slower than expected kinetics of nucleotide incorporation might be the reason for these discrepancies, since in MGC the target transcript in placed under a modified controllable promoter and decay is measured post-shutoff, and the delay from nucleotide incorporation is eliminated (Wada and Becskei 2017). Regardless of the difference in the absolute values, the high correlation of the methods (hence ability to appropriately group transcripts as stable and unstable) can be appreciated since all those techniques are methodically independent. Moreover, if indeed the delayed incorporation kinetics has an effect on half-life calculation, future correction can be made using a factor calculated based on Baudrimont et al. (2017), given that data for more transcripts becomes available.

We report new NMD substrates based on changes in RNA half-life, rather than steady-state abundance, and were able to identify 225 new transcripts, whereas the remaining 61.2% overlapped with previous studies (Table 1). We find it likely that a large fraction of the 225 newly identified NMD targets in this work, using altered decay rates in upf3Δ mutants as criteria, escaped detection in earlier studies using mRNA steady-state levels. The highest overlap with a single previous study (50%) was with Celik et al. (2017), where RNA-seq was utilized to identify NMD targets in high resolution, indicating that SLAM-seq agrees well with more recent data sets that aimed to identify NMD targets with the added benefit of producing a quantitative value for the half-lives. We successfully capture major classes of transcripts expected to be NMD targets, with the exception of the inefficiently spliced transcripts class (see above). We found an enrichment for spliceosomal factors in agreement with what has been found in mammalian cells (Saltzman et al. 2008). We envision that SLAM-seq would be highly useful for studies looking to identify changes in transcripts stabilities under various conditions such as those looking at the effect of particular drugs on the cells or metabolic changes under different nutritional conditions such as nitrogen sources (Miller et al. 2018), in addition to the application to NMD pathway substrate detection as we describe here.

MATERIALS AND METHODS

Strains and media

BY4741 (MATa his3Δ1 leu2Δ0 met15Δ0 ura3Δ0) was used as the control strain (WT), and an isogenic upf3Δ strain was used as the tester. YPD (2% glucose, 2% peptone, and 1% yeast extract) was used for routine culturing of strains. Low uracil medium (LUM; 0.19% yeast nitrogen base, 0.5% ammonium sulfate, 2% glucose and 0.077% Complete Supplement Mixture without uracil [-URA CSM; ForMedium], supplemented with 0.1 mM uracil) was used for labeling conditions. The chase medium used to remove the RNA label was identical to LUM except for being supplemented with 20 mM uracil. Liquid cultures were grown in a rotary shaker at 30°C at 200 rpm.

4tU metabolic labeling

Overnight cultures were made in LUM and diluted to OD600 nm = 0.1 on the next day in a total volume of 200 mL, then grown to OD600 nm = 0.4. Subsequently, the cultures were split into two 100 mL cultures, and 4tU (Sigma) was added to one culture to a final concentration of 0.2 mM, and DMSO was added to the negative control culture to a final concentration of 0.02%. Both cultures were then grown for an additional 2 h. To chase the 4tU label, the cultures were then briefly centrifuged, the 4tU-containing media removed and replaced with an equivalent volume of chase medium and incubated for 2 min at the same growth condition used for labeling. Time 0 was taken after the 2 min incubation with the chase media and at 1, 4, 8, 18, and 58 min after time 0. For each time point, a 10 mL aliquot was taken and added to an equal volume of 95% ethanol prechilled at −70°C, and then kept on dry ice until the end of the sampling. Only one time point was sampled for the DMSO-treated cultures to calculate the background T-to-C conversion rate. This treatment was repeated in biological triplicates for both WT and upf3Δ.

RNA extraction

Total RNA was extracted from the sampled cell pellets using the Quick-RNA Fungal/Bacterial Kit (Zymo Research) as per manufacturer instructions, with the modification of adding reducing agent (RA) from the SLAM-seq Kinetics Kit—Catabolic Kinetics Module (Lexogen) to the isolation, wash, and elution buffers of the RNA extraction kit at a ratio of 1:1000, 1:1000, and 1:100, respectively. Five micrograms of the extracted total RNA was subjected to alkylation with iodoacetamide (IAA) at a final concentration of 10 mM. Each reaction consisted of 5 µL of freshly prepared 100 mM IAA dissolved in 100% ethanol, 25 µL organic solvent (OS), and 5 µL of sodium phosphate (NP) in a final volume of 50 µL (all reagents were from the SLAM-seq Kinetics Kit—Catabolic Kinetics Module [Lexogen]). The reaction was run for 15 min at 50°C, then stopped using 1 µL of the stopping reagent (SR) from the same kit. RNA was repurified by ethanol precipitation according to the SLAM-seq Explorer and Kinetics Kits User Guide (Lexogen). The entire procedure starting from the extraction was done while protecting the RNA from white light by wrapping the tubes in foil and working under red light until alkylation of the samples was complete.

Library preparation

An amount of 200 ng of alkylated total RNA was used as input for generating 3′-end mRNA sequencing libraries using the QuantSeq 3′ mRNA-seq Library Prep Kit FWD for Illumina (Lexogen) and the PCR Add-on Kit for Illumina (Lexogen) according to the QuantSeq 3′ mRNA-seq Library Prep Kit User Guide (Lexogen). Library preparation, sequencing in SR75 mode, and read demultiplexing were carried out by Lexogen Services.

Bioinformatics analysis and half-life analysis

Demultiplexed reads were analyzed using SLAM DUNK (https://t-neumann.github.io/slamdunk/) (Neumann et al. 2019) through the nf core SLAM-seq processing and analysis pipeline (Ewels et al. 2020). Default settings were used with the exception of designating the read length manually. A custom BED file containing the position of the ORFs 3′-UTRs was created based on three annotations (Xu et al. 2009; Roy and Chanfreau 2020) and our own analysis of the TIF-seq data described previously (Pelechano et al. 2013) and used for defining the windows for 3′-UTRs (Supplemental Files S9, S10). The resulting count files were filtered to include genes that had at least 1 CPM in their corresponding time series. Subsequently, T-to-C conversion rates were corrected by subtraction of the background conversion rate of the corresponding DMSO-treated samples and normalized to the onset of the chase (time 0). Half-lives were calculated using the minpack.lm package as described previously (Herzog et al. 2017), and values that had an R2 value of less than 0.6 were filtered out from the resulting calculated half-lives; exact sampling times were used to reduce error associated with slight deviations between theoretical and actual sampling time points, with actual time points being the moment the sample is mixed with the chilled ethanol (see Supplemental File S11 for exact time points). Values of half-lives were corrected for growth during the chase using the half-life formula in Munchel et al. (2011). Transcripts that passed the more stringent criteria n = 3033 (all replicates had values and a coefficient of variation less than 0.35) were tested using an unpaired t-test with a single pooled variance assumption after log10 transformation of the values, and q values were calculated using the two-stage linear step-up procedure described previously (Benjamini et al. 2006). Transcripts with a q value of less than 0.05 and a fold change in a half-life of 1.5 or more were considered positives. For transcripts that were not part of the above set, a greater than or equal to twofold half-life fold change was used as a criterion, and each transcript needed to have at least two measurements in one strain and three measurements in the other. Enrichment analysis was carried out using Yeastmine (https://yeastmine.yeastgenome.org/yeastmine/). Graphing of RNA degradation curves was done using GraphPad Prism 9. The SLAM-seq data sets are available at the Gene Expression Omnibus (GEO) with accession number GSE196690.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.

COMPETING INTEREST STATEMENT

J.A.Z.-M. is an employee of Lexogen GmbH. Lexogen declares no conflicts of interest.