Cells alter their tRNA abundance to selectively regulate protein synthesis during stress conditions

Quantifying changes in tRNA abundance under diverse stress conditions

We first quantified changes in abundance of each of the 42 nuclear-encoded tRNAs during adaptation to different stress conditions in the yeast S. cerevisiae using reverse transcription quantitative polymerase chain reaction (RT-qPCR) (Fig. 1A; table S1; fig. S1, A and B; and data file S1). There are several approaches for quantifying tRNA abundance, which include tRNA microarrays (25), Northern blot (26), and sequencing-based methods (22). Each of these approaches has their strengths and limitations based on various considerations such as detection sensitivity, scalability, and the ability to resolve the identity of tRNAs and discriminate cleaved tRNA fragments from mature tRNA, as well as effects of tRNA modifications. Although the efficiency of reverse transcription may vary due to nucleotide modifications in tRNAs, this approach has been validated as a reliable method to quantify mature tRNAs (22, 27).

Measurement of the relative tRNA abundance profiles under four different stress conditions (oxidative stress, osmotic stress, temperature stress, and diauxic shift) at three different time points (20, 60, and 120 min) revealed that tRNA abundances changed substantially (about 2 to 5 log2 fold change; Fig. 1B) in a reproducible manner (fig. S2). Analysis of the changes in relative abundance of the individual tRNAs immediately upon stress revealed that decreasing the abundance of existing tRNA molecules (possibly through rapid degradation) could be a mechanism that changed relative tRNA abundance under different stress conditions (except under temperature stress) (Fig. 1C). The abundance of some tRNA molecules increased for all stress conditions except during oxidative stress, suggesting that active transcription could be another mechanism that regulates tRNA abundance during stress. At 120 min after stress, a higher proportion of the tRNA molecules showed decreased abundance, suggesting that repression of transcription- or degradation-based mechanisms might be a prevalent mechanism to alter tRNA abundance upon prolonged exposure to different stress conditions (Fig. 1C).

Patterns of changes in tRNA abundance during stress

Using t-distributed stochastic neighborhood embedding [t-SNE (28)] and K-means clustering, we analyzed the patterns of tRNA expression across the 12 conditions and time points and found that tRNAs can be segregated into three clusters (C1, C2, and C3) (Fig. 2A). The distribution of tRNAs among the three clusters corresponded to key functional features. C1 contained four of five tRNAs that are coded by single essential genes and six of seven tRNAs that are unique acceptors of their amino acid (Fig. 2A). Almost all tRNAs in C3 (7 of 8) and more than half of the tRNAs in C1 (9 of 16) contained nonoptimal anticodons (namely, those that can form a wobble codon-anticodon pair with low affinity). In contrast, C2 was depleted of tRNAs that carry nonoptimal anticodons (4 of 14). These observations indicate that under stress conditions, the tRNA pool is rearranged to adopt a complex structure that may influence translation.

On one hand, the abundance of C1 tRNAs either increased or remained stable during stress compared to normal conditions (Fig. 2B). This result is consistent with the notion that a decrease in the abundance of these essential tRNAs would likely negatively affect cellular fitness (21). On the other hand, the abundance of C2 tRNAs marginally decreased or remained stable and that of C3 tRNAs decreased under all stress conditions and at all time points (Fig. 2B). Because C1 and C3 tRNAs primarily decode nonoptimal codons, we reasoned that the differences in their abundance might have a role in protein production by controlling the rate at which transcripts with nonoptimal codons are translated. The tRNAs coding for Glu, Cys, and Gly showed reduced abundance under all stress conditions (Fig. 1B). These are the three amino acids that are required for the nonribosomal synthesis of the antioxidant glutathione (GSH), which is required for adaptation to stress conditions (29). One possible explanation is that a feedback mechanism between the nonribosomal GSH biosynthesis pathway and ribosomal protein translation could have evolved to ensure that the precursor amino acids are available in higher abundance for GSH production under stress conditions.

From a kinetic perspective, our results showed that after 20 min of exposure to stress, changes in tRNA abundances were most different within and across stress conditions, as indicated by their average correlation, which was calculated as the average Pearson correlation on the off-diagonal elements (Fig. 2C). After 60 min of exposure, we found a better correlation (Fig. 2C); however, the correlation was the highest for prolonged stress (t = 120 min) (Fig. 2C). The observed changes in tRNA fold change during stress suggested a biphasic behavior during adaptation to stress (fig. S3): An immediate transient response (at 20 min) with stress-specific variations, followed by a long-term adapted response (at 120 min) in which the tRNA pool is remodeled to a similar extent under all stress conditions but altered relative to the nonstress condition.

