Modeling of mRNA deadenylation rates reveal a complex relationship between mRNA deadenylation and decay

Abstract

Complete cytoplasmic polyadenosine tail (polyA-tail) deadenylation is thought to be essential for initiating mRNA decapping and subsequent degradation. To investigate this prevalent model, we conducted direct RNA sequencing of S. cerevisiae mRNAs derived from chase experiments under steady-state and stress condition. Subsequently, we developed a numerical model based on a modified gamma distribution function, which estimated the transcriptomic deadenylation rate at 10 A/min. A simplified independent method, based on the delineation of quantile polyA-tail values, showed a correlation between the decay and deadenylation rates of individual mRNAs, which appeared consistent within functional transcript groups and associated with codon optimality. Notably, these rates varied during the stress response. Detailed analysis of ribosomal protein-coding mRNAs (RPG mRNAs), constituting 40% of the transcriptome, singled out this transcript group. While deadenylation and decay of RPG mRNAs accelerated under heat stress, their degradation could proceed even when deadenylation was blocked, depending entirely on ongoing nuclear export. Our findings support the general primary function of deadenylation in dictating the onset of decapping, while also demonstrating complex relations between these processes.

Keywords: Ccr4-NOT and Pan2/3 Deadenylases; Dcp2 Decapping and Xrn1 Degradation; ONT Nanopore Direct RNA Sequencing (DRS); Pab1; mRNA Deadenylation and Degradation.

Figure 1. mRNA nuclear export reveals cytoplasmic deadenylation and decay dynamics.

(A) Schematic showing the theoretical experimental set-up. (B) Overview of collected time points and their relevance for decay or deadenylation modeling. (C) Time-dependent abundance change of HHF1 or RPL18B mRNAs following Mex67 depletion in DRS datasets (dot) or by reverse transcription coupled with qPCR (ring). A fitted line represents the decay factor calculated from DRS data. (D) mRNA and ncRNA half-life distribution calculated from the Mex67-depletion time course. (E) Comparison of mRNA half-life estimations from various studies visualized as a matrix of Spearman rho coefficients (Miller et al, ; Neymotin et al, ; Presnyak et al, ; Chan et al, 2018). (F) Correlation between half-life values in Mex67-depleted sample and half-life values from Miller et al, . Source data are available online for this figure.

Figure 2. mRNA groups exhibit various deadenylation rates.

(A) Global distribution of mRNA pA-tail lengths during the Mex67-chase experiment. DRS data are presented in a density plot (left) or violin plot (right). Replicates, shown separately in Fig. EV1E, were merged. The number of transcripts (pA-tail estimates) in each density plot is given on the panel. As in Fig. 1B, the violin plot also highlights in gray the timepoints used for deadenylation rate modeling. The latter timepoints were discarded due to the occurrence of mRNA hyperadenylation; a phenotype specific to very few newly made mRNAs in nuclear export-block conditions (Jensen et al, 2001). These scarce species only significantly impact the overall pA-tail distribution after most cytoplasmic mRNA have been deadenylated and degraded. (B) Autoradiogram depicting the global distribution of pA-tail lengths in selected Mex67-chase samples used to construct the sequencing library. (C) Log2 mRNA abundance compared to mean pA-tail length for control Mex67-AA cells, with highly abundant and RPG mRNAs as gold and blue dots, respectively. A bar plot summarizes the read fraction for each transcript category. (D) Global distribution of pA-tail lengths of low-abundance mRNAs. The number of transcripts in each density plot is given on the panel. (E) Time-dependent changes in median pA-tail length across the entire coding transcriptome and four large mRNA groups: high-low abundance, and RPG-non-RPG. Local regression trend lines are shown with a 95% confidence interval. The number of pA-tail estimates used to calculate the median for each point is given in (A, D) and Fig. EV1E,H,J. Replicates were merged to produce nine median estimates for various time points. Source data are available online for this figure.

Figure 3. Pab1 controls deadenylation and decapping.

(A) Global distribution of pA-tail lengths of RNAs isolated from control, Dcp2- or Xrn1-depleted cells. Gray areas delineate the pA-tail lengths underestimated in the DRS library, as previously reported (Tudek et al, 2021). (B) pA-tail length distribution of HCH1 and YMR122W-A mRNAs in control and Dcp2- or Xrn1-depleted cells. (C) pA-tail length distribution of DED81, GAS1, HHF1, HTA1, MCH5, and DCP2 mRNAs in control cells compared to strains depleted of Pab1 using the AID system for 1 or 2 h. The mRNAs are ranked by half-life time. The number of reads contributing to each distribution is indicated in the panel with n = ‘control’/‘1 h depletion’/‘2 h depletion’. The gray boxes show the range of adenosines added de novo in the nucleus by the polyA-polymerase Pap1 on pre-mRNAs (Turtola et al, 2021). (D) Comparison of log2 fold change in mRNA abundance to absolute change in mean pA-tail length for 2 h Pab1-depleted cells compared to control. (E) Schematic illustrating that one Pab1 can bind from 20 to 30 adenosines, therefore newly synthesized mRNA pA-tail bear di- or trimers of Pab1. Depending on the number of Pab1 attached to the pA-tail, the mRNA is either susceptible to deadenylation of decapping. Source data are available online for this figure.

Figure 4. The modified gamma model predicts mRNA pA-tail length distributions.

(A) The probability of protection against decapping in function of pA-tail length is represented by the red line. Gray dots show experimental pA-tail length distributions from all Mex67-chase time points. The number of reads summing up to form the distribution is indicated in Fig. 2A. The presumed number of Pab1 subunits on a pA-tail of a given length is displayed. (B) Conceptual scheme illustrating the modified gamma distribution parameters (see main text). γ_rate corresponds to the exponential coefficient of the right arm of the pA-tail distribution, whereas γ_shape + 1 represents the number of critical events leading to complete mRNA decay. (C) Using gray dots, three density plots show experimental pA-tail length distribution for whole mRNA transcriptome in control cells for replicate A and in control samples for RPL4A and HHF1 mRNAs. Black lines represent fitted modified gamma distributions. Estimation of new mRNA production (by log-fitting) is shown using blue lines. The red line represents the sum of the modified gamma and transcription estimate distributions. the n number indicates the number od pA-tail estimations (reads) used to produce the experimental distribution. (D, E). Changes in the value of γ_shape (D) and γ_rate (E) parameters of the modified gamma probability distribution over time for all mRNAs or those of low and high abundance separately (RPG and non-RPG are omitted for clarity). The parameters are given as full-colored dots supplemented with vertical standard error bars. Each estimate was derived from distributions shown in Figs. 2A, D and EV1A, H–J (refer to those panels for the number of reads). Continuous functions, which were fitted are shown with the equations given. The number of dataset points considered for the estimation of the functions’ equation is indicated on each panel. Those functions were used in subsequent modeling (F, G and Fig. EV2F–J). (F, G) Time-dependent evolution of whole-transcriptome pA-tail length distribution (E) or median (G) based on functions described in (D, E). (H) Modified gamma model of whole-transcriptome deadenylation over time in three-dimension shown as a continuous surface colored according to population density. Black dots show experimental pA-tail distributions. The number of reads summing up to form the distribution is indicated in Fig. 2A. (I) As red triangles with vertical error bars are plotted terminal adenosine half-lives calculated independently for each position relative to their rank in the pA-tail, calculated from distributions adjusted for mRNA decapping (Fig. EV3B). A continuous exponential function is fitted with a 0.95 confidence interval. (J) Experimental pA-tail distributions (same as in (A)) displayed on a normal logarithmic scale (gray dots). The linear distribution between pA-tail lengths 25–60 (yellow box) is highlighted. The slope of the function represents the correction factor used to calculate microscopic enzymatic deadenylation speed from the apparent one. (K) Table summarizing the steps of the deadenylation modeling performed using the modified gamma distribution. Source data are available online for this figure.

Figure 5. At steady-state, deadenylation and decay rates correlate.

(A) Time-dependent change in quantile pA-tail lengths (as color-coded by quantile dots: upper—75–80–85–90–95^th; median—50^th; and lower—15–10–5^th) for the entire coding transcriptome distribution in each replicate. Fitted continuous lines represent the quantile deadenylation coefficient calculated based on the upper quantiles. (B) Distribution of adenosine half-lives relative to pA-tail length in each quantile. The terminal adenosine half-life of 39 + /− 7 s was calculated by averaging the half-lives obtained from quantiles 75–95^th (blue area). (C) Table summarizes key differences between the modified gamma and quantile methods for calculating deadenylation rate. (D) Histogram displaying the distribution of single mRNA terminal adenosine half-lives for the entire coding transcriptome (light-red bars) or mRNAs of high abundance (gray bars). Vertical lines indicate the group median. (E, F) Comparison of decay and quantile deadenylation rates (E) or log2-scaled mRNA half-life to terminal adenosine half-life (F). Density plot represents all transcripts, whereas mRNAs with at least 70 reads in replicate A are shown as individual dots with a blue regression line. Spearman’s rho correlations and p-values were calculated separately for each set using the rstatix package in R. (G) Schematic highlighting the strong, potentially causal link between deadenylation and decapping inferred from data presented in (E, F). (H) Comparison of decay and quantile deadenylation rates for various gene ontology groups (see also Fig. EV4D). (I) Percentage of optimal codons compared to: log2 mRNA abundance, mean mRNA pA-tail length, log2 decay, or quantile deadenylation rates. Source data are available online for this figure.

Figure 6. Deadenylation stimulates decay but is not a prerequisite for RPG mRNA decapping under stress conditions.

(A) Time points collected for the thiolutin 25 °C and heat stress 37 °C chase experiments. Colored surfaces designate the utility of each time frame. Dot contours represent data points excluded from the analysis. (B, C) Distribution of mRNA half-lives (B) or terminal adenosine half-lives (C) in Mex67 depletion, thiolutin, and heat stress chase datasets for RPG and non-RPG mRNAs transcriptionally down-regulated during heat stress according to Vinayachandran et al (2018). The number of estimates in each group is indicated next to the figure legend. (D) Time-dependent change in levels of RPG mRNAs in Mex67 depletion, thiolutin, and heat stress time course normalized to control. Fitted lines show the group trend with a 95% confidence interval in gray. Data points normalized to corresponding controls obtained for 13 min heat stress of a double ccr4Δ pan2Δ mutant (black triangle) or Mex67-AA strain treated or not with rapamycin 3 min prior to heat stress for 13 min at 37 °C. For these data, the median is shown along with the data-point density contour. The experiment legend is listed next to panel 6E. (E) Time-dependent changes in pA-tail length values of the 95–90–85–80–75^th quantiles normalized to the control sample for the same samples as in (D). Normalization allows direct comparison of all changes within the group without grouping transcripts by control pA-tail length. (F, G) Northern blot showing RPS5 levels during a Mex67-depletion time course (F) or heat stress at 37 °C and thiolutin treatment at 25 °C (G). (H) Normalized pA-tail density plot showing the whole-transcriptome adenylation profiles in wild-type cells compared to single ccr4Δ or pan2Δ mutants and a double ccr4Δ pan2Δ strain. The gray box shows the range of adenylation produced de novo in the nucleus by the polyA-polymerase Pap1 on pre-mRNAs (Turtola et al, 2021). (I) Spot tests comparing wild-type cells to ccr4Δ, pan2Δ, and double ccr4Δ pan2Δ mutant at indicated temperatures. (J) Scheme illustrating the role of deadenylation in mRNA decapping and decay. For most mRNAs, deadenylation is a rate-limiting factor that dictates the onset of decapping. RPG transcripts are a special group for which deadenylation can accelerate decay in conditions such as heat stress, but ultimately, decapping and decay are activated by an unknown factor linked to nuclear export. Source data are available online for this figure.

Figure EV1. Following nuclear decay of pre-mRNAs in export-block conditions the cytoplasmic deadenylation and decay of mRNAs is revealed.

(A) Time-dependent change in mRNA abundance relative to control samples in Mex67-depleted cells shown in the form of a series of boxplots. The boxplot central line is the median. The box edges are the 25th and 75th percentiles. Whiskers extend to 1.5 times the IQR (Inter Quartile Range). All three replicates were combined into time point ranges. For each time-slot the number of transcripts assessed in each replicate is given (3875, 3989 and 3032). Sum of 10896 mRNAs is assessed in each time-slot. Note that those are largely overlapping). (B) Time-dependent abundance of selected mRNAs in DRS datasets from Mex67-depleted cells. Fitted lines demonstrate the decay factor calculated for each mRNA based on DRS data. (C) Scatterplot comparing half-life times in Mex67-depleted DRS datasets to those reported by Chan et al (2018). (D) Scatterplot comparing half-life times reported by Miller et al (2011) to those by Chan et al (2018). (E) Global distributions of pA-tail lengths of total mRNA in three replicates of Mex67-AA chase experiments used to model deadenylation. The number of transcripts (pA-tail estimates) in each density plot is given on the panel. (F) Scatterplot comparing single mRNA log2 abundance to half-life. Abundant mRNAs and RPG mRNAs are highlighted with gold and blue dots. (G) Scatterplot comparing mean pA-tail length to mRNA half-life. (H) Global distribution of pA-tail lengths of highly abundant mRNAs isolated from control samples and at various times of Mex67 depletion, presented as a density plot (left panel) and violin plot (right panel). Replicates were merged. The number of transcripts (pA-tail estimates) in each density plot is given on the panel. (I) Global distribution of pA-tail length of low-abundance mRNAs. Replicates were merged. (J) Global pA-tail distribution of RPG and non-RPG mRNAs. Replicates were merged. The number of transcripts (pA-tail estimates) in each density plot is given on the panel. (K) Global pA-tail distribution of individual RPL36A, GAS1, RPL4A, and HHF1 mRNAs. Stars mark long pA-tailed mRNAs accumulating in later Mex67-depletion time points, indicative of hyperadenylation (Jensen et al, 2001). The number of transcripts (pA-tail estimates) in each density plot is given on the panel.

Figure EV2. A modified gamma distribution accurately describes experimental yeast pA-tail distributions.

(A) Plot displaying modified gamma distributions for γ_shape and γ_rate parameters within experimental values. The green line represents a function that effectively describes the pA-tail distribution of the entire transcriptome in control cells (γ_shape = 4; γ_rate = 0.05). (B) Fitting of the modified gamma distribution to experimental datasets. Figure complements examples shown in Fig. 4C. (C) Comparison of the experimental whole-transcriptome mean pA-tail length to the mean pA-tail length of the fitted modified gamma model. A regression line with a 95% confidence interval was fitted to the data points for comparison with the diagonal. The modified gamma distribution slightly overestimates the mean due to discrepancies in estimating very long pA-tails, as seen in Fig. EV2B for pA-tails of 50 and greater. (D) Graph comparing the experimental whole-transcriptome pA-tail length variance length to the variance of the fitted modified gamma model. A regression line with a 95% confidence interval was fitted. (E) Graph comparing changes in the γ_shape and γ_rate parameters of the modified gamma distribution fitted to Mex67-depletion chase time points for all mRNAs and those of high and low abundance. The graph complements plots shown in Fig. 4D, E. (F, G) Graphs showing time-dependent changes in the value of γ_shape (F) and γ_rate (G) parameters of the modified gamma probability distribution fitted into experimental data for selected mRNAs: GAS1, HHF1, RPL36A, and RPL4A. The parameters are given as full-colored dots supplemented with vertical standard error bars. Each estimate was derived from distributions shown in Fig. EV1K (refer to the panel for the number of reads). Equations and continuous lines describe functions fitted to modified gamma parameters to predict the evolution of theoretical distributions (Fig. EV2J). (H) Density plots showing time-dependent evolution of modeled pA-tail length distributions for mRNAs of low and high abundance. The evolution of the median of those distributions is shown in Fig. 4G. (I) Graphs comparing the time-dependent evolution of the modeled modified gamma distribution mean and variance to the experimental one for the whole coding transcriptome. Lines were fitted to both datasets and are displayed with a 95% confidence interval. (J) Density plots showing predicted GAS1, HHF1, RPL36A, and RPL4A pA-tail length distributions obtained by modeling changes in modified gamma parameters (γ_shape and γ_rate), as shown in Fig. EV3F,G. These should be compared to experimental distributions shown in Fig. EV1K.

Figure EV3. The deadenylation process can be reconstituted in silico using the modified gamma distribution.

(A) Graphs illustrating the time-dependent evolution of the γ_shape and γ_rate parameters of the 3-dimensional modified gamma distributions presented in Fig. 4H. (B) Graph depicting time-dependent change in pA⁺ RNA fraction recovery from total using beads coated with oligo-dT₍₂₅₎. A linear function was fitted to predict changes in whole-transcriptome pA⁺ RNA levels and was used to normalize pA-tail levels for adenosine half-life calculations shown in Fig. EV3C. (C) Series of graphs showing absolute change in levels of adenosines at specific pA-tail positions for the whole transcriptome. Half-lives of adenosines at each position were calculated and plotted in Fig. 4I to derive a transcriptomic apparent deadenylation rate. (D) Graph comparing experimental pA-tail density plots for the coding transcriptome in Mex67-depletion chase replicate A (black lines) to distributions resulting from consecutive deadenylation simulation cycles (series of blue-green lines). The constant dictating each deadenylation simulation interval was defined by the dimensionless parameter α (see “Methods”). This was carried out to ensure marked differences between consecutive distributions. All distributions were normalized to the distribution peak (distribution maximum). The in silico distribution that best overlapped with the experimental distributions was selected, assigning a specific value of α. (E) Graph comparing experimental deadenylation times [min] with the ordinal number of the dimensionless α parameter (dimensionless scaling factor) describing the number of in silico deadenylation steps. (F) Graph displaying the time-dependent decrease in pA-tail length values in each quantile normalized to control. The dots designating the upper quantiles (orange dots) cluster around a steeper slope than those designating lower quantiles (blue and maroon dots). (G) Graph showing the number of adenosines removed per minute as a function of pA-tail length. This graph complements the one shown in Fig. 5B.

Figure EV4. A simplified method based on changes to distribution quantile values can measure deadenylation rates.

(A) Series of graphs depicting as color-coded dots the changes in pA-tail length in selected quantiles (upper—75–80–85–90–95^th, median—50^th, and lower—15–10–5^th) over time for selected single mRNA examples. The continuous lines represent the quantile deadenylation coefficients calculated from the upper quantiles. (B, C) Scatterplots comparing terminal adenosine half-life to log2 mRNA abundance (B) or mean mRNA pA-tail length (C). Spearman’s rho correlations and p-values were calculated separately for each set using the rstatix package in R. (D) Comparison of decay to deadenylation rates for various gene ontology groups, which are also shown in Fig. 5H. (E) Correlation matrix comparing decay, deadenylation rates, mRNA abundance, and mean pA-tail length derived from the Mex67-depletion time course to estimates of protein abundance from Ho et al (2018), various translation rate parameters (Siwiak and Zielenkiewicz (2010)), and percent optimal or non-optimal codons.

Figure EV5. RPG and non-RPG mRNAs have altered deadenylation and decay rates during response to heat stress.

(A) Series of graphs showing, in rows, the time-dependent and normalized to control changes in: (top) mRNA abundance by DRS, and (bottom) pA-tail length in upper quantiles (75–95^th) for single RPG mRNAs: NSR1, YEF 3, NOP3, and EGD1 in Mex67 depletion, heat stress 37 °C, and thiolutin 25 °C chase sequencing data. (B) Mean pA-tail lengths of individual mRNAs in wild-type cells were subtracted from the mean lengths observed in deadenylase mutants (ccr4Δ or pan2Δ, or pop2Δ) and displayed in the form of a density plot. The mean change in each strain is marked as a vertical line and specified in the figure legend. (C) Scatterplots showing the absolute change in mean pA-tail length in ccr4Δ - WT on the x-axis compared to pop2Δ – WT (top graph) or pan2Δ – WT (bottom graph) on the y-axis. Transcripts of high abundance or RPG mRNAs are highlighted in blue and red, respectively. Spearman’s rho coefficients and the number of mRNAs compared (n) are listed in each panel legend. (D) Plots displaying the raw pA-tail length density distribution of two example upper quantiles in the wild-type strain heat stress chase compared to the control and 13 min heat shock of a double pan2Δ ccr4Δ mutant strain. The plots present raw data also shown in Fig. 6E. The fitted lines are for orientation purposes only and connect the local maxima either with local or linear regression. (E) Abundance of selected RPG mRNAs normalized to 25S rRNA in control or ski2Δ, pop2Δ, ccr4Δ, pan2Δ, lsm1Δ, and dhh1Δ cells at 25 °C compared to 13 min heat shock at 37 °C determined using reverse transcription coupled to qPCR. Single dots show biological replicate values used to calculate the mean. (F) Abundance of selected RPG mRNAs normalized to 25S rRNA in a double Dcp2-AID and Mex67-AA strain under steady-state and after 13 min heat stress determined using reverse transcription coupled to qPCR. Prior to heat stress the strains were either treated with auxin for 2 h to deplete Dcp2 or with rapamycin for 3 min to deplete Mex67 prior (25 °C). Heat stress was conducted at 37 °C for 13 min. Single dots represent biological replicates. (G) Graphs showing for RPL4A (top) mRNA abundance by DRS, and (bottom) change in pA-tail length in upper quantiles (75–95th) normalized to control for Mex67 depletion, heat stress 37 °C, and thiolutin 25 °C chase sequencing data, along with point heat stress of ccr4Δ pan2Δ double mutant and cells Mex67-depleted 3 min prior to heat stress. (H) A theoretical pA-tail distribution of median 75As and low variance (red line) subjected to deadenylation simulation (series of blue and green lines). Artificial deadenylation was performed with and without inducing decapping, displayed with or without normalization to each distribution maximum value. Normalization to the maximum is shown to better display the systematic change in distribution variance and location of the maximum.