High initiation rates decrease protein output of stall-containing S. cerevisiae mRNAs

To measure the effects of initiation rate variation on protein expression, we used a fluorescent reporter system in the budding yeast S. cerevisiae (Fig 1A). Our reporters consist of the PGK1 gene from S. cerevisiae, which has been extensively used in previous studies of mRNA translation and stability [5, 13]. We fused the PGK1 coding sequence to a 3×FLAG-encoding sequence at the 5′ end and the yellow fluorescent protein (YFP) gene at the 3′ end. We integrated this 2-kb protein-coding region, under the control of the GPD promoter, the GPD 5′ UTR, and the CYC1 3′ UTR, as a single copy into the genome (see Strain and plasmid construction). We expressed a red fluorescent protein (mKate2) gene with identical regulatory sequences as our reporter from a different genomic locus for normalizing the measured YFP fluorescence.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. High initiation rates decrease protein output of stall-containing S. cerevisiae mRNAs.

(A) Schematic of 3×FLAG-PGK1*-YFP reporters used in (C) and (D). The −4 to −1 nt region preceding the ATG start codon had one of eight different Kozak sequences as indicated. One of three different stall sequences or their respective controls were inserted at the end of PGK1*. PGK1* was modified from wild-type PGK1 sequence based on a previous study [11] (see Strain and plasmid construction). (B) Protein levels of 3×FLAG-PGK1-YFP control reporters with no insert and wild-type PGK1 sequence. (C) Protein levels of the constructs shown in (A). Protein levels are quantified as the mean YFP fluorescence of 10,000 cells for each strain as measured using flow cytometry. Protein levels are expressed as a.u. relative to the mean RFP levels from a constitutively expressed mKate2 control. Error bars in (B) and (C) show standard error of the mean over three or four independent transformants. Many error bars are smaller than data markers. Two of the total 192 strains were clear outliers and were removed after manual inspection. (D) Western blots of reporters with low (CTGC), medium (CCAC), or high (CAAA) initiation rates and with indicated stall sequences or controls using antibody against the FLAG epitope at the N terminus. Histone H3 level is shown as loading control. Numbers for each lane indicate the ratio of the FLAG signal to the H3 signal and are normalized to a maximum of 10 within each blot. The underlying data for panels B and C can be found at https://github.com/rasilab/ribosome_collisions_yeast. a.u., arbitrary unit; chrI, chromosome I; nt, nucleotide; YFP, yellow fluorescent protein.

https://doi.org/10.1371/journal.pbio.3000396.g001

To confirm our reporter system's utility for studying both normal and stall-containing mRNAs, we introduced five tandem arginine AGA codons (5×AGA) at five different locations along the PGK1 gene (henceforth referred to as PGK1*, S1 Fig panel A), following earlier work [11]. We individually mutated each of the five 5×AGA repeats to their synonymous CGG rare arginine codons and measured their effect on YFP expression using flow cytometry. The 5×CGG rare codons reduce YFP expression by 40%–70%, with the stalls away from the ATG start codon having a stronger effect (S1 Fig panel B).

To vary the initiation rate of our reporters, we designed eight different 5′ UTR variants in which we mutated the 4 nucleotides preceding the ATG start codon. These mutations vary 3×FLAG-PGK1-YFP expression over a 5-fold range (Fig 1B). This hierarchy of YFP expression among the 5′ UTR variants is concordant with previous measurements [38] (S1 Fig panel C). Because the 5′ UTR mutations are located over 30 nucleotides away from the 5′ end of the mRNA, we expect them to have minimal effects on the binding of translation initiation factors and mRNA decay factors to the 5′ cap.

We then combined the 5′ UTR mutations with various stall and control sequences to measure their combinatorial effect on protein expression (Fig 1A). We used the 5×CGG stall sequence from above, as well as 8×CCG rare proline codon repeats and 10×AAG lysine codon repeats, which are known to trigger ribosome stalling [39–41]. As controls, we used 5×AGA, 8×CCA, and 10×AGA repeats, respectively.

Strikingly, the three 5′ UTR variants (CCAA, CAAA, AAAA) that have up to 2-fold-higher protein level than the median for normal mRNAs (Fig 1B) instead cause up to a 4-fold reduction in protein level when the mRNA contains any of the three ribosome stalls (orange triangles, Fig 1C). This decrease in protein expression at high initiation rates is largely absent when the stalls are replaced by control sequences of synonymous common codons (5×AGA and 8×CCA, Fig 1C).

Similar to our above results with rare arginine and proline codon-mediated stalls, replacing the 10×AAG lysine repeat stall by a 10×AGA arginine control abrogates the decrease in protein expression at high initiation rates (Fig 1C, left panel). Importantly, the 10×AAG repeat differs from the 10×AGA control only by a single-nucleotide frameshift, and it is located over 1.3 kb away from the initiation codon. Hence, changes in mRNA secondary structure or other subtle physical interactions between the 5′ UTR mutations and the stall mutations are unlikely to account for the dramatically different behaviors of the 10×AAG and 10×AGA inserts at high initiation rates. Furthermore, because AAG is a common lysine codon with abundant cellular tRNA [42, 43] and our reporter is integrated as a single copy in the S. cerevisiae genome, the observed effects on protein expression at high initiation rates are unlikely to arise from titration of rare tRNAs by the stall-encoding repeats.

Western blotting against the N terminus 3×FLAG epitope qualitatively reproduced the differences in protein expression as measured by flow cytometry (Fig 1D). These results confirm that high initiation rates have distinct effects on expression of the full-length protein depending on the presence of ribosome stall sequences in the corresponding mRNA.

High initiation rates decrease lifetime of stall-containing S. cerevisiae mRNAs

We reasoned that the change in steady-state protein expression with initiation rate in our experiments could arise from change in stability of mRNAs, change in synthesis rate of the full-length protein from each mRNA, or a combination of both. To isolate the effect of initiation rate variation on mRNA stability, we developed a high-throughput sequencing–based assay to quantify mRNA levels of our PGK1-YFP reporters at high resolution (Fig 2A). In our assay, each combination of 5′ UTR mutation and codon repeat is tagged with four distinct 8-nucleotide barcodes in their 3′ UTR. We pooled these barcoded reporters and integrated them into the genome of S. cerevisiae. Reporter mRNA from the pooled strains was reverse transcribed to cDNA. The barcode region from both cDNA and genomic DNA was then amplified by PCR and counted by sequencing. The relative mRNA level of each reporter variant was obtained by normalizing the cDNA count of the corresponding barcode by its genomic DNA count.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. High initiation rates decrease stability of stall-containing S. cerevisiae mRNAs.

(A) Schematic of deep-sequencing assay used for measuring mRNA levels. Reporters with different 5′ UTR mutations and stall or control inserts were tagged with four unique 8-nt barcodes in their 3′ UTR. The reporters were pooled and transformed into S. cerevisiae. cDNA and gDNA from the pooled reporters were amplified, and the barcodes were counted by high-throughput sequencing. Barcodes that had fewer than 100 counts for either cDNA or gDNA were discarded during analysis. mRNA level was calculated as cDNA counts normalized by gDNA counts for each barcode and was median normalized across all reporters in each pool. (B) mRNA levels of 3×FLAG-PGK1-YFP control reporters with no insert and wild-type PGK1 sequence. The leftmost point is for a reporter with the ATG start codon mutated to CTG. (C) mRNA levels of 3×FLAG-PGK1*-YFP constructs with stall or control inserts. Error bars in (B) and (C) show standard error over 3 or 4 distinct barcodes for each reporter variant. One reporter variant (10×AGA, ACGC) had only two barcodes with read counts higher than the threshold for discarding. Most error bars are smaller than data markers. The underlying data for panels B and C can be found at https://github.com/rasilab/ribosome_collisions_yeast. a.u., arbitrary unit; gDNA, genomic DNA; nt, nucleotide; YFP, yellow fluorescent protein.

https://doi.org/10.1371/journal.pbio.3000396.g002

We tested our experimental strategy by measuring steady-state mRNA levels of the 3×FLAG-PGK1-YFP reporters with no ribosome stall sequences and eight different 5′ UTR mutations. mRNA levels differ less than 25% between the eight 5′ UTR variants (Fig 2B), even though protein levels differ more than 5-fold between these variants (Fig 1B). Because steady-state mRNA levels reflect the balance between mRNA synthesis and decay, we parsimoniously interpret these measurements to imply that our 5′ UTR mutations change only the translation initiation rate at the ATG start codon of our reporters without affecting their mRNA stability or transcription rate. As a positive control, mutating the ATG start codon of our reporter to CTG decreases mRNA levels by 4-fold (Fig 2B). We expect this decrease if absence of the main ATG start codon results in translation initiation at one of the four downstream out-of-frame start codons in the 3×FLAG-encoding region, which then destabilizes the mRNA by nonsense-mediated mRNA decay.

We then applied our experimental strategy to measure mRNA levels of the 5′ UTR variants of our 3×FLAG-PGK1*-YFP reporters with various stall or control inserts (Fig 2C). The 8×CCG and 10×AAG stall sequences result in up to a 2-fold decrease in mRNA levels, specifically in the high initiation rate regime (yellow triangles, Fig 2C). By contrast, 3×FLAG-PGK1*-YFP reporters with the 8×CCA and 10×AGA control inserts show little or no decrease in mRNA levels at high initiation rates (gray circles, Fig 2C), which is similar to the 3×FLAG-PGK1-YFP constructs without inserts (Fig 2B). We interpret the decrease in mRNA levels of the 8×CCG and 10×AAG stall-containing reporters at high initiation rates as changes in their mRNA lifetime. This interpretation is justified given that the 5′ UTR mutations, on their own, do not affect mRNA levels (Fig 2B).

Comparison between the protein and mRNA measurements from our PGK1-YFP reporters show that protein levels at high initiation rate decrease up to 4-fold from their peak value, whereas mRNA levels decrease only 2-fold (CCGA versus AAAA 5′ UTR mutations in Figs 1C and 2C). Thus, we conclude that the decrease in protein expression at high initiation rates arises from both steady-state change in mRNA levels and change in the synthesis rate of the full-length protein from each mRNA.

Collision-stimulated abortive termination model predicts reduced protein output at high initiation rates

We hypothesized that computational modeling can provide mechanistic insight into how high initiation can result in lower protein expression and reduced stability of stall-containing eukaryotic mRNAs. Because kinetic models that consider both translation and stability of eukaryotic mRNAs have not been formulated so far to our knowledge, we first defined a joint model for normal translation and canonical mRNA decay in eukaryotes (see Kinetic modeling). We modeled normal translation as a three-step process composed of initiation, elongation, and termination, similar to previous studies [44, 45]. Each ribosome occupies a footprint of 10 codons on mRNAs [46], and initiation and elongation occur only if the required footprint is not blocked by another leading ribosome. We modeled canonical mRNA decay in eukaryotes as a three-step process of deadenylation, decapping, and exonucleolysis [47]. Even though a more detailed kinetic model of canonical mRNA decay is available [47], because we aimed to predict overall mRNA lifetime, we considered only the three main steps of canonical mRNA decay. We chose the rate constants for these decay steps [47] (S3 Table) to result in a mean mRNA lifetime of 35 min, corresponding to that of the PGK1 mRNA in S. cerevisiae [5]. We used a rule-based modeling approach and performed exact stochastic simulations using an agent-based framework to predict mRNA lifetime and protein expression (see Kinetic modeling).

We then systematically investigated the effect of quality control at elongation stalls on protein expression. We first studied kinetic models of how stalled ribosomes undergo abortive termination by setting the mRNA decay rate to zero (S3 Table). In our model, abortive termination includes all molecular processes that terminate polypeptide synthesis including splitting of the stalled ribosome, degradation of the nascent peptide, and hydrolysis of the peptidyl-tRNA. We did not distinguish between these different molecular events because all of them reduce synthesis rate of the full-length protein, which is the quantity that we measure in our experiments.

We considered three kinetic models of how abortive termination might occur at stalled ribosomes from intact mRNAs (Fig 3A, Kinetic modeling). In the simple abortive termination (SAT) model, stalled ribosomes are subject to kinetic competition between abortive termination and normal elongation (Fig 3A). This SAT model has been considered in recent modeling studies [35, 48, 49]. Because recent experiments found a role for ribosome collisions in quality control [26, 36, 37], we then considered two possible effects of ribosome collisions on abortive termination. In the collision-stimulated abortive termination (CSAT) model, only stalled ribosomes that get hit by trailing ribosomes undergo abortive termination (Fig 3A). Conversely, in the collide and abortive termination (CAT) model, only trailing ribosomes that run into leading stalled ribosomes undergo abortive termination (Fig 3A). Finally, as a control, we also considered the TJ model, in which ribosomes do not undergo abortive termination (Fig 3A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Collision-stimulated abortive termination model predicts reduced protein output at high initiation rates.

(A) Schematic of different abortive termination models simulated in (B). k_a denotes the nonzero rate of abortive termination of ribosomes from indicated configurations. In the TJ model, ribosomes do not abort translation. In the SAT model, only ribosomes that have not undergone collision from the mRNA entry side abort. In the CAT model, only ribosomes that have undergone collision from the mRNA entry side (“trailing” ribosomes in a collision) abort. In the CSAT model, only ribosomes that have undergone collision from the mRNA exit side (“leading” ribosomes in a collision) abort. (B) Average protein synthesis rate as a function of initiation rate predicted using stochastic simulations of the four models in (A). The simulations were of a 650-codon mRNA corresponding to the 3×FLAG-PGK1*-YFP reporters in our experiments. Ribosome stalls were simulated as a consecutive stretch of six slowly translated codons with a net elongation rate of 0.1 s⁻¹. All other codons had an elongation rate of 10 s⁻¹. The stalls were located after 400 codons, corresponding to their approximate location in our experiments. (C) Effect of varying the number of tandem (consecutive) slowly translated codons encoding the stall in the CSAT model. The elongation rate of each slowly translated codon was scaled to maintain a net elongation rate of 0.1 s⁻¹ across the stall. mRNA decay rate was set to zero in these simulations to isolate the effect of abortive termination. All other model parameters are listed in S3 Table. Standard errors from repeating stochastic simulations with different initial random seeds are smaller than data markers in (B) and (C). The underlying data for panels B and C can be found at https://github.com/rasilab/ribosome_collisions_yeast. CAT, collide and abortive termination; CSAT, collision-stimulated abortive termination; SAT, simple abortive termination; TJ, traffic jam; YFP, yellow fluorescent protein.

https://doi.org/10.1371/journal.pbio.3000396.g003

We simulated the above-described kinetic models of abortive termination with varying initiation rate and analyzed the predicted protein synthesis rate (defined as the number of full proteins produced per second). We simulated elongation stalls by introducing a stretch of six tandem poorly translated codons after 400 codons in a 650-codon mRNA, which is similar in length to the PGK1-YFP reporters that we use in our experiments. In our simulations, the control TJ model without abortive termination exhibits a linear increase in protein synthesis rate with initiation rate until the elongation rate at the stall becomes rate limiting (gray triangles, Fig 3B; S2 Fig panel A). This observation is consistent with previous studies of this model [50, 51]. The SAT model exhibits a similar behavior to the TJ model, but with lower protein synthesis rate and a higher initiation rate at which protein synthesis rate reaches saturation (yellow triangles, Fig 3B, S2 Fig panel B, C).

Unexpectedly, the two kinetic models in which ribosome collisions cause abortive termination of either the trailing or the leading ribosome (Fig 3A) exhibit different behaviors as initiation rate is varied. In the CAT model in which the trailing ribosome abortively terminates upon collision, protein synthesis rate increases monotonically as the initiation rate is increased (green diamonds, Fig 3B). By contrast, in the CSAT model, in which the leading ribosome abortively terminates upon collision, protein synthesis rate increases initially, then reaches a maximum, and surprisingly decreases with further increase in initiation rate (blue circles, Fig 3B).

The above dichotomy between the CSAT and CAT models of abortive termination at ribosome collisions can be intuitively understood as follows: In the CAT model, even though the trailing ribosomes at collisions undergo abortive termination, every ribosome that reaches the stall eventually finishes protein synthesis at a rate determined by the elongation rate at the stall. In the CSAT model, the leading stalled ribosome is stimulated to undergo abortive termination upon collision; thus, at sufficiently high initiation rates and abortive termination rates, very few ribosomes will get past the stall. Stated differently, the difference between the CSAT and CAT models arises because all ribosomes that finish protein synthesis stall at the slow codons, whereas not all ribosomes collide with a leading ribosome before they finish protein synthesis.

The decrease in protein expression at high initiation rates in the CSAT model is counterintuitive, and to further probe its origin, we systematically varied the elongation rate at the stall and the number of tandem (consecutive) codons encoding the stall. The decrease in protein synthesis rate requires the initiation rate to exceed the total elongation rate past the stall (S2 Fig panel D) because ribosome collisions occur only in this regime. More surprisingly, the decrease in protein synthesis occurs only if ribosomes stall at multiple consecutive codons (Fig 3C). As a concrete example, the CSAT model predicts that two tandem stall codons with a 5 s ribosome dwell time at each codon decreases protein synthesis rate at high initiation rates, but a single stall codon with a 10 s ribosome dwell time does not. This nonintuitive prediction can be understood by observing that the number of kinetic partitions between normal elongation and abortive termination changes dynamically with initiation rate in the CSAT model with n consecutive stall codons. At low initiation rates (k_init≪k_stall), there is no kinetic partitioning toward abortive termination, similar to the TJ model. At high initiation rates (k_init≫k_stall), each ribosome that finishes protein synthesis undergoes n kinetic partitions toward normal elongation, similar to the SAT model. However, at intermediate initiation rates (k_init~k_stall), the average number of kinetic partitions increases from 0 to n as the initiation rate is increased. In this sense, the CSAT model “interpolates” between the TJ model and the SAT model.

Thus, our simulations of different kinetic models of abortive termination reveal distinct signatures of initiation rate variation on protein expression. The CSAT model, unlike the other models, predicts a decrease in protein expression at high initiation rate. Importantly, this decrease is not simply a consequence of ribosome collisions stimulating quality control, but it has two other essential ingredients in our simulations (Fig 3B and 3C): the leading ribosome in the collision undergoes abortive termination, and the stall itself is composed of multiple kinetic steps.

Collision-stimulated endonucleolytic cleavage model predicts decrease in mRNA lifetime at high initiation rates

Because endonucleolytic cleavage has been shown to reduce mRNA lifetime of stall-containing mRNAs [15], we next considered different kinetic models of how endonucleolytic mRNA cleavage at ribosome stalls decreases mRNA lifetime as a function of initiation rate (Fig 4A, Kinetic modeling). We modeled endonucleolytic mRNA cleavage as occurring on the 5′ side of the stalled ribosome [52, 53]. In our modeling, once mRNAs are endonucleolytically cleaved, translation initiation is immediately repressed through decapping, and ribosomes that have already initiated on the 5′ mRNA fragment are efficiently recycled once they reach the truncated 3′ end [19, 20]. We simulated ribosome stalling by introducing a stretch of six poorly translated codons after 400 codons in a 650-codon mRNA, which is similar to our PGK1-YFP reporters. To isolate the effect of endonucleolytic cleavage on our predictions, we set the abortive termination rate from intact mRNAs to be zero. We monitored the mean lifetime of mRNAs, which we define as the time interval between the end of transcription and the start of 5′–3′ exonucleolytic decay after decapping. As before, we also monitored protein synthesis rate as we varied the initiation rate of mRNAs in our simulations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. CSEC model predicts decrease in mRNA lifetime at high initiation rates.

(A) Schematic of SEC and CSEC mRNA cleavage models simulated in (B) and (C). (B, C) Average mRNA lifetime and protein synthesis rate over 10⁶ s as a function of initiation rate predicted using stochastic simulations of the models in (A). Reporters were simulated as in Fig 3, but with a nonzero endonucleolytic cleavage rate of 0.001 s⁻¹ and a net elongation rate of 0.1 s⁻¹ across the stall. Canonical mRNA decay was allowed, whereas abortive termination rate was set to zero. All other model parameters are listed in S3 Table. Standard errors from repeating stochastic simulations with different initial random seeds are smaller than data markers in (B) and (C). The underlying data for panels B and C can be found at https://github.com/rasilab/ribosome_collisions_yeast. CSEC, collision-stimulated endonucleolytic cleavage; SEC, simple endonucleolytic cleavage.

https://doi.org/10.1371/journal.pbio.3000396.g004

We considered two distinct kinetic models of how endonucleolytic mRNA cleavage might occur at ribosome stalls (Fig 4A, Kinetic modeling). In the simple endonucleolytic cleavage (SEC) model, cleavage occurs through simple kinetic competition with normal elongation [54]. In the collision-stimulated endonucleolytic cleavage (CSEC) model, cleavage occurs only between two collided ribosomes.

Both the simple and collision-stimulated models of endonucleolytic cleavage predict a decrease in mRNA lifetime at high initiation rates (Fig 4B). However, whereas the SEC model predicts a gradual decrease in mRNA lifetime (yellow triangles, Fig 4B), the CSEC model predicts a sharp decrease in mRNA lifetime as the initiation rate is increased (blue circles, Fig 4B). The dependence of mRNA lifetime on initiation rate in the SEC model is surprising at first glance because there is no direct role for initiation rate in this model. However, even in the SEC model, the endonucleolytic cleavage probability of a given mRNA increases with the frequency of ribosome stalling, which in turn increases with the frequency of initiation. Thus, in both our kinetic models of endonucleolytic cleavage, mRNAs transition from being primarily degraded through the canonical decay pathway at low initiation rates to being endonucleolytically cleaved at high initiation rates.

Even though mRNA lifetimes decrease at high initiation rates in both the simple and collision-stimulated cleavage models, protein synthesis rates in these models display strikingly different behaviors as a function of initiation rate (Fig 4C). In the SEC model, protein synthesis rate increases monotonically with initiation rate (yellow triangles, Fig 4C): at low initiation rates, it increases linearly where mRNAs are degraded primarily through the canonical decay pathway, and it saturates at high initiation rates at a value determined by the endonucleolytic cleavage rate at stalls (S3 Fig panel A). In the high initiation rate regime of the SEC model, each mRNA is translated by a fixed number of ribosomes (on average) before it undergoes endonucleolytic cleavage. By contrast, in the CSEC model, protein synthesis rate exhibits a nonmonotonic behavior (blue circles, Fig 4C): it increases linearly until the initiation rate matches the elongation rate at the stall, at which point it decreases sharply and saturates at the same value as the SEC model. Unlike the CSAT model (blue circles, Fig 3C), the behavior of the CSEC model does not depend sensitively on the number of tandem slowly translated codons that the ribosome transits through at each stall (S3 Fig panel B).

Thus, our simulations of different kinetic models of endonucleolytic mRNA cleavage reveal distinct signatures of initiation rate variation on mRNA lifetime and protein expression. The CSEC model, unlike the SEC model, predicts the decrease in mRNA stability observed in our experiments only at high initiation rates (Fig 2C). The CSEC model also correctly predicts the observed decrease in protein expression at high initiation rates (Fig 1C), which is independent of the predicted effect of abortive termination on protein expression in the CSAT model.

Hel2 and Asc1 attenuate gene expression from stall-containing mRNAs only at high initiation rates

We then sought to test the prediction from our modeling that collision-stimulated quality control (CSEC and CSAT) directly contributes to the decrease in protein expression and mRNA stability at high initiation rates from stall-containing mRNAs. Toward this, we focused on several factors that have been previously implicated in quality control at ribosome stalls but whose effect on gene expression as a function of initiation rate has not been characterized. The ribosome-ubiquitin E3 ligase Hel2 (ZNF598 in humans) is activated by collided ribosomes [26, 36, 37]. The ribosome-associated protein Asc1 (RACK1 in humans) is located at the interface between collided ribosomes [36, 37] and, along with Hel2, couples translational arrest to nascent chain degradation [55–58]. We also consider the ribosome recycling factor Dom34 (PELO in humans), which plays a critical role in no-go decay [15, 54, 59]. Nascent chain degradation is mediated by the E3 ligase Ltn1 (LTN1 in humans), which is part of the RQC complex [16, 18, 21, 22, 60]. We measured protein expression and mRNA levels from our 3×FLAG-PGK1*-YFP reporters in S. cerevisiae strains in which HEL2, ASC1, DOM34, and LTN1 were individually deleted (Fig 5A–5C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Hel2 and Asc1 attenuate gene expression from stall-containing mRNAs only at high initiation rates.

(A) Protein levels of 3×FLAG-PGK1*-YFP reporters (see Fig 1A) with low (CTGC), medium (CCAC), or high (CAAA) initiation rates and with stall (8×CCG) or control (8×CCA) repeats in different deletion backgrounds. (B) Western blots of reporters from (A) performed and quantified as in Fig 1. (C) mRNA levels of 3×FLAG-PGK1*-YFP reporters. Mean mRNA levels and their standard errors quantified as in Fig 2. The ΔLTN1-8×CCA-AAAA variant has only two barcodes. (D) Protein levels of the 8×CCG stall reporter expressed in ΔHEL2 (top) or ΔASC1 (bottom) strain and complemented with the indicated Hel2 or Asc1 variant. Error bars in (A) and (D) calculated as in Fig 1. ** in (D) denotes P < 0.01 (two-tailed t test) for comparison between indicated Hel2/Asc1 variant and the parent strain (“None”) for the same 5′ UTR sequence. The underlying data for panels A, C, and D can be found at https://github.com/rasilab/ribosome_collisions_yeast. a.u., arbitrary unit; nt, nucleotide; WT, wild-type; YFP, yellow fluorescent protein.

https://doi.org/10.1371/journal.pbio.3000396.g005

Deletion of ASC1 completely rescues protein expression from the 8×CCG stall-containing reporters (ΔASC1 in Fig 5A). In fact, the 8×CCG reporter has slightly higher expression than the 8×CCA control reporter at high initiation rate (bottom versus top panel, ΔASC1 in Fig 5A). Similarly, deletion of HEL2 increases protein expression of stall-containing reporters by 3-fold at high initiation rate compared with the wild-type (WT) strains but has no effect at low initiation rate or in the control reporters (ΔHEL2 in Fig 5A). Deletion of either DOM34 or LTN1 has little to no effect on protein expression from either the stall-containing or control reporters (ΔDOM34 and ΔLTN1 in Fig 5A). However, deletion of LTN1, the E3 ligase necessary for nascent chain degradation, stabilizes the peptide resulting from abortive termination at the 8×CCG stalls only at high initiation rates (S4 Fig panel C). Western blotting against the N terminus 3×FLAG epitope qualitatively matches the above measurements of protein expression using flow cytometry (Fig 5B). The rescue of protein expression at high initiation rates in the ASC1 and HEL2 deletion strains is also observed with the weaker 5×CGG stalls (S4 Fig panel A).

We then used our pooled sequencing–based strategy to quantify the mRNA levels of our reporters (Fig 2A) in different deletion backgrounds. Deletion of either ASC1 or HEL2 selectively rescues the decreased mRNA levels of the 8×CCG stall-containing reporters at high initiation rates (ΔASC1 and ΔHEL2 in Fig 5C). By contrast, deletion of neither DOM34 nor LTN1 has any effect on the decreased mRNA levels at high initiation rates (ΔDOM34 and ΔLTN1 in Fig 5C).

We checked whether the observed rescue of protein expression in the ASC1 and HEL2 deletion strains can be reversed by constitutive expression of the respective proteins in trans from the HIS3 locus (Fig 5D). Complementing by the WT Asc1 protein fully reverses the protein expression rescue at high initiation rates from 8×CCG stall-containing mRNAs in ΔASC1 strains (WT, Fig 5D, lower panel). This reversal is abrogated either partially or completely with Asc1 mutants that are defective in translational arrest [61] (16HNG18AAA, 38RDK40AAA, and 85WDK87AAA, Fig 5D, lower panel). Similarly, complementing by the WT Hel2 protein partially reverses the protein expression rescue at high initiation rates from stall-containing mRNAs in ΔHEL2 strains (WT, Fig 5D, upper panel). This reversal is not observed with a Hel2 mutant that cannot bind the interacting E2 enzyme Ubc4 [58] (ΔRING, Fig 5D, upper panel). In contrast to the 8×CCG stall reporters, complementing with various Asc1 and Hel2 mutants had little effect on protein expression from the 8×CCA-containing control reporters (S4 Fig panel B).

Based on the above measurements, we conclude that Asc1 and Hel2, which have been previously associated with ribosome collision–stimulated quality control [26, 36, 37], are necessary to reduce protein expression and mRNA stability at high initiation rates from stall-containing mRNAs. By contrast, neither Dom34 nor Ltn1 have a role in regulating gene expression at high initiation rates.

Endogenous mRNAs with stall sequences show signatures of inefficient translation initiation

Stall-inducing polybasic tracts decrease protein expression and mRNA stability at high initiation rates in our reporter-based experiments; do they also shape the translation of endogenous S. cerevisiae mRNAs? S. cerevisiae coding sequences that contain stretches of 6 or more lysine and arginine codons within a 10-codon window lead to ribosome stalling as measured using ribosome profiling [16]. Similarly, tandem repeats of 2 or more prolines also induce ribosome pausing [39, 40]. Over 1,250 S. cerevisiae protein-coding sequences contain either 10-codon stretches with at least 6 lysine and arginine codons or 10-codon stretches with at least 6 proline codons. Using a high-quality S. cerevisiae ribosome profiling dataset [62], we recapitulated previous observations [16] of increased ribosome density about 8 codons into the stall (Fig 6A). We did not observe a similar increase around control sequences that are enriched for glutamate or aspartate codons (S5 Fig panel A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Endogenous mRNAs with stall sequences show signatures of inefficient translation initiation.

(A) Mean ribosome density around stall-encoding sequences in S. cerevisiae mRNAs using data from Weinberg and colleagues [62]. Stall-encoding sequences are defined as 10-codon windows that have either a minimum of 6 lysine and arginine codons or a minimum of 6 proline codons. The first nucleotide of the 10-codon window is at distance 1 nt. In total, 1,251 S. cerevisiae mRNAs have at least one stall-encoding sequence. Arrow indicates peak in ribosome density at +24 nt that is consistent with Brandman and colleagues [16]. (B) Effect of endogenous stalls on protein expression from reporter mRNAs. Protein levels measured by flow cytometry and their mean and standard error quantified as in Fig 1. ***, *, N.S. indicate P < 0.001, P < 0.05, P > 0.05 (one-tailed t test) for higher ΔASC1/WT protein-level ratio of the stall-encoding reporters in comparison with the no-stall control. (C) TE of S. cerevisiae mRNAs that either contain or do not contain stall-encoding sequences. TE is defined as the normalized ratio of ribosome footprint counts to total mRNA counts as measured by Weinberg and colleagues [62]. (D) Fold change in mRNA levels [57] between both ΔHEL2 and ΔASC1 strains and WT strain. (C, D) Box plots shows mean and standard deviation within each gene group; violin plot shows the gene density at each y-axis value; P values calculated using the two-sided Wilcoxon rank-sum test. The underlying data for panels A, B, C, and D can be found at https://github.com/rasilab/ribosome_collisions_yeast. a.u., arbitrary unit; nt, nucleotide; N.S., not significant; TE, translation efficiency; WT, wild-type.

https://doi.org/10.1371/journal.pbio.3000396.g006

To test whether ribosome stall sequences in endogenous genes can cause collision-stimulated decrease in gene expression, we cloned the 10-codon stall-encoding sequences from six endogenous mRNAs into our PGK1-YFP reporter system with an optimal initiation context (Fig 6B, left panel). Two of the six tested endogenous stall sequences have a significantly lower protein expression compared with a no-stall control (Fig 6B, right panel, yellow circles). However, this comparison does not account for the different amino acid sequences among the stall-containing reporters, which could alter the measured fluorescence. Because Asc1 is required for collision-stimulated decrease in gene expression (Fig 5A), we expect that protein expression from stall-containing reporters should increase upon deletion of ASC1. Consistent with this expectation, five of the six tested endogenous stall sequences have a significantly higher ratio of protein expression between the ΔASC1 and WT strains when compared to the no-stall control (Fig 6B, right panel, yellow versus blue circles).

The widespread nature of stall-encoding sequences in endogenous S. cerevisiae mRNAs suggests two distinct possibilities in light of the collision-stimulated decrease in gene expression uncovered in our work. First, an endogenous mRNA containing a stall could have a sufficiently high initiation rate such that it is constantly subject to collision-stimulated quality control. Indeed such a regulatory mechanism was proposed for RQC1, which encodes one of the components of the RQC complex [16]. Alternatively, a stall might be present because of functional requirements on the protein that are unrelated to collision-stimulated quality control. In this case, the presence of the stall might necessitate that the corresponding mRNA have sufficiently low initiation rate to prevent collision-stimulated decrease in mRNA stability and protein expression.

To distinguish between the above two possibilities, we tested whether stall-containing mRNAs have significant differences in initiation rate compared with other endogenous mRNAs in S. cerevisiae. We used TE as quantified by the ratio of ribosome footprints to mRNA abundance as a proxy for initiation rate [62]. We found that the TE of S. cerevisiae mRNAs containing stall-encoding sequences (as defined above, N = 1,030) is significantly lower than that of mRNAs (N = 3,477) that do not contain such sequences (P < 10⁻¹⁵, Wilcoxon 2-sided test, Δlog₂ TE = 0.23, Fig 6C). One potential confounding factor in this analysis is that abortive termination at stalls can selectively decrease the apparent TE of stall-containing mRNAs. However, we observed a similar difference in TE between stall-containing and control mRNAs, even when TE was calculated using only mRNA regions 5′ to stalls (S5 Fig panel B). Another possible confounding factor is that the stall-encoding sequences cause an increase in local ribosome occupancy on the corresponding mRNA, thus possibly affecting the measured TE. However, this increase will result in a higher apparent TE, thereby lowering the measured difference between stall-containing and control mRNAs.

Based on our observation that stall-encoding S. cerevisiae mRNAs have lower TE, we hypothesize that inefficient translation initiation on these mRNAs allows them to avoid collision-stimulated reduction in protein expression and mRNA stability. This hypothesis predicts that deletion of HEL2 or ASC1 should have no effect on the levels of these endogenous mRNAs, unlike the increase in levels of our reporter mRNAs at high initiation rates (Fig 5C). To verify this prediction, we compared the levels of S. cerevisiae mRNAs between ΔHEL2 or ΔASC1 strains and WT strains, as measured previously by RNA sequencing (RNA-seq) [57]. Because deletion of these factors can have indirect and pleiotropic effects unrelated to quality control [63], we looked for consistent differential regulation of mRNA levels in both ΔHEL2 and ΔASC1 strains compared with the WT strain. We find that, on average, stall-containing S. cerevisiae mRNAs do not show a significantly higher up-regulation compared with control mRNAs in the ΔHEL2 and ΔASC1 strains (P > 0.7, Wilcoxon 2-sided test, Fig 6D). This observation suggests that stall-containing S. cerevisiae mRNAs have been evolutionarily selected to avoid collision-stimulated reduction in gene expression through a decrease in their translation initiation rates. An alternate possibility is that the stalls are not sufficiently strong to trigger quality control on many of the endogenous mRNAs. The latter possibility, however, does not explain why these group of stall-containing endogenous mRNAs should have significantly lower TE.

Strain and plasmid construction

S1 and S2 Tables list the plasmid cloning vectors and the S. cerevisiae strains used in this work. DNA_ sequences.fasta (https://github.com/rasilab/ribosome_collisions_yeast/blob/master/data/dna_sequences.fasta) contains the sequences of the reference plasmid vectors and inserts used in this work. Plasmids and strains will be sent upon request. PCR oligo sequences are available upon request.

The BY4741 S. cerevisiae strain background (S288C, MATa HIS3Δ1 LEU2Δ0 MET15Δ0 URA3Δ0, Thermo Fisher) was used for all experiments in this study. The pHPSC16 vector containing mKate2 and the LEU2 marker was integrated as a single copy at the HO locus. The resulting strain (scHP15) was used as a parent for inserting all PGK1-YFP reporters and for deleting quality control–associated genes. The p41894 vector and its derivatives containing PGK1-YFP variants and the URA3 marker were inserted as a single copy at a ChrI intergenic locus. The pSB2273 vector containing HEL2 or ASC1 variants and the HIS3 marker were inserted as a single copy at the HIS3 locus. Integration into the S. cerevisiae genome was performed by transforming 0.5–2 μg of the linearized (NotI or PmeI restriction) plasmid vector according to the LiAc/SS carrier DNA/PEG method [73]. Single yeast colonies were selected on synthetic complete (SC) agar plates lacking either leucine, uracil, or histidine.

LTN1, DOM34, HEL2, or ASC1 was deleted from the S. cerevisiae genome using PCR-mediated gene disruption [74]. Deletion cassettes with flanking homology arms are provided in DNA_sequences.fasta. Primers for amplifying KAN or NAT resistance markers with 40–50-bp homology arms were designed using the sequences from the S. cerevisiae genome deletion project webpage (http://www-sequence.stanford.edu/group/yeast_deletion_project/downloads.html).

All plasmids were cloned by restricting the parent vectors and inserting 2–4 PCR-amplified fragments with 20–30-bp homology arms using isothermal assembly [75]. The inserted sequences were confirmed by Sanger sequencing. The 5′ UTR mutations for varying initiation rate and the ribosome stall/control repeats were introduced into the PCR primers used for amplifying inserts prior to isothermal assembly. These 5′ UTR and insert sequences are provided in DNA_sequences.fasta. Standard molecular biology procedures [76] were followed for all other steps of plasmid cloning. For cloning the endogenous stall inserts in Fig 6B, the last lysine codon in PGK1 was removed so that it did not artificially amplify the observed effect of the endogenous stalls on protein expression.

The 3×FLAG-PGK1-YFP reporter genes were inserted into the p41894 vector [77] between the GPD (TDH3) promoter + 5′ UTR and the CYC1 3′ UTR + transcriptional terminator using the XbaI and XhoI restriction sites. The sequence around the XbaI site of p41894 was modified to retain the native 5′ UTR sequence of TDH3 while mutating the –12 to –7 nucleotide from ATG to XbaI (TAAACA to TCTAGA). PGK1* protein-coding sequence was generated from WT S. cerevisiae PGK1 by introducing 5×AGA repeats at the five locations used in a previous study [11]. The S. cerevisiae-optimized YFP (mCitrineV163A variant) sequence was taken from a previous work [78] and was fused with PGK1 or PGK1* with an intermediate BamHI site. The sequence of the resulting p41894-3×FLAG-PGK1*-YFP vector, pHPSC57, is provided in DNA_sequences.fasta. The 3×FLAG-PGK1-YFP insert sequence between XbaI and XhoI of pHPSC417 is also provided.

The HO3-LEU2-pGPD-mKate2-Cyc1t vector was constructed by inserting two 400-bp sequences from the HO locus of the S. cerevisiae genome into a pUC19 backbone. The LEU2 gene from pRS315 [79] and a GPDp-mKate2-CYC1t cassette were inserted between the two 400-bp HO sequences. The mKate2 coding sequence was taken from previous work [80] and was synonymously mutated to the most frequent codon for each amino acid in the S. cerevisiae genome. The GPDp-Cyc1t flanking regions are identical to the ones in our reporter vector, described above. The sequence of the resulting vector, pHPSC16, is provided in DNA_sequences.fasta.

The HEL2 and ASC1 WT or mutant coding sequences were inserted into pSB2273 between the GPD promoter and the ADH1 terminator between the XhoI and BamHI restriction sites. The HEL2 coding sequence was PCR-amplified from the S. cerevisiae genome, whereas the ASC1 coding sequence was PCR-amplified from S. cerevisiae cDNA to exclude the intron. Sequence of the parent vector, pSB2273, is provided in DNA_sequences.fasta. All the HEL2 and ASC1 insert sequences cloned into pSB2273 are also provided.

For creating the barcoded reporter strains for mRNA measurements, each of the PGK1-YFP reporter coding sequences was PCR-amplified separately with four distinct barcodes in the 3′ UTR. The ATG → CTG control reporters were amplified with three distinct barcodes. The PCR products were separately inserted into pHPSC57 between the XbaI and XhoI restriction sites using isothermal assembly. The isothermal assembly products were pooled and transformed at high efficiency into electro-competent Escherichia coli. After overnight selection on LB-agar plates with carbenicillin, the bacterial lawn was scraped to extract the pooled plasmids. The pooled plasmids were integrated into either scHP15 WT strain or one of the four deletion strain backgrounds at the ChrI intergenic locus using the same transformation protocol as above. A minimum of 500 S. cerevisiae colonies were pooled for each strain and stored as glycerol stocks. The pooled reporter plasmids and their respective barcodes are provided as part of the data analysis code at https://github.com/rasilab/ribosome_collisions_yeast#high-throughput-sequencing.

Flow cytometry

In total, 3–8 single S. cerevisiae colonies from transformation were inoculated into separate wells of 96-well plates containing 150 μl of yeast extract peptone with 2% dextrose (YEPD) medium in each well and grown overnight for 16 h at 30°C with shaking at 800 rpm. The saturated cultures were diluted 100-fold into 150 μl of fresh YEPD medium. After growing for 6 h under the same conditions as overnight, the plates were placed on ice and analyzed using the 96-well attachment of a flow cytometer (BD FACS Aria or Symphony). Forward scatter (FSC), side scatter (SSC), YFP fluorescence (FITC), and mKate2 fluorescence (PE.Texas.Red) were measured for 10,000 cells in each well. The resulting data in individual .fcs files for each well were combined into a single tab-delimited text file. Mean YFP and mKate2 fluorescence in each well was calculated, and the YFP background from the scHP15 strain and the mKate2 background from the BY4741 strain were subtracted. The ratio of the background-subtracted YFP to mKate2 signal was normalized by that of a standard internal control strain included in all experiments. The average and standard error of this ratio were calculated over all transformant replicate wells and is shown in Figs 1, 5 and 6. Analysis code for reproducing the figures from raw flow cytometry data are available at https://github.com/rasilab/ribosome_collisions_yeast. In total, 2–5 outliers among several hundred wells that were not fluorescent were discarded in some experiments, and they are indicated as such in the analysis code for that experiment. The mKate2 channel measurement was not steady during the flow cytometry measurement in Fig 5D, and hence the YFP fluorescence for this experiment is plotted without normalizing by mKate2 signal. In general, normalizing by mKate2 slightly reduced the fluctuations between biological replicates but did not have any qualitative effect in any experiment. Full analysis code for flow cytometry data is provided in https://github.com/rasilab/ribosome_collisions_yeast#flow-cytometry. The P values in Fig 5D are calculated using the t.test function in R. The analysis code for calculating these P values is provided in https://github.com/rasilab/ribosome_collisions_yeast/blob/master/analysis/flow/hel2_asc1_mutants.md.

Deep-sequencing quantification of mRNA levels

S. cerevisiae glycerol stocks containing pooled reporter strains (> 10⁷ cells) were thawed and grown overnight in 3 ml of YEPD at 30°C in roller drums. The saturated cultures were diluted 200-fold into 50 ml of fresh YEPD medium and grown for 4.5 h at 30°C with shaking at 200 rpm. Each culture was split into two 25-ml aliquots (one for extracting RNA and the other for extracting genomic DNA) into prechilled 50-ml conical tubes and spun down at 3000g, 4°C, 2 min. The supernatant was removed, and the cell pellets were flash frozen in a dry ice/ethanol bath and stored at −80°C.

For RNA extraction, the cell pellets were resuspended in 400 μl of Trizol (15596-026, Thermo Fisher) in a 1.5-ml tube and vortexed with 500 μl of glass beads (G8772, Sigma) at 4°C for 15 min. RNA was extracted from the resulting lysate using the Direct-zol RNA Miniprep Kit (R2070, Zymo) following manufacturer's instructions. For genomic DNA extraction, the cell pellets were processed using the YeaStar Genomic DNA Kit (D2002, Zymo) following manufacturer's instructions.

After RNA extraction, 1 μg of total RNA was reverse transcribed into cDNA using SuperScript IV (18090-010, Thermo Fisher) reverse transcriptase and a primer annealing to the 3′ UTR (oMF321: ACACTCTTTCCCTACACGACGCTCTTCCGATCTGCGTGACATAACTAATTACAT) following manufacturer's instructions. A 252-nucleotide region surrounding the 8-nucleotide barcode was PCR-amplified from both the cDNA and genomic DNA in two rounds. Round 1 PCR was carried out for 10 cycles with oMF321 and oHP139 (GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTAAGGACCCAAACGAAAAG) primers using Phusion polymerase (F530, Thermo Fisher) following manufacturer's instructions. Round 2 PCR was carried out for 10 cycles for genomic DNA samples and 10 or 14 cycles for cDNA samples with a common reverse primer (oAS111: AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTC) and indexed forward primers for pooled high-throughput sequencing of different samples (CAAGCAGAAGACGGCATACGAGATxxxxxxGTGACTGGAGTTCAGACGTGTGCTC, xxxxxx refers to the 6-nucleotide sample barcode and is provided in a table as part of our analysis script). The pooled libraries were sequenced in a HiSeq 2000 (Illumina) in a 50-nucleotide single-end rapid-run mode.

The high-throughput sequencing data were obtained as demultiplexed .fastq files. The barcode corresponding to each reporter was identified in each 50-nucleotide read, and the barcode counts for each reporter were tallied within each sample. The log2 barcode counts for the cDNA were normalized by those of the genomic DNA, after applying a minimum threshold of 100 counts per barcode within each sample. The average and standard error of this log2 ratio was calculated over all barcodes that crossed the 100-count threshold. The average log2 ratio was median-normalized within each set of 5′ UTR variants for a given coding sequence and sample and is shown in Figs 2 and 5. Analysis code for reproducing the figures from raw high-throughput sequencing data is available at https://github.com/rasilab/ribosome_collisions_yeast#high-throughput-sequencing.

Western blotting

Overnight cultures were grown from single S. cerevisiae colonies in 3 ml of YEPD at 30°C in roller drums. The saturated cultures were diluted 100-fold into 3 ml of fresh YEPD medium and grown for 5 h at 30°C in roller drums. The cultures were quickly chilled in an ice-water bath and spun down in 1.5-ml tubes. The cell pellet was washed in 500 μl of water and incubated in 200 μl of 0.1 M NaOH for 5 min at room temperature. The pellet was resuspended in 50 μl of 1× Laemmli buffer, and western blots were performed using standard molecular biology procedures [76]. The anti-FLAG antibody (F1804, Sigma) and the anti-H3 antibody (ab1791, Abcam) were visualized using IRDye antibodies (926-68072 and 926-32211, Licor) on an Odyssey imager. The raw signal in each lane was quantified using ImageJ using the rectangle-select tool followed by the Analyze-Measure tool. For each lane, the FLAG signal was divided by the H3 signal and normalized to a maximum of 10 within each blot. Uncropped images of blots and their quantification are provided in https://github.com/rasilab/ribosome_collisions_yeast/tree/master/data/western. The normalization of lanes is done in https://github.com/rasilab/ribosome_collisions_yeast/blob/master/analysis/western/western_analysis.md.

Kinetic modeling

All scripts and files mentioned below are available at https://github.com/rasilab/ribosome_collisions_yeast/modeling. We specify our kinetic models using the PySB interface [81] to the BioNetGen modeling language [82]. The resulting Python script tasep.py is parsed using the BioNetGen parser into the tasep.bngl file and converted into the tasep.xml file for input to the agent-based stochastic simulator, NFsim [83].

The molecules and the reactions in our model, along with the parameters and initial conditions used in our simulations are described below. A complete list of parameters in our model is provided in S3 Table. The tasep.bngl file provides a machine-readable summary of the description below using BioNetGen language.

Kinetic reactions.

Below, we describe each type of kinetic reaction in our model of quality control during eukaryotic translation. We use a syntax similar to that of BioNetGen [82] for illustrating kinetic reactions.

Transcription produces capped mRNAs without any ribosomes (Eq 1). All mRNA backbone sites and poly-A tail sites are intact. The transcription rate constant k_{transcription} is set to 0.001 s⁻¹. The value of this parameter is arbitrary and is chosen to maintain a pool of 1–2 translatable mRNAs at any given time in the simulation and to produce around 1,000 mRNAs during a typical simulation run for 10⁶ s.

(1)

Translation initiation occurs on capped mRNAs when the start region is free, the start codon is not bonded, and the mRNA backbone at the start codon is intact (Eq 2). It results in the start region being blocked and the ribosome A site being bonded to the start codon (“!1” in Eq 2). The translation initiation rate k_init is systematically varied in our simulations. We use a default value of 0.1 s⁻¹ to match previous estimates [44, 84]. The absolute value of this parameter is not critical in our modeling. Rather, the dimensionless ratio of the translation initiation rate to elongation rate is the important parameter in our model because it determines the frequency of ribosome collisions at the stall.

(2)

Elongation of a ribosome results in the ribosome A site moving from codon c_i to codon c_i+1 (Eq 3). Elongation requires that there is no leading ribosome within a footprint distance at c_i+10. Additionally, when a ribosome elongates at the ninth codon from the initiation codon, the mRNA start region is cleared for a new round of initiation (Eq 4). Finally, elongation of ribosomes that have been hit from the back results in loss of the molecular interaction with the trailing ribosome (“!3” in Eq 5). The elongation rate in all these cases, k_elong, is set to 10 s⁻¹ to match experimental estimates [85, 94].

(3)(4)(5)

Termination results in dissociation of the ribosome–mRNA complex along with production of a protein molecule (Eq 6). Termination of ribosomes that have been hit from the back results in loss of the molecular interaction with the trailing ribosome (Eq 7). The termination rate, k_term, is set to 1 s⁻¹ to be lower than the normal elongation rate of 10 s⁻¹. This choice reflects the observation that ribosome density at stop codons tends to be several-fold higher than at sense codons [40]. The absolute value of k_term is not important in our modeling as long as it is greater than the total rate of initiation and elongation across the mRNA.

(6)(7)

Collision between two ribosomes requires them to be separated by exactly a footprint distance on an intact mRNA, and it results in a bond (“!3” in Eq 8) between the f site of the trailing ribosome and the b site of the leading ribosome. The collision rate constant is set equal to the elongation rate at the A site codon of the trailing ribosome. This choice reflects our assumption that collision occurs when the trailing ribosome tries to translocate to the next codon on the mRNA but is sterically blocked by the leading ribosome.

(8)

Deadenylation of the poly-A tail occurs in the 3′–5′ direction. Deadenylation at position j requires that the A at position j+1 has been removed (Eq 9). This constraint does not apply to the 3′ end of the poly-A tail (at position Lp), whose removal begins the process of deadenylation (Eq 10). We set the deadenylation rate per nucleotide k_{deadenylation} to be 0.03 s⁻¹ and the poly-A tail to be 60 nucleotides to match previous estimates for the PGK1 mRNA [47]. The deadenylation rate sets the overall rate of canonical mRNA decay in our modeling. Our model of deadenylation can be refined further based on recent biochemical and genome-wide studies of poly-A tail metabolism [14, 86]. We do not pursue this here because we do not monitor poly-A tail length, and our focus is on modeling cotranslational quality control.

(9)(10)

Decapping during canonical mRNA decay occurs after the poly-A tail has been fully deadenylated (Eq 11). We set the decapping rate constant k_decapping to be 0.01 s⁻¹, similar in magnitude to previous work [47]. Because we do not know the initiation rate of deadenylated but capped mRNAs, we assume that mRNAs are initiated at their normal efficiency right until the cap is removed (Eq 2).

(11)

5′–3′ exonucleolysis at position i requires that the mRNA backbone at position i−1 has been removed (Eq 12) and that the position is not blocked by a ribosome. 5′–3′ exonucleolysis can initiate during canonical mRNA decay only after the mRNA has been decapped (Eq 13). We set the 5′–3′ exonucleolysis rate constant k_{exo_53} to be 1 s⁻¹ [47], such that it is faster than the total rate of deadenylation and decapping.

(12)(13)

3′–5′ exonucleolysis at position i requires that the mRNA backbone at position i+1 has been removed (Eq 14) and that the position i is not blocked by a ribosome. The 3′–5′ exonucleolysis can initiate during canonical mRNA decay only after the mRNA has been fully deadenylated (Eq 15). We set the 3′–5′ exonucleolysis rate constant k_{exo_35} to be 0. This choice is consistent with the much slower observed rate of 3′–5′ exonucleolysis in comparison with 5′–3′ exonucleolysis in S. cerevisiae [47]. Furthermore, any mRNA undergoing 3′–5′ exonucleolysis will not contribute to the pool of full-length proteins or mRNA molecules, which are the experimental observables of interest in this work.

(14)(15)

Abortive (premature) termination results in dissociation of the aborting ribosome from the mRNA. Abortive termination of a ribosome also results in the loss of molecular interaction with any collided ribosome in the front or the back. In our modeling of collision-stimulated quality control from intact mRNAs, ribosomes have different rates of abortive termination depending on whether they have undergone collision from the front (k_{abort_front_hit}, Eq 18), back (k_{abort_back_hit}, Eq 17), or both (k_{abort_both_hit}, Eq 19) sides. Ribosomes that have not undergone collisions also abortively terminate (k_{abort_no_hit}, Eq 16). The different models of abortive termination in our work (Fig 3A) correspond to different combinations of the four types of abortive termination. In the SAT model (Eq 8 in Fig 3A), k_{abort_no_hit} = k_{abort_back_hit}>0 and k_{abort_front_hit} = k_{abort_both_hit} = 0. This choice implies that collisions do not have any stimulatory effect on abortive termination and that ribosomes that are stacked behind the leading stalled ribosome do not undergo abortive termination. In the CSAT model (Fig 3A), k_{abort_back_hit} = k_{abort_both_hit}>0 and k_{abort_no_hit} = k_{abort_front_hit} = 0. This choice implies that collisions from the back are necessary for abortive termination, but ribosomes that are stalled either on their own or because of a leading stacked ribosome do not undergo abortive termination. In the CAT model (Fig 3A), k_{abort_front_hit} = k_{abort_both_hit}>0 and k_{abort_no_hit} = k_{abort_back_hit} = 0. This choice implies that collisions from the front are necessary for abortive termination, but ribosomes that are stalled without a leading stacked ribosome do not undergo abortive termination. The in vivo abortive termination rates are not known. In each of the above three models, we chose the respective nonzero abortive termination rates to be such that they would decrease protein expression appreciably because of stalling.

(16)(17)(18)(19)

We also modeled abortive termination of ribosomes with a truncated mRNA in its A site using the same rules as above but with r_i~endocleaved instead of r_i~intact. This is necessary for releasing all the ribosomes that are 5′ to the endonucleolytic cleavage site. We assume that ribosomes abort with a uniformly high rate of 1 s⁻¹ from mRNAs with cleaved A sites independent of whether they have undergone collision. This is consistent with previous in vitro experiments showing that Dom34/Hbs1 can efficiently recycle ribosomes that are stalled on a truncated mRNA [19, 20].

Endonucleolytic mRNA cleavage on translating ribosomes results in the mRNA being cleaved at a distance Lc 5′ to the A site. We further assume that the endonucleolytic cleavage triggers immediate decapping. This assumption is consistent with increased rates of degradation of mRNAs that have undergone endonucleolytic cleavage [59]. This assumption does not affect any of the conclusions in this work because further rounds of initiation of cleaved mRNAs will not result in the production of full-length proteins, the observable of interest to us. In our modeling of CSEC, we consider different rates of endonucleolytic mRNA cleavage depending on whether ribosomes have undergone collision from the front (k_{cleave_front_hit}, Eq 21), back (k_{cleave_back_hit}, Eq 22), or both (k_{cleave_both_hit}, Eq 23) sides. Ribosomes that have not undergone collisions also trigger endonucleolytic mRNA cleavage (k_{cleave_no_hit}, Eq 20). The different models of cotranslational endonucleolytic cleavage in our work (Fig 4A) correspond to different combinations of the four types of endonucleolytic cleavage. In the SEC model (Fig 4A), k_{cleave_no_hit} = k_{cleave_back_hit}>0 and k_{cleave_front_hit} = k_{cleave_both_hit} = 0. This choice implies that ribosome collisions do not have any stimulatory effect on endonucleolytic mRNA cleavage and that ribosomes that are stacked behind the leading stalled ribosome do not cause endonucleolytic mRNA cleavage. In the CSEC model (Fig 4A), k_{cleave_back_hit} = k_{cleave_both_hit}>0 and k_{cleave_no_hit} = k_{cleave_front_hit} = 0. This choice implies that collisions from the back are necessary for stimulating endonucleolytic mRNA cleavage, but ribosomes that are stalled either on their own or because of a leading stacked ribosome do not result in endonucleolytic cleavage. Unlike the case of abortive termination, we do not separately model the effect of the trailing versus leading ribosome in a collision stimulating endonucleolytic cleavage, because both these models will result in a cleaved mRNA.

(20)(21)(22)(23)

We set Lc = 10 codons in our simulations, even though experiments indicate that endonucleolytic mRNA cleavage occurs 14–15 nucleotides behind the stalled ribosome [52, 53]. This choice of Lc to be equal to the ribosome footprint size simplifies the specification of reaction rules but does not affect the predictions for either protein output or mRNA lifetime. This can be set to the more physiological value of 5 in case the model predictions for cleaved fragments are examined (such as when comparing against short ribosome footprints [52]). The in vivo endonucleolytic cleavage rates are not known. In the above two models, we chose the respective nonzero endonucleolytic mRNA cleavage rates to be such that they would decrease mRNA lifetimes appreciably because of stalling.

Analysis of FACS-Seq data from Dvir and colleagues

S1 Table was downloaded from http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1222534110/-/DCSupplemental/sd01.xlsx. The exp column in this table was averaged over all sequence_ variant rows that share the same last 4 nucleotides. The mean and standard error of this average are shown in S1 Fig panel C. This analysis is described at https://github.com/rasilab/ribosome_collisions_yeast/#analysis-of-dvir-et-al-2013-data.

Ribosome profiling analysis

Raw .fastq files from a previous study [62] were downloaded for the SRA experiments, SRX391789 (RNA-seq) and SRX391790 (Ribo-seq). The adapters were removed using cutadapt [88] with the arguments --cut=8 --adapter=TCGTATGCCGTCTTCTGCTTG --minimum-length=15. The trimmed reads were first depleted of rRNA contaminant reads and then aligned to the S. cerevisiae transcriptome and genome (sacCer3) using tophat 2.0.14 [89] with the arguments --bowtie1 --num-threads 8 --no-novel-juncs --library-type fr-unstranded --keep-tmp --read-mismatches 2. The resulting .bam files were sorted and indexed using samtools [90]. The aligned reads were assigned to the 13^th nucleotide from the 5′ end, and this was used to calculate coverage across the sacCer3 genome.

To identify stall sequences, the protein-coding sequences not marked as “Dubious” in the sacCer3 genome (annotation set: saccharomyces_cerevisiae_R64-1-1) were scanned for 10-codon windows that contained at least 6 lysine and arginine codons or 6 proline codons. The control sequences were similarly identified by scanning for windows that contained at least 6 glutamate and aspartate codons. This analysis is described at https://github.com/rasilab/ribosome_collisions_yeast/#identification-of-putative-rqc-stalls.

The Ribo-seq coverage in the 150-nucleotide window for each stall or control sequence was normalized by the mean coverage within that window. This normalized coverage was then averaged at each position in the window across all stall or control sequences and shown in Fig 6A, S5 Fig panel A.

To calculate TE of S. cerevisiae genes, the Ribo-seq and RNA-seq rpkm values for each gene were downloaded from ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE75nnn/GSE75897/suppl/. TE for each gene was calculated as the ratio of the Ribo-seq rpkm to RNA-seq rpkm. Only genes with a minimum of 5 rpkm in each sample were used for the TE calculation.

To calculate TE of mRNA regions 5′ to stalls, we restricted the Ribo-seq and RNA-seq coverage to this region for stall-containing mRNAs and calculated TE for mRNAs that contained at least 100 reads in this region. For control mRNAs, we restricted the coverage to regions 5′ to the median-stall location on stall-containing mRNAs (215 codons). We used this restricted region of comparison to account for the 5′– 3′ skew in ribosome density along all S. cerevisiae mRNAs [62]. This analysis is described at https://github.com/rasilab/ribosome_collisions_yeast/#analysis-of-weinberg-et-al-2016-data.

RNA-seq analysis

Raw .fastq files from a previous study [57] were used for differential gene expression analysis between ΔHEL2, ΔASC1, and WT strains. Only RNA-seq data from strains containing control reporters were used for the analysis. Any ambiguous terminal nucleotides were trimmed using cutadapt [88] with the argument --trim-n. The trimmed reads were first depleted of rRNA contaminant reads and then aligned to the S. cerevisiae genome (sacCer3) using bowtie 1.1.1 [91] with the arguments -v 2 --un /dev/null --threads 8 --sam. The resulting .sam files were converted to the binary .bam format, sorted, and indexed using samtools [90]. The alignments were used to calculate the number of reads aligning to each transcript in the TxDb.Scerevisiae.UCSC.sacCer3.sgdGene R package using the findOverlaps function from GenomicAlignments R package [92]. Transcripts that had a minimum of 100 alignments in each sample were used for differential gene expression analysis using DESeq2 with default parameters [93]. The ΔHEL2 and ΔASC1 samples were treated as replicates for DESeq2 input and compared with the WT strain. The log2FoldChange output from the results function in DESeq2 was used for plotting Fig 6D. This analysis is described at https://github.com/rasilab/ribosome_collisions_yeast/#analysis-of-sitron-et-al-2017-data.