Genome-scale analysis of adaptation to tRNA abundance changes

Given that the tRNA abundance influences the decoding rate of codons during translation (7, 30), our observations implied that the rate of protein synthesis of individual genes might be selectively affected during stress. To investigate this notion, we first measured the extent to which the codon usage of all yeast genes was adapted to the tRNA abundance under normal and each of the four stress conditions. This was quantified by computing the tRNA adaptation index under normal conditions (n-tAI) and a newly developed metric, stress-adjusted tAI for each of the four stress conditions (s-tAI) (Fig. 3A, and data files S2 and S3), using the experimental measurements described above. For each stress condition, we identified genes that were better adapted, not affected, and less adapted to the experimentally measured tRNA pool by obtaining a z score of the rank change of their n-tAI and s-tAI values [Fig. 3A (orange, gray, and red data points) and data file S3].

We found that the genes with better codon adaptation to the tRNA pool under the different stress conditions were specifically enriched for functions related to external stimulus responses (Fig. 3A and fig. S4). Genes whose codon adaptation was not substantially altered were enriched in ribosomal proteins and translation machinery. Finally, the ones whose codon usage patterns were less well adapted were enriched in anabolic functions such as amino acid and lipid biosynthesis enzymes and tricarboxylic acid cycle (Fig. 3A, fig. S4, and data file S4). Although this observation suggested that the translation efficiency (TE) of certain genes might be selectively altered, we still observed a global positive correlation between tRNA supply and codon demand under all stress conditions, suggesting that overall translation was unlikely to be substantially compromised (fig. S5, A and B).

Alteration in translation rate and protein abundance during stress

To assess whether the genes that were identified as better adapted to the tRNA pool during stress, based on the changes in s-tAI ranks (Fig. 3A), showed a gain in TE, we analyzed experimentally derived ribosome footprinting data measured under oxidative stress (31, 32). The median log2 TE (or the amount of footprint normalized to mRNA abundance) of genes with better codon adaptation under oxidative stress was substantially higher compared to that of genes whose adaptation remains unaffected (Fig. 3B). Genes that were less well adapted had a lower median log2 TE compared to genes whose adaptation remained stable during oxidative stress (Fig. 3B).

Furthermore, investigation of protein abundance data during oxidative stress, as measured using confocal microscopy (33), revealed that genes that were better adapted to stress showed a significant increase in protein abundance compared to those that were less well adapted, although their mRNA abundance were comparable (Fig. 3C). These observations collectively suggest that changes in tRNA abundances might be an independent mechanism that can selectively increase the abundance of proteins encoded by certain transcripts that are better adapted (in terms of their codon usage) to the tRNA pool under stress condition.

Measuring protein production rates of designed sequences during stress

To experimentally determine whether changes in the tRNA pool can influence protein production rates, we designed monomeric enhanced green fluorescent protein (mEGFP) variants with identical amino acid sequence but different codon usage. Because most changes in the tRNA abundance occur immediately after stress (t = 20 min), we used the tRNA fold change data in this time window to design mEGFP variants. We examined each codon in the wild-type mEGFP sequence and substituted it with a synonymous codon if the fold change of the relevant isoacceptor tRNA increased more (or decreased less) than its native cognate tRNA (Fig. 3D and table S2). If the isoacceptor tRNA could recognize more than one codon, then the new codon was picked at random.

Because mRNA structure can affect translation (34, 35), we ensured that the designed transcripts did not contain any unusual secondary structure and had comparable free energies of folding using the ViennaRNA package (fig. S6). To avoid variations in mRNA abundance due to position effect, or due to plasmid copy number, we integrated each mEGFP variant in the TRP1 region of the yeast genome under a GAL1 promoter and obtained four distinct yeast strains that express the mEGFP gene (fig. S7). After confirming that all strains showed similar mEGFP mRNA levels (fig. S8), we measured the increase in fluorescence over time (as a proxy for protein production rate) of each strain subjected to the four stress conditions. We observed a positive relationship between the computed s-tAI values for the designed transcripts and the respective protein production rates (jointly and independently) for the different stress conditions (Fig. 3D). Although the relationship between these variables may be nonlinear, it was best modeled by a linear fit. This suggests that the changes in s-tAI affect the protein production rates of different mEGFP variants.

Together with the observations from the genomic analysis (Fig. 3, A to C), these results suggest that changes in tRNA abundance in response to stress might globally influence the proteome by selectively altering the protein production rates of transcripts under different stress conditions. Thus, transcripts with codons that are better adapted to the tRNA pool under a given stress condition are likely to be selectively translated with higher efficiency, leading to an increased abundance even after controlling for mRNA abundance and structure. Our data also show that the TE of the same transcript may be substantially altered when the tRNA pool is altered during stress, although transcript abundance remains the same.

Inferring the impact of s-tAI on protein production during stress

To characterize the possible impact of tRNA abundance changes on the S. cerevisiae proteome during stress, we used a stochastic model developed by Shah and co-workers (36) and simulated translation based on the experimentally measured changes in mRNA (37) and tRNA abundance (Fig. 4A and fig. S9). Consistent with our observation that stress-related genes showed better adaptation to the tRNA pool and increased protein abundance (Fig. 3C), we observed that the genes that are specifically up-regulated transcriptionally during the environmental stress response (ESR) showed a higher s-tAI and increased protein production rate in our simulations compared to other genes (Fig. 4B, left and middle). Thus, the combined effects of increased mRNA abundance and a higher s-tAI act in favor of increasing protein production rates of ESR up-regulated genes during stress. In line with the simulation results, the experimentally observed change in protein abundance (with respect to nonstressed condition) is significantly higher for the ESR genes during oxidative stress compared to the other genes (Fig. 4B, right).

We observed that a predicted increase in protein production rate was coupled with a slight decrease in overall steady-state protein abundance during oxidative stress (Fig. 4B). This effect may be because translation is globally reduced during oxidative stress (37), and most changes in tRNA abundance are due to reduced abundance rather than increased abundance (Fig. 1C). Despite the global reduction in translation, better adaptation of stress genes to the tRNA pool suggested that they were relatively rapidly translated compared to the other genes. The decrease in overall protein abundance could also be due to different gene sets used for the simulations and for which abundance measurements were available. Finally, the difference may also be explained by an increased degradation rate of proteins under stress, which might skew the overall protein abundance despite an increase in protein production rate.

To quantify the relative contribution of different factors such as changes in transcript abundance and codon adaptation to tRNA pools on the protein production rate, we performed an unbiased analysis using Bayesian statistics on the data obtained from the simulation. The most likely model that best explains the data showed that changes in protein abundance are highly dependent on changes in tRNA abundance, which contribute substantially to protein production rates, and are themselves explained by changes in tAI of the genes (Fig. 4C). In line with previous findings, regulation of translation initiation and changes in mRNA abundance also greatly influence protein production rate (34, 36, 37). Consistent with the observations from the experimental data (Fig. 3, C and D), we also found that changes in mRNA abundance were independent from tRNA-related effects (Fig. 4C), namely, changes in mRNA abundance and tRNA abundance are two distinct levers that regulate protein production. This suggests that changes in tRNA abundance can alter the relative translation rates of transcripts throughout the transcriptome and selectively influence the production rates, and hence abundances, of proteins during stress. Therefore, codon adaptation to the tRNA pool is an independent mechanism that can fine-tune protein abundance and complements other mechanisms such as transcriptional regulation, transcript degradation, and regulation of translation initiation.

Yeast strains and culture

All strains used in this work were based on Y03157 (BY4741; Mat α; his3Δ1; leu2Δ0; met15Δ0; ura3Δ0; YBR020w::kanMX4) obtained from Euroscarf deletion collection. Yeast strains were grown in YPD medium at 30°C (nonstress control). For stress conditions, changes were introduced as follows: (i) temperature was increased from 30° to 37°C (temperature stress); (ii) cells were gently centrifuged at 3000g for 2 min, and the medium was replaced by YPD containing 0.5 mM H2O2 (oxidative stress); (iii) cells were gently centrifuged at 3000g for 2 min, and the medium was replaced by YPD containing 1 M sorbitol (osmotic stress); and (iv) cells were gently centrifuged at 3000g for 2 min, and the medium was replaced by YPD containing 2% ethanol instead of 2% glucose as carbon source (diauxic shift). Cultures were followed for 120 min, and aliquots for further analysis were obtained at 0, 20, 60, and 120 min after stress.

RNA extraction

RNA was extracted from yeast cells using TRIzol phenol-chloroform extraction. Briefly, 5 ml of yeast culture were centrifuged at 3000g for 2 min. Yeast cells were resuspended in 150 μl of lysis buffer (0.1 M lithium acetate and 0.5% SDS) and heated at 70°C for 5 min. After that, 450 μl of TRIzol LS reagent was added and mixed for 15 s. Then, 150 μl of chloroform was added, mixed for 15 s, and incubated at room temperature for 5 min. Samples were centrifuged at 12,000g during 30 min, and the aqueous phase was obtained. RNA was then precipitated with 450 μl of isopropanol. The RNA pellet was recovered by centrifugation at 10,000g for 20 min and then cleaned twice with 75% ethanol. Finally, RNA was dissolved in 25 μl of sterile RNase-free water (Life Technologies) and immediately used for cDNA synthesis. The samples were quantified using a NanoDrop ND-1000 (Thermo Fisher Scientific). All samples had an A260/A230 ratio of >2 and A260/A280 ratio of >2.

cDNA synthesis

To ensure that the isolated RNA was free from DNA, samples were incubated with RNase-free DNase I (Sigma-Aldrich) for 30 min according to the manufacturer’s instructions. The reaction was stopped by adding 5 mM EDTA (Sigma-Aldrich) and heating the sample at 70°C for 10 min. Then, 1 μg of RNA was reverse-transcribed using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific). To minimize the effect of secondary structure, the reaction was conducted at 60°C (instead of 45°C) and extended for 30 min (instead of 15 min) to reduce the effect of transcription pauses due to tRNA modifications.

Quantification of individual tRNAs

Individual tRNAs were quantified by qPCR on an ECO Real-Time qPCR thermocycler (Illumina). One microliter of a 1:1000 dilution of cDNA from the reverse transcription reaction and 175 μM of the corresponding primers were used for the PCR reaction. Primers were designed using Primer-BLAST (https://ncbi.nlm.nih.gov/tools/primer-blast/) to ensure that each sequence is specific for the template and does not hybridize with any other region of the S. cerevisiae genome (table S1 and fig. S1, A and B). The PCR reaction was conducted using the SYBR Green PCR Master Mix (Life Technologies), and the results obtained were processed using the ΔΔCt method. The geometric mean of four housekeeping genes (ALG9, TAF10, TFC1, and UBC6) was used as reference (53, 54). No bias was observed in qPCR efficiency for the 42 primer pairs designed (table S1 and fig. S1, A and B). All quantifications were done using biological triplicates.

Design of mEGFP genes

Four mEGFP sequences were designed to test the effect of codon usage and altered tRNA abundance on protein translation. The sequences were designed using the Visual Gene Developer software (55) using the following criteria: All codons were substituted by the most up-regulated codon under the particular stress condition. In those cases where the same tRNA translated two codons, we randomly chose the codon. We also tested that all transcripts derived from our genes have similar RNA folding energy (fig. S6). Modified mEGFP sequences are in table S2.

Construction of GFP genes and cloning into yeast

A pMA synthetic vector carrying a multicloning site flanked by 60– base pair regions that is homologous to TRP1 was purchased from GeneArt (Life Technologies). TRP1 gene was chosen as the position for insertion due to its close proximity to the centromere and the high expression levels of genes in this region; these properties will likely minimize the frequency of gene silencing and will ensure a steady-state expression level through several generations. All GFP versions were synthesized de novo (GeneArt, Life Technologies) and inserted in the pMA vector between pGal1 and TEF terminator using the In-Fusion HD Cloning Kit (Clontech). A plasmid scheme is displayed in fig. S7. The constructs were propagated in Escherichia coli and purified using a QIAprep kit (Qiagen). The insertion sequence was amplified by PCR, and the product was purified by agarose gel electrophoresis using a MinElute kit (Qiagen). Finally, the gene was introduced in the S. cerevisiae genome by homologous recombination using the lithium acetate protocol (56). To ensure that the gene was inserted in the correct position, it was amplified using PCR after recombination and sequenced.

mRNA expression levels

Fresh cells (5 ml) were grown with 1% galactose, centrifuged and suspended with 1 ml of 0.2 M lithium acetate and 1% SDS solution, and boiled at 70°C for 5 min. Then, RNA was extracted as detailed before. qPCR was performed using the SYBR Green PCR Master Mix (Life Technologies) as described before using the primers listed in table S3.

Fluorescence measurements

S. cerevisiae cells expressing the four modified mEGFP constructs were grown in a plate reader (Tecan) monitoring both mEGFP fluorescence and cell growth [optical density at 600 nm (OD600)] over time. Fluorescence measurements were normalized by cell number (measured as fluorescence/OD600) and corrected for autofluorescence. We also measured the transcript abundance for the different sequences and did not observe any difference in terms of the mRNA expression levels (fig. S8). The slope of normalized fluorescence over 3 hours (▲A.U./s) was taken as a measure of protein production.

Computational analysis and stochastic simulations

Details on computational analysis, simulation, and statistical evaluation of the data presented can be found in the Supplementary Materials.