Coding regions affect mRNA stability in human cells

1Transcript destruction is key for controlling gene expression. Messenger RNA (mRNA) decay not only helps set the overall expression level of genes, but it also determines the dynamics of their expression. Unsurprisingly, proper control of mRNA decay is critical in diverse biological processes ranging from the maternal-to-zygotic transition to the inflammatory response. A major line of research has thus been to understand both the elements within an RNA that control its stability as well as the mechanisms underlying this regulation.

2For decades, researchers have focused on the 3′ untranslated regions (3′UTRs) for several important reasons. First, especially in humans, 3′UTRs contain extensive sequence space that has large regulatory potential while also being devoid of other constraints (such as encoding a protein). Second, regulatory factors, such as RNA binding proteins (RBP), can stably bind the 3′UTR, which is not exposed to the translating ribosome, unlike the 5′UTR or open reading frame (ORF). Thus, elements within the 3′UTR, such as those bound by microRNAs (miRNAs), provide our best understanding of the mechanisms and regulation of mRNA decay. Briefly, sites in the 3′UTR are recognized by regulatory factors, such as RBP or miRNAs. These regulatory factors both remove stabilizing protein components and also recruit deadenylase complexes, which then shorten the poly(A) tail found on nearly all eukaryotic transcripts. Following deadenylation, the decapping enzyme then removes the 5′ cap, thereby exposing the body of the transcript to the major cytoplasmic 5′ → 3′ exonuclease, Xrn1.

3In contrast, outside of surveillance pathways, our understanding of how the coding region influences mRNA stability lags behind that of 3′UTR-based regulation. Nonetheless, over the past few years, research in a number of model organisms has demonstrated that coding regions play a major role in mRNA stability. Best defined by work in Saccharomyces cerevisiae, the primary mechanism by which ORFs regulate mRNA stability is through codon usage. Although synonymous codons encode the same amino acid, they are not equivalent to the elongating ribosome. Due to differences in tRNA abundance, demand, and wobble interactions, some synonymous codons are decoded faster than others. These differences in elongation speed are sensed by the decay machinery, resulting in recruitment of the CCR4-NOT deadenylase complex, the decapping activator Dhh1p, and eventual decapping. In other words, mRNAs with low optimality codons (i.e., those decoded slowly) are more unstable than those with high optimality codons. Although initially made in budding yeast, these observations have now been extended to Escherichia coli, Schizosaccharomyces pombe, trypanosomes, and metazoans, likely pointing to a conserved evolutionary mechanism.

4An important unanswered question has been the extent to which these observations hold in humans. This question has been raised in part because human 3′UTRs are often longer than 1000 nt, which has led to the model that 3′UTRs are the dominant regulatory force for post-transcriptional regulation in humans. Importantly, all elements of a gene, including the coding region and UTRs, coevolve, and so exploring the contribution of codons to endogenous mRNA stability is confounded by UTR effects. Another barrier has been a lack in the consensus of which codons mediate fast elongation speeds in human cells, and so calling codons as “optimal” or “nonoptimal” has been challenging. The lack of consensus has come about for at least three reasons. First, unlike in yeast, only slight relationships have been observed between translation elongation speeds and standard metrics of codon optimality, such as tRNA abundance and codon adaptivity index. Second, tRNA abundance is notoriously hard to measure due to extensive RNA modifications and secondary structure. Highly processive reverse transcriptases, like TGIRT, and other high-throughput sequencing methods have enabled better quantitation, but measurements from different labs still fail to correlate well, even when performed on the same cell type. Third, tRNA abundance also differs between different cell types and states, leaving open the possibility that some variation in tRNA measurements and elongation speeds reflect biological differences.

5To understand how the coding region could impact mRNA stability in humans, we turned to the human ORFeome collection, which contains full-length human ORFs surrounded by invariant UTRs. By using this collection as a parallel reporter system, we isolated the effects of the coding region on mRNA stability. We found that different ORFs mediate a wide range of mRNA stabilities and that the majority of their impact depends on translation. Similar to results from model organisms and a parallel study in humans, our results point to codon usage as a key feature of the ORF that alters mRNA stability. Moreover, our analyses revealed that codons associated with instability have longer A-site dwell times, thus providing a link between elongation speeds and mRNA decay in humans.

6Given recent reports on the influence of codon use on mRNA stability, we wanted to explore the potential of the ORF to impact mRNA stability in human cells. However, codon use and other ORF features coevolve with mRNA characteristics (such as length, translational efficiency, 3′UTR regulation), all of which have the potential to impact mRNA stability. We reasoned that these confounding factors would significantly impact genome-wide analyses, making it challenging to disentangle correlative and causative features. To explore this question, we thus turned to the human ORFeome collection because it contains thousands of coding regions surrounded by invariant UTRs as lentiviral constructs.

7We decided to use the ORFeome collection as a parallel reporter system and split the collection into six pools, which we then used to generate six corresponding sets of pooled stable HEK293T cell lines through lentiviral infection. Western blotting against the V5 tag at the carboxy terminus of each ORFeome construct confirmed that many ORFs were expressed after stable line generation (Supplemental Fig. S1A). We then used approach-to-equilibrium 4-thiouracil (4SU) metabolic labeling followed by biotinylation and streptavidin pull-down to measure mRNA stability transcriptome-wide.

8We did not use ORFeome-targeted RNA sequencing, and so we simultaneously captured stabilities for both ORFeome-derived and endogenous mRNAs. However, because very few reads mapped to regions that could be unambiguously defined as ORFeome-derived (i.e., mapping to the junctions between the UTRs and coding region), these reads could not be used to generate half-lives. To enable half-life calculation of the ORFeome-transcripts, we developed a computational pipeline to classify each mRNA (and corresponding half-life) as “ORFeome” or “endogenous” (Fig. 1B). Briefly, reads were mapped to the RefSeq coding regions. To classify an mRNA as “ORFeome,” we required three times more reads in the corresponding cell line in our steady-state libraries than in the paired cell line. We first analyzed the experiment derived from cell lines expressing pools 1 and 4. Of the 2211 pool 1 ORFs that were detected in cell line 1, only 359 were expressed three-times as much in cell line 1 than in cell line 4 (presumably because some were expressed endogenously, while others were not expressed at high enough levels). Similarly, 302 pool 4 ORFs passed this threshold (Fig. 1B). As expected based on this criterion, pool 1 mRNAs were significantly more highly expressed in pool 1 cells, and pool 4 mRNAs, in the pool 4 cells (Supplemental Fig. S1B,C). Coding regions contained in the ORFeome pool 1 that did not cross this threshold were excluded from downstream analysis. We then calculated half-lives for three sets of mRNAs (pool 1, pool 4, and endogenous). Some ORFeome mRNAs, despite passing our initial cut-off, were still lowly expressed, and thus we could not calculate half-lives for all of the ORFeome mRNAs. In total, we obtained half-lives for 357 ORFeome mRNAs (221 from pool 1 and 136 from pool 4), and 10,208 endogenous mRNAs (Supplemental Tables S1, S2).

9Having thus obtained half-lives for 357 ORFeome reporters, we next asked whether changing the coding region was sufficient to change mRNA stability. In general, ORFeome mRNAs were more stable than endogenous mRNAs (median half-life: 8.2 h for endogenous mRNAs versus 13.8 h for ORFeome mRNAs; P < 10−15, Kolmogorov–Smirnov test; Fig. 1C), likely due to the WPRE stabilizing element included in the invariant 3′UTR of the ORFeome collection. Strikingly, despite containing invariant UTRs, the ORFeome mRNAs showed similar variation in stability as endogenous transcripts: half-lives ranged from 0.25 h to >100 h for endogenous mRNAs, and from 0.56 h to >100 h for ORFeome mRNAs (Fig. 1C,D). This result held even when we subsampled endogenous mRNAs to match sample sizes between the two groups (data not shown).

10To see whether this result held more broadly, we turned to our other cell lines, obtaining half-lives for (i) 357 ORFeome mRNAs from pools 2 and 3 and (ii) 384 ORFeome mRNAs from pools 5 and 6 (Supplemental Tables S1, S3, S4). Similar to what we observed with pools 1 and 4, the ORFeome mRNAs had similar variation in stabilities as endogenous ones (Supplemental Fig. S1E). Finally, we analyzed a complementary ORFeome data set from ). This data set allowed us to directly compare the variability in decay rates for the same ORFs expressed as endogenous or ORFeome-derived transcripts (Supplemental Fig. S1F). Importantly, there was a similar range in stabilities for these matched ORFeome and endogenous transcripts. Thus, we conclude that the coding region can regulate mRNA stability in human cells. One surprising implication of this result is that merely changing the coding region can give a similar range of stabilities seen with endogenous mRNAs, which differ not only in their coding regions, but also in their UTRs.

11The primary function of the coding region of a transcript is to be read by the ribosome, and it has become increasingly clear that translation impacts mRNA stability not just through surveillance pathways, but also through normal decay. Given the connection between translation, the coding region, and mRNA stability, we next asked whether the effects of the coding region on mRNA stability required ribosome loading. To do so, we repeated our ORFeome and endogenous stability measurements (in cell lines 1 and 4), but this time cells were treated either with DMSO or the translation inhibitor 4EGI-1, which disrupts eIF4F and so blocks initiation.

12We first confirmed that 4EGI-1 reduced translation by analyzing polysomes by sucrose gradient fractionation and measuring puromycin incorporation followed by western blotting (Fig. 2A,B). As expected, 4EGI-1 treatment robustly reduced translation. (Note that complete inhibition of translation was incompatible with mRNA stability measurements.) 4EGI-1 treatment had broad effects on mRNA stability, both for endogenous and ORFeome mRNAs. Endogenous mRNAs were globally less stable in the presence of the inhibitor (Fig. 2C), consistent with previous studies where translational inhibition destabilized mRNA. While there was strong correlation for endogenous mRNAs measured in our original experiment and in the presence of DMSO (Spearman r [rs] = 0.67, P < 10−15, Supplemental Fig. S2), there was only a modest correlation with those measured in the presence of 4EGI-1 (rs = 0.34, P < 10−15, Supplemental Fig. S2). Surprisingly, the variation in mRNA decay rates was significantly, albeit modestly, reduced in the presence of 4EGI-1 (σ2 4EGI-1/σ2 DMSO = 0.65, P < 10−15, Fig. 2D).

13We saw similar results with the ORFeome mRNAs, although the effects were more pronounced here. Like the endogenous mRNAs, ORFeome mRNAs were also less stable in the presence of 4EGI-1 (P < 10−10, Kolmogorov–Smirnov test; Fig. 2E). More strikingly, there was substantially less variation in the stabilities of ORFeome mRNAs in the presence of 4EGI-1 (σ2 4EGI-1/σ2 DMSO = 0.40, P < 10−15, Fig. 2F). These results thus suggest that coding region predominantly requires translation to impact stability. It may be that the reduced variability in the stabilities of endogenous mRNAs reflects the contribution of coding regions to endogenous transcript stability. In addition, because ribosome loading was incompletely blocked and inhibiting translation also exposes previously inaccessible regulatory sites in the coding region (see below), the differences observed with 4EGI-1 are likely a lower bound for the influence of translation of the coding region on mRNA decay.

14We next wanted to identify elements within the coding region that influenced RNA stability. Many different ORF features have been associated with mRNA stability: transcript or ORF length; local secondary structure; 3′UTR-like regulation; and codon usage. For those ORF elements directly regulating mRNA stability, we expected that their signal would be strengthened in the ORFeome data set (compared with the endogenous one) because confounding UTR effects are removed in the ORFeome measurements. In contrast, for elements that coevolve with true causative features, but do not directly affect stability, we expected that their signals would be weaker in the ORFeome data set than in the endogenous one.

15We first investigated ORF length in our data from cell lines 1 and 4. Consistent with many previous observations, for endogenous mRNAs, there was a negative correlation between ORF length and mRNA stability such that transcripts with longer coding regions were more unstable (rs = −0.14, P < 10−15, Supplemental Fig. S3A). In contrast, we observed no correlation between ORF length and mRNA stability for ORFeome transcripts (rs = 0.01, P = 0.9, Supplemental Fig. S3). We observed similar results in our other ORFeome pools and the complementary published ORFeome data, where the correlation between length and stability was stronger for endogenous transcripts than for ORFeome-derived ones; (endogenous rs: −0.30 to −0.15 versus ORFeome rs: −0.22 to 0.01; Supplemental Fig. S3B). Thus, we conclude that ORF length, at most, only weakly affects mRNA stability.

16We next looked at the local secondary structure in the ORF. To do so, we calculated the folding energy within a 100-base sliding window and took the lowest value for each coding region. For both endogenous and ORFeome mRNAs, we found no general relationship between transcript stability and secondary structure (endogenous rs: −0.10 to 0.04 versus ORFeome rs: −0.06 to 0.07, Supplemental Fig. S3C,D). Although it is possible that strong secondary structures influence stabilities for only some transcripts, it appears unable to provide a general explanation for the range of stabilities seen with the ORFeome mRNAs.

17We next examined the potential of coding regions to modulate mRNA stability through 3′UTR-like regulation. We first divided ORFeome mRNAs into those containing or lacking sites for the top five expressed microRNA (miRNA) families. Consistent with individual miRNA sites in the ORF having little impact, half-lives for these two sets of ORFeome transcripts failed to be significant for any of the ORFeome data sets (P = 0.06 to 0.99, Supplemental Fig. S3E,F). Finally, we compared stabilities for ORFeome mRNAs containing and lacking AU-rich elements, and only observed modest differences, although these were significant for two of the four ORFeome data sets (P = 10−6 to 0.9, Supplemental Fig. S3G,H). Thus, although many coding regions contain sites that have the potential to be recognized by 3′UTR regulatory factors, these sites cannot explain the range of ORFeome half-lives.

18This conclusion is consistent with previous reports showing that miRNA sites in the ORF are rarely effective, because of the translating ribosome. This result also suggests that some of the residual variation in ORFeome mRNA stability upon 4EGI-1 treatment (Fig. 2F) may be due to RBP-based regulation that now has an opportunity to impact stability.

19We next turned to the role of codon use in mRNA stability. To do so, we used a previously defined metric, the codon stability coefficient or CSC. For each nonstop codon, we calculated its frequency in each ORF and measured the Spearman correlation between the codon frequency with the associated mRNA half-life across the collection of endogenous or ORFeome mRNAs (giving rise to endogenous and ORFeome CSC values, respectively). We calculated these scores from each of our three ORFeome/endogenous stability data sets, as well as those generated during DMSO treatment (Fig. 2), giving us four independent measurements for each CSC value.

20We calculated the average score of the four values for each codon (Fig. 3A,B; Supplemental Table S5). Overall, the average endogenous and ORFeome CSC values correlated with each other (rs = 0.70, P < 10−15; Fig. 3C). Consistent with these effects being mediated by translation, ORFeome values calculated based on computationally +1 or +2 frameshifted codons were significantly less correlated with endogenous CSC values (rs = 0.32 and 0.22, respectively; P < 0.0002, Fishers r-to-z transformation). Similarly, there was not a significant correlation between the ORFeome values and those from yeast, an unrelated species (rs = 0.22, P = 0.07).

21In investigating the effects of codons on mRNA stability, much of the focus has been on the effects of synonymous codons because of the model that the abundance of cognate tRNAs differs, leading to differences in decoding of synonymous codons and thus mRNA stability. However, we noted that, surprisingly, some synonymous codons tended to have similar relationships with stability, especially for the “unstable” ORFeome CSCs (Supplemental Fig. S4A). To investigate this issue more thoroughly, as has been done before, we calculated an “amino acid stabilization coefficient” (AASC) for the 20 standard amino acids by averaging the CSCs for the corresponding codons (Supplemental Fig. S4B; Supplemental Table S6). As with CSCs, average endogenous and ORFeome AASCs were correlated (rs = 0.60, P = 0.007). Although codons for several amino acids, such as leucine, showed a wide range in stability scores, codons for other amino acids (such as lysine and glutamate) had similar values. Taken together, these data demonstrate that codon usage and mRNA stability are related in human cells, as has been seen in other organisms, and that, for some amino acids, synonymous codons appear to behave similarly.

22Based on studies in model organisms, the current model explaining the relationship between codon usage and mRNA stability is that some codons are translated slowly, and slow elongation in turn triggers mRNA degradation. We thus asked whether this model held in human cells by investigating the relationship between stability coefficients and elongation rates. To do so, we used a published ribosome profiling data set from HeLa cells and calculated the enrichment of each codon (when in the inferred A site) in ribosome-protected fragments compared with its background abundance in the pool of translated transcripts (see Materials and Methods, Supplemental Table S5). That is, high enrichment indicates long dwell times and thus slow relative elongation; low enrichment, short dwell times and fast relative elongation. As expected, we observed a strong enrichment for stop codons in the A site, and depletion in the P site (data not shown). We then compared these relative dwell times with the endogenous and ORFeome CSCs. When we examined the codon scores from cell lines 1 and 4, there was no relationship between endogenous CSCs and dwell times (rs = 0.02, P = 0.83), but there was a significant correlation with the ORFeome values (rs = 0.33, P = 0.008) such that codons translated more slowly were associated with instability (Fig. 4A,B).

23To explore the relationship between codon stability measurements and elongation speeds more fully, we calculated codon-level dwell times from four other ribosome profiling data sets: a biological replicate from HeLa cells; HEK293T cells; U2OS cells; and NIH3T3 cells. We next compared our four sets of ORFeome CSCs with the calculated dwell times from both HeLa replicates, as well as HEK293T, U2OS, and NIH3T3 experiments, given 20 total comparisons (Fig. 4C). Consistent with the known differences in codon usage and tRNA abundance between mouse and human cells, the ORFeome CSCs were significantly better correlated with human ribosome speeds than NIH3T3 ribosome speeds (human line rs values: 0.16 to 0.34; NIH3T3 rs values: 0.01 to 0.10; P = 0.003, Wilcoxon test; Fig. 4C). Importantly, when we repeated this analysis with the corresponding endogenous CSC values, the ORFeome values were significantly better correlated with human ribosome dwell times than the endogenous values (P = 0.0005, paired Wilcoxon test; Fig. 4D).

24We next compared the ORFeome CSC values determined in the contemporaneous study with A-site dwell times. As with our ORFeome values, these ORFeome CSCs were significantly more correlated with human ribosome dwell times than those in NIH3T3 cells (P = 0.04, Wilcoxon test; Fig. 4E). Similarly, their ORFeome CSCs performed significantly better than the corresponding endogenous values (P = 0.0001, Wilcoxon test; Fig. 4F). Interestingly, despite their CSC values being derived from many more ORFs than our values, there was no difference in the correlations of either set with ribosome elongation speeds (P = 1, Wilcoxon test). Taken together, these results indicate that codons associated with instability—but only when determined by ORFeome, not endogenous, measurements—are translated more slowly.

25Although the A-site codon is thought to predominantly influence elongation speeds, other steps in elongation, such as peptide bond formation or translocation, can also influence ribosome speed and so could conceivably also impact mRNA stability. To explore this possibility, we calculated dwell times as before, but this time for codons in the inferred P and E sites, and compared the four sets of human dwell times for codons in the A, P, and E sites with our four sets of ORFeome CSC measurements.

26Although correlations between CSCs and A-site speeds were tightly clustered, there was a broader range with the P- and E-site speeds (Fig. 4G). Overall, however, the effect size of the correlations was negligible for P- and E-site speeds (median correlation: 0.29, 0.01, and 0.01 for A, P, and E sites, respectively), and A-site speeds correlated significantly better with ORFeome CSC values than P- or E-site speeds (P = 0.003 and 0.01, respectively, paired Wilcoxon test). These results suggest that the A-site speed is likely the major contributor to ORF-mediated regulation of mRNA stability.

27Inspired by reports of coding regions affecting mRNA stability in model organisms, we set out to investigate this phenomenon in human cells. By using the ORFeome collection, where human coding regions are surrounded by invariant UTRs, we showed that ORFs impact mRNA stability in human cells and that the majority of the effect depends on translation. Through our ORFeome stability measurements, we were able to explore features that affect stability, and, in contrast to previous suggestions, we found little evidence that ORF length, local secondary structure, and 3′UTR-like regulation generally impact mRNA stability. Some features (e.g., miRNA-mediated regulation) likely regulate specific mRNAs; but others (e.g., length) may coevolve with direct effectors of RNA stability, explaining results from studies that relied on measurements of endogenous transcript stability.

28Instead, coding regions seem to affect mRNA stability through differences in codon usage. Based on our results, experiments in model organisms, and an independent study in human cells, we propose that the relationship between codon use and mRNA stability is causative, and that instability codons stimulate mRNA decay by slowing translation elongation.

29One of the strongest lines of evidence for the connection between mRNA stability and translation elongation speeds comes from a comparison of codon stability scores from ORFeome and endogenous mRNAs with elongation speeds. While A-site dwell times did not correlate with endogenous scores (which derive from a mix of ORF- and UTR-based regulation), we saw a reproducible, significant correlation with ORFeome stability scores (which derive from only ORF-based regulation).

30Because we found the strongest association between instability-associated codons and amino acids with A-site speeds, we favor a model that the elongation speeds (and thus mRNA stability) are predominantly determined by decoding rates, although there may be some amino acids (such as proline) that slow elongation by other mechanisms (such as peptide bond formation). In turn, decoding rates are probably determined by a combination of tRNA abundance and aminoacylation. Interestingly, we found that synonymous codons tend to have similar A-site dwell times (Supplemental Fig. S4B), and such “amino acid” effects have now been seen in several studies in mammalian cells. Moreover, we found that A-site dwell times for different amino acids differ between HeLa/HEK293T and U2OS cells (Supplemental Fig. S4C,D). Thus, we propose that the combination of tRNA abundance and charging affect decoding rates and mRNA stability in human cells.

31Nonetheless, an unanswered question is why the codon stability scores do not correlate better with A-site dwell times. There are several possible technical and biological explanations. For instance, this result may be a technical artifact of ribosome profiling because it is known that ribosomes can carry out one or two rounds of elongation during the polysome preparation, which could blur the resolution of elongation speeds. Arguing against this possibility is that the P site and A site elongation speeds are very different. Similarly, it is unlikely that this result is due to only using a few hundred transcripts for the ORFeome scores: Corresponding scores from a parallel paper, which were based on thousands of transcripts, did not perform any better.

32Instead, we suspect there may be several biological explanations. First, within endogenous coding regions, there is a co-usage of codons and amino acids, which may partially confound the stability scores for individual codons and amino acids. Second, in yeast, di- and tricodons have been shown to dramatically affect protein output. Such multicodon effects are not captured by either our stability scores or current elongation speed measurements. Finally, although our data indicate that coding regions primarily affect mRNA stability through translation, translation-independent effects remain unaccounted for in our current analysis and could explain some of the stability variation observed in the ORFeome.

33One final issue remains: What is the overall contribution of the ORF to mRNA stability? Although our work and that from others indicate that coding regions impact mRNA stability, the amount to which they do so remains unclear. Based on what we have already found, determining the overall contribution of the coding region to mRNA stability will be a challenge, and the answer probably differs for different transcripts and between different cell types. That is, because the effect of a coding region depends upon ribosome loading, endogenous mRNAs, which have different translation initiation rates, will differ in their opportunity to be regulated by elongation speed. Indeed, it is likely that those mRNAs most affected are those both with high initiation rates and containing destabilizing codons. Another important consideration is that tRNA abundance and charging change between different cell states, and so the identities of mRNAs with destabilizing coding regions may also change between different cell types. One important wrinkle is the potential for interplay between ORF and 3′UTR regulation. For instance, a recent report indicated that, in Drosophila S2 cells, the coding region can influence the amount of regulation exerted by 3′UTR sites, but it is unknown whether the same interaction occurs in human cells. Although integrating the contributions of different parts of an mRNA on stability requires additional work, it is the critical next step for understanding how gene expression is controlled.

34Human HEK293T cells were cultured in DMEM (Lonza) supplemented with 10% FBS (VWR Seradigm) and 1% penicillin–streptomycin solution. Cell lines were cultured at 37°C in a humidified incubator with 5% CO2.

35Drosophila melanogaster Schneider 2 (S2) cells (Cat #R69007, Thermo Fisher Scientific) were cultured in ExpressFive SFM media (Thermo Fisher Scientific) supplemented with 10% heat-inactivated FBS (Wisent) and 20mM l-Glutamine (Thermo Fisher Scientific) at 28°C.

36S. cerevisiae USY006 was grown in YPD liquid or plates at 30°C. Cultures were obtained from Dr. John Rubinstein (The Hospital for Sick Children).

37The human ORFeome collection version 8.1 (ccsbBroad304) cloned into lentiviral vector pLX304 was obtained from Dr. Jason Moffat (University of Toronto) as a series of 96-well overnight bacterial cultures. Equal volumes of bacterial cultures were pooled into 36 pools comprising ∼576 clones each. Plasmid DNA was isolated using the GeneJET Plasmid Midiprep Kit (Thermo Fisher Scientific) as per manufacturer's instructions, yielding an average of ∼70 µg plasmid DNA per pool. The 36 isolated pools were further combined into six unique pools for downstream cell line generation.

38Each of the six unique virus pools were packaged by transfection into HEK293T cells using Lipofectamine 2000 (Thermo Fisher Scientific), according to manufacturer's instruction. Cells were transfected with the lentivirus pLX304 pool, psPAX2 packaging vector, and pVSV-G envelope vector. After 8 h, transfection media was then removed and switched to harvest media (DMEM + 10% FBS + 1.1 g/100 mL BSA [7.5% solution, Thermo Fisher Scientific]). Cells were left for 2 d to complete virus production. Media was then collected from the plate and filtered through a 0.45 µm filter by syringe. Harvested viruses were aliquoted.

39Freshly thawed HEK293T cells were grown in 10 cm dishes to reach ∼30%–50% confluence for the day of infection. Media was removed and 9 mL prewarmed infection media (DMEM + 10% FBS + 8 µg/mL Polybrene) was added to cells. A total of 2 mL of freshly harvested virus pool was added to one plate of HEK293T cells each and incubated overnight. Cells were then trypsinized and expanded into 15 cm dishes. Cells were selected using selection media (DMEM + 10% FBS + 6 µg/mL Blasticidin [BioShop]) for 6 d. Selection media was changed every day. Cells were then frozen in cell freezing medium (Sigma-Aldrich) and stored in liquid nitrogen.

40Cells were harvested by trypsinization and pelleted by centrifugation at 1000g for 2 min at 4°C. Cell pellets were resuspended in 500 µL lysis buffer (100 mM KCl, 0.1 mM EDTA, 20 mM HEPES-KOH pH 7.6, 0.4% NP-40, 10% glycerol, 1 mM DTT, complete mini EDTA-free protease inhibitors [Roche]) and clarified at 21,000g for 5 min at 4°C. A total of 250 µL supernatant was mixed with 20 µL 4× Bolt LDS sample buffer (Thermo Fisher Scientific), 8 µL 10× Bolt sample reducing agent (Invitrogen) and proteins were denatured at 75°C for 10 min. Protein samples were loaded into Bolt 4%–12% Bis-TRIS Plus gels (Thermo Fisher Scientific) and run at 160 V for ∼1 h. The gel was transferred onto an Amersham Hybond PVDF membrane (GE Healthcare), according to manufacturer's instructions. Primary antibodies were added at 1:10,000 concentration for α-V5 antibodies (Cat #V8012, Sigma-Aldrich), 1: 10,000 for α-puromycin antibodies (Cat #3RH11, Cedarlane), and 1: 5000 for α-tubulin antibodies (Cat #T5168 Sigma-Aldrich). α-mouse secondary antibody (Cat #7076, New England Biolabs) was used at 1:10,000. Blots were imaged using ECL Prime Western Blotting Detection Reagent (GE Healthcare) and exposed on Amersham Hyperfilm (GE Healthcare).

41hORF cell line 1 was grown for 24 h in the presence of either DMSO or 100 µM 4EGI-1 (Cedarlane). Cells were treated with 100 µg/mL cycloheximide (CHX) (BioShop) for 10 min. Cells were harvested on ice by washing 2× with ice-cold PBS containing 100 µg/mL CHX, and lysing with 500 µL ice-cold filter-sterilized lysis buffer (10 mM Tris-HCl [pH 7.4], 5 mM MgCl2, 100 mM KCl, 1% Triton X-100, 2 mM DTT, 500 U/mL RNasin [Promega], 100 µg/mL CHX, Protease inhibitor [1× complete, EDTA-free, Roche]). Cells were scraped off the dish into tubes and sheared gently 4× with a 26-gauge needle. Lysed cells were centrifuged at 1300g for 10 min at 4°C, and clarified supernatant was isolated.

42A 10/50% sucrose gradient was created by combining heavy and light solutions on a BioComp Gradient Master. Heavy and light solutions consisted of 20 mM HEPES-KOH (pH 7.4), 5 mM MgCl2, 100 mM KCl, 2 mM DTT, 100 µg/mL CHX, and 20 U/mL SUPERaseIn, and either 10% or 50% sucrose (w/v), respectively. An amount of 300 µL of samples was layered on sucrose gradients and centrifuged in a precooled Beckman Ultracentrifuge L-90K using SW41 rotor at 36,000 RPM (221632.5G) for 2 h at 4°C. The gradient was fractionated using the BioComp Piston Gradient Fractionator and absorbance measurements were made using an Econo EM-1 UV Monitor (BioRad).

43hORF cell line 1 was grown for 1, 8, or 24 h in the presence of either DMSO, 100 μM 4EGI-1, or 5 μg/mL cycloheximide. Cells were pulsed with 1.5 μg/mL puromycin dihydrochloride (Thermo Fisher Scientific) for 10 min at 37°C. Cells were then harvested and lysed as above.

44Two sets of spike-in RNA were generated. An unlabeled S. cerevisiae spike-in is used to determine the enrichment of 4SU-labeled RNA over unlabeled RNA, as described previously. S. cerevisiae strain USY006 was grown in YPD liquid culture at 30°C, and RNA was isolated using a hot acidic phenol method. A 4SU-labeled D. melanogaster spike-in was also generated by supplementing S2 culture media with 100 µM 4SU for 24 h prior to harvesting. RNA was extracted using TRI-reagent (Molecular Research Center) as per manufacturer's instructions.

45Freshly thawed HEK293T hORFeome cell lines were cultured for three to four passages and seeded into DMEM + FBS culture media in 15 cm dishes such that they attained ∼50% confluence on the day of the time course. Media was replaced with DMEM + 10% FBS + 100 µM 4SU (Sigma-Aldrich) reconstituted in DMSO. Cells were harvested at 1, 2, 4, 8, 12, and 24 h after addition of 4SU. Harvesting was performed by dislodging cells off the plate during two washes with cold 1× PBS followed by spinning at 1000g for 5 min at 4°C. Cell pellets were resuspended in 1 mL TRI-Reagent (Molecular Research Center) and extracted according to manufacturer instructions.

46For translation inhibition experiments, hORF cell line 1 or cell line 4 cells growing in 10 mL DMEM + FBS in 10 cm dishes were pretreated with either 0.1% DMSO or 100 µM 4EGI-1 (Cedarlane) dissolved in DMSO for 1 h. Following this, 100 µM 4SU was added to media for all plates, and the time course was performed as described above.

47RNA was labeled as described previously. Briefly, 100 µg of total hORF RNA was mixed with 10 µg unlabeled S. cerevisiae RNA (i.e., 10% w/w) and 10 µg 4SU-labeled S2 D. melanogaster RNA (i.e., 10% w/w). Water was added to bring the volume up to 120 µL. An amount of 1 mg/mL HPDP-biotin (Thermo Fisher Scientific) was reconstituted in dimethylformamide by shaking at 37°C for 30 min at 300 RPM. 120 µL of 2.5× Citrate buffer (25 mM citrate, pH 4.5, 2.5 mL EDTA) and 60 µL of 1 mg/mL HPDP-biotin were added to the premixed RNA sample for each time point. This solution was incubated at 37°C for 2 h at 300 RPM on an Eppendorf ThermoMixer F1.5 in the dark. Samples were extracted twice with acid phenol, pH 4.5 (Invitrogen), and once with chloroform. RNA was precipitated with 18 µL 5M NaCl, 750 µL 100% ethanol, and 2 µL GlycoBlue (Invitrogen) overnight at −20°C. Precipitated RNA was pelleted for 30 min at 21,000g at 4°C. The RNA pellet was resuspended in 200 µL of 1× wash buffer (10 mM Tris-HCl, pH 7.4, 50 mM NaCl, 1 mM EDTA).

48Biotinylated RNA was purified using the µMACS Streptavidin microbeads system (Miltenyi Biotec). An amount of 50 µL Miltenyi beads per sample were preblocked with 48 µL 1× wash buffer and 2 µL yeast tRNA (Invitrogen), rotating for 20 min at room temperature. µMACS microcolumns were washed 1× with 100 µL nucleic acid equilibration buffer (Miltenyi Biotec), followed by 5× washes with 100 µL 1× wash buffer. Beads were applied to microcolumns in 100 µL aliquots, and again washed 5× with 100 µL 1× wash buffer. Beads were demagnetized and eluted off the column with 2× 100 µL 1× wash buffer, and columns were placed back on the magnetic stand. A total of 200 µL beads was mixed with each sample of biotinylated RNA and rotated at room temperature for 20 min.

49Samples were then applied to the microcolumns in 100 µL aliquots, washed 3× with 400 µL wash A buffer (10 mM Tris-HCl, pH 7.4, 6 M urea, 10 mM EDTA) prewarmed to 65°C, and then washed 3× with 400 µL wash B buffer (10 mM Tris-HCl, pH 7.4, 1 M NaCl, 10 mM EDTA). RNA was eluted with 5× 100 µL of 1× wash buffer supplemented with 0.1 M DTT, and flow through was collected in a tube. Purified RNA was precipitated with 30 µL 5 M NaCl, 2 µL GlycoBlue, and 1 mL 100% ethanol, incubated at −20°C overnight. Samples were then spun at 21,000g at 4°C for 30 min and resuspended in 20 µL water. RNA quality was assessed by running 3 µL of samples on a ∼1.5% agarose gel.

50A total of 10 µL of purified 4SU-labeled RNA or unpurified total RNA from the 24-h time point was used to prepare RNA-seq libraries using the TruSeq Stranded mRNA Sample Preparation Kit (Illumina), according to manufacturer's instructions. Adapter-ligated fragments were enriched with 14× PCR cycles. Approximately 16 to 22 samples were multiplexed on a single lane in an Illumina HiSeq 2500 at The Centre for Applied Genomics (The Hospital for Sick Children, University of Toronto) to obtain ∼10 million 50 bp single-end reads per sample. The data is available from GEO, accession GSE123165.

51Human (hg38), D. melanogaster (dm6), and S. cerevisiae (sacCer3) genomes were obtained in 2bit format using the UCSC Table Browser. 2bit files were converted to FASTA using the kentUtils command twoBitToFa, and GTF annotations were downloaded using the kentUtils command genePredToGtf. The three genomes were combined using custom bash scripts to make a hg38 + dm6 + sacCer3 genome.

52Library quality was assessed using FastQC v0.11.5 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Reads were trimmed and clipped for Illumina adapters using Trimmomatic v0.36 using the following settings: -phred33 ILLUMINACLIP: TruSeq3-SE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36.

53Trimmed reads were aligned to the indexed hg38 + dm6 + sacCer3 genome using STAR version 2.5.2 with the following nondefault settings: –outFilterMultimapNmax 10 –outFilterMismatchNoverLmax 0.05 –outFilterScoreMinOverLread 0.75 –outFilterMatchNminOverLread 0.85 –alignIntronMax 1 –outFilterIntronMotifs RemoveNoncanonical –outSAMtype BAM SortedByCoordinate –quantMode GeneCounts.

54HTSeq version 0.6.1 was then used to quantify gene counts from aligned BAM files using the following settings: –order=pos –stranded=reverse –minaqual=10 –mode=intersection-strict. Note that counting of features from human CDS and intronic GTF files were performed separately.

55Gene counts were loaded into RStudio version 3.3.1. Dplyr package (https://CRAN.R-project.org/package=dplyrnn) was used for all data manipulation and filtering. Steady state RNA-sequencing counts mapping to human CDS features were obtained from each cell line and normalized to library size to allow comparisons across samples. For a given cell line X, genes were described as “detectable hORF genes” if they met each of the following conditions: They were in the list of hORFs infected into cell line X;Normalized steady state RNA-sequencing for cell line X was greater than threefold that in the matched cell line;Normalized steady state RNA-sequencing for cell line X was greater than four reads.Genes that were not infected into cell line X were described as endogenous genes.

56They were in the list of hORFs infected into cell line X;

57Normalized steady state RNA-sequencing for cell line X was greater than threefold that in the matched cell line;

58Normalized steady state RNA-sequencing for cell line X was greater than four reads.

59All half-life calculations were performed in RStudio version 3.3.1 as described previously. Briefly, read counts for mature human mRNAs were filtered such that each gene had at least one read mapped to its CDS at each time point, and at least five reads mapped at any (at least one of six) time point. CDS-mapping reads for each gene at each time point were then normalized to the sum of all corresponding D. melanogaster mapping reads.

60Half-lives were calculated by fitting these normalized read counts at each time point to a bounded growth equation using weighted nonlinear least squares. The bounded growth equation has been previously described. Briefly, the equation states: where y(t) is the amount of a given transcript remaining at time t, yeq is the amount of that transcript at steady state, and k is the transcript-specific decay constant.

61The nls function in the stats package was used to fit the time points to the equation above, with settings equivalent to the following: start = c(yeq = max(y), k = −0.5)algorithm = “port”weights = 1/y(t)lower = c(yeq = 0, k = −Inf), upper = c(yeq = Inf, k = 0)

62start = c(yeq = max(y), k = −0.5)

63algorithm = “port”

64weights = 1/y(t)

65lower = c(yeq = 0, k = −Inf), upper = c(yeq = Inf, k = 0)

66If the data did not converge, a value of NA was returned. The half-life of each transcript is then obtained using the following equation:

67To measure local secondary structure, each gene's coding sequence was assayed by sliding 100 bp windows, each starting 3 bp apart. Each window was folded using ViennaRNA version 2.2.8 RNAfold function using default parameters. Minimum folding energy (MFE) for each 100 bp sequence was extracted from output files using custom bash scripts. Median and minimum MFEs across each CDS were determined using group_by and summarize functions in the dplyr package in RStudio.

68Codon usage frequency was calculated from each gene's coding sequence using the seqinr package's uco function. For frame shift controls, codon and amino acid usage were calculated after shifting the frame by +1 (removing positions 1, n − 2, and n − 1 from CDS of length n) and +2 (removing positions 1, 2, and n − 1 from CDS of length n).

69As described previously, codon stability coefficients (CSC) from a given half-life data set were calculated by determining the Spearman correlation between the codon frequency for each codon in a transcript with the measured half-lives of that transcript. Stop codons were excluded from CSC calculations. AASCs were calculated by taking the arithmetic mean of the CSC values for the corresponding codon(s).

70Ribosome profiling data sets for HeLa cells were obtained from GEO, accession GSE102113, and for HEK293T, U2OS, and NIH3T3 cells from GEO, accession GSE83616. To identify codon enrichments in the A, P, and E sites of translating ribosomes from ribosome profiling data, reads were first trimmed of their adapters using cutadapt v2.1. Trimmed reads were then aligned to the human genome (hg38, Gencode version 28) or the mouse genome (mm10, Gencode version 17) using STAR v2.5.2. Only uniquely aligned reads of length 28, 29, or 30 were then used for further analyses. The A, P, and E sites of these reads were defined as positions 17–19, 14–16, and 11–13, respectively. If these sites were not in-frame in a coding sequence, the read was ignored. If they were, the codons in the A, P, and E sites were recorded, allowing the calculation of the frequency with which each codon appeared within each site. This frequency was then compared to the null expected frequency. The null expected frequency for a codon was calculated as its frequency in the longest ORF of each gene weighted by the abundance of the gene in the ribosome profiling data set, using only 28–30 nt, in-frame reads. Amino acid pause scores were calculated by taking the arithmetic mean of the codon scores for the corresponding codon(s).

71Number of replicates, statistical tests used, and P-values are specified in the figures and figure legends.

72The accession number for the raw data files reported in this paper is GEO GSE123165. The accession number for the HeLa ribosome sequencing data used in this paper is GEO GSE102113. The accession number for the HEK293T, U2OS, and NIH3T3 ribosome sequencing data used in this paper is GEO GSE60426.

73Supplemental material is available for this article.