Highly Predictive Reprogramming of tRNA Modifications Is Linked to Selective Expression of Codon-Biased Genes.

Cells respond to stress by controlling gene expression at several levels, with little known about the role of translation. Here, we demonstrate a coordinated translational stress response system involving stress-specific reprogramming of tRNA wobble modifications that leads to selective translation of codon-biased mRNAs representing different classes of critical response proteins. In budding yeast exposed to four oxidants and five alkylating agents, tRNA modification patterns accurately distinguished among chemically similar stressors, with 14 modified ribonucleosides forming the basis for a data-driven model that predicts toxicant chemistry with >80% sensitivity and specificity. tRNA modification subpatterns also distinguish SN1 from SN2 alkylating agents, with SN2-induced increases in m(3)C in tRNA mechanistically linked to selective translation of threonine-rich membrane proteins from genes enriched with ACC and ACT degenerate codons for threonine. These results establish tRNA modifications as predictive biomarkers of exposure and illustrate a novel regulatory mechanism for translational control of cell stress response.

Introduction

Cells respond to environmental stressors and xenobiotic exposures by controlling gene expression with multilayered, complex regulatory networks. The emergence of tools to quantify the molecular changes associated with these networks has led to an appreciation for the mechanistic insights and predictive power derived from transcriptional profiling, proteomics, and metabolomics.[1−5] Using a novel bioanalytical platform, we recently described an example of another network that regulates the cellular response to xenobiotic exposures and other stresses: stress-induced reprogramming of a system of two-dozen post-transcriptional modifications on tRNA (tRNA), which promotes selective translation of codon-biased mRNAs for critical response proteins.[6−8] Most forms of RNA contain modified ribonucleosides, with more than 120 different chemical structures across all organisms and 2–3 dozen in any one organism.[9,10] tRNA is the most extensively modified RNA species, with the presence of specific ribonucleoside structures affecting the rate and fidelity of translation,[11,12] tRNA stability,[13,14] cellular stress responses,[6,15,16] and cell growth.[17] To define a systems-level behavior for tRNA modifications, we developed a chromatography-coupled tandem quadrupole mass spectrometry (LC-MS/MS) method to quantify stress-induced changes in the system of modified ribonucleosides in a population of cellular tRNA molecules or in individual, purified tRNA species.[7,8,18−20] We applied this approach to budding yeast exposed to four mechanistically different toxicants (hydrogen peroxide, H2O2; methylmethanesulfonate, MMS; sodium arsenite NaAsO3; and sodium hypochlorite, NaOCl) and observed agent-specific changes in the relative quantities of 23 modified ribonucleosides in tRNA.[7] These results pointed to the spectrum of modified ribonucleosides in tRNA as a predictive biomarker of exposure. Furthermore, stress-specific reprogramming of tRNA wobble modifications was found to result in selective translation of mRNAs containing biased use of the codons corresponding to these tRNAs, with the mRNAs representing critical stress response proteins.[6−8]

These observations were made with four chemically distinct toxicants, which raised questions about the specificity of the tRNA reprogramming signatures for classes of chemical toxicants (e.g., all oxidants) and the generality of the translational response mechanism. To test the hypothesis that cells respond to mechanistically similar toxicants with similar patterns of tRNA modification reprogramming and that the patterns are predictive of the chemical class of the toxicant, we analyzed changes in the levels of 23 modified ribonucleosides in tRNA from Saccharomyces cerevisiae exposed to four oxidizing agents (hydrogen peroxide (H2O2), t-butyl hydroperoxide (TBHP), γ-radiation (γ-rad), and peroxynitrite (ONOO–)) and five alkylating agents (methylmethanesulfonate (MMS), ethylmethanesulfonate (EMS), isopropyl methanesulfonate (IMS), N-methyl-N′-nitro-N-nitrosoguanidine (MNNG), and N-nitro-N-methylurea (NMU)). Multivariate statistical analysis and data-driven modeling of the RNA modification patterns proved to be highly predictive of toxicant chemistry, with bioinformatic and proteomic analyses establishing a mechanistic link among exposure, tRNA modifications, and codon-biased translation of response proteins. The results support a 30 year old model for one of the mechanisms by which SN2 alkylating agents promote cell death, by damaging membrane proteins, and establish a novel regulatory mechanism linking exposure to translational response in cells, which is summarized in Figure 1.

Figure 1

Proposed model for stress-specific reprogramming of tRNA modifications leading to selective translation of codon-biased response proteins. The specific model shown here is based on exposure to the SN2 alkylating agent, methylmethanesulfonate (MMS), which leads to upregulation of the modified ribonucleoside m3C in tRNA species coding for threonine and serine, with their cognate codons enriched (ACC, ACT) in membrane-associated proteins. This model is consistent with broad mechanisms for MMS-induced toxicity proposed by Smith and Grisham.[50]

Experimental Procedures

Materials

All chemicals and reagents were of the highest purity available and were used without further purification. MMS, EMS, NMU, TBHP, H2O2, RNase A, and alkaline phosphatase were obtained from Sigma Chemical Co. (St. Louis, MO). IMS was obtained from Pfaltz & Bauer, Inc. (Waterbury, CT). MNNG was obtained from TCI America (Portland, OR). Sodium peroxynitrite was obtained from Cayman Chemical Co. (Ann Arbor, MI). Nuclease P1 was obtained from Roche Diagnostic Corp. (Indianapolis, IN). Phosphodiesterase I was obtained from USB (Cleveland, OH). Yeast extract and peptone were obtained from Biomed Diagnostics, Inc. (White City, OR). Micron YM10 filters were obtained from PALL Corp. (Port Washington, NY). HPLC-grade water and acetonitrile were obtained from Mallinckrodt Baker (Phillipsburg, NJ). S. cerevisiae BY4741 was obtained from American Type Culture Collections (Manassas, VA).

S. cerevisiae Cytotoxicity Dose–Response Assays

S. cerevisiae BY4741 was cultured in YPD (yeast extract-peptone-dextrose) media with 200 μg/mL of geneticine at 30 °C with shaking at 220 rpm. Each culture was grown to mid log phase (OD660 ∼ 0.6–0.8) followed by exposure to the following toxicants and concentrations: H2O2 (0, 2, 3.5, 5, 10, 15, and 20 mM), TBHP (0, 0.7, 2, 4, 7, 14, 22, 25, and 29 mM), ONOO– (0, 0.3, 0.5, 0.8, 1.0, 1.5, and 2.0 mM), γ-rad (116 Gy/min) (0, 21.3, 168, 327, 513, and 606 Gy), MMS (0, 1.2, 6, 12, 24, 36, and 48 mM), EMS (0, 0.19, 0.29, 0.39, 0.49, and 0.58 M), IMS (0, 8, 17, 33, 50, and 66 mM), MNNG (0, 41, 61, 82, 102, and 136 mM), and NMU (0, 1.3, 2.3, 3.2, and 4.2 mM). After 1 h at 30 °C, these cultures were diluted 104-fold with YPD medium, and 50 μL was plated on YPD agar. Survival rates of exposed cells were determined by comparing colony counts for untreated and treated samples after 2 days of growth.

Cell Exposure and tRNA Isolation

Cultures of S. cerevisiae at mid log phase (OD660 ∼ 0.6; ∼2 × 107 cells) were treated with 5 mM H2O2, 25 mM TBHP, 0.8 mM ONOO–, 500 Gy γ-rad (116 Gy/min), 24 mM MMS, 190 mM EMS, 50 mM IMS, 82 mM MNNG, or 3.2 mM NMU and incubated for 1 h at 30 °C alongside untreated control cultures. Cells were pelleted by centrifugation and resuspended in TRIzol reagent with the addition of antioxidants (0.1 mM desferrioxamine and 0.1 mM butylated hydroxytoluene) and deaminase inhibitors (5 μg/mL coformycin and 50 μg/mL tetrahydrouridine). Cells in this solution were lysed by three cycles of bead beating in a Thermo FP120 bead beater at 6.5 m/s for 20 s each cycle, with 1 min of cooling on ice between cycles. The lysates was mixed with one-volume of chloroform, and phases were separated by centrifugation. Small RNA in the aqueous phase was isolated using using the PureLink miRNA isolation kit according to manufacturer’s instructions. Each sample yielded ∼6 μg of RNA composed of ∼90% tRNA, as judged by Bioanalyzer analysis (Agilent Corporation).

Quantification of tRNA Modifications

tRNA modifications were quantified using an established LC-MS/MS method,[7] which starts with hydrolysis of tRNA (6 μg) in a solution (50 μL; pH 6.8) containing 30 mM sodium acetate, 2 mM ZnCl2, 0.02 unit/μL of nuclease P1, 0.1 unit/μL of RNase A, 5 μg/mL coformycin, 50 μg/mL tetrahydrouridine, 0.1 mM deferoxamine mesylate, 0.1 mM butylated hydroxytoluene, and 6 pmol of [15N]5-2′-deoxyadenosine ([15N]5-dA) as an internal standard. Following incubation at 37 °C for 3 h, the solution was augmented (additional 50 μL volume) with an additional 30 mM sodium acetate, 0.2 unit/μL of alkaline phosphatase, and 0.01 unit/μL of phosphodiesterase I at a final pH of 7.8 and incubated 37 °C overnight. Proteins were removed by centrifugal ultrafiltration using a Microcon YM-10 filter. An aliquot of filtrate containing ∼0.4 μg of ribonucleosides was loaded on a Thermo Scientific Hypersil GOLD aQ reverse-phase HPLC column (150 × 2.1 mm, 3 μm particle size), and ribonucleosides were eluted with the following gradient of acetonitrile in 8 mM ammonium acetate at a flow rate of 0.3 mL/min at 36 °C: 0–18 min, 0%; 18–23 min, 0–1%; 23–28 min, 1–6%; 28–30 min, 6%; 30–40 min, 6–100%; 40–50 min, 100%. The HPLC column was coupled to an Agilent 6410 Triple Quadrupole LC/MS mass spectrometer with an electrospray ionization source operated in positive ion mode with the following parameters: gas temperature, 350 °C; gas flow, 10 L/min; nebulizer, 20 psi; and capillary voltage, 3500 V. The first and third quadrupoles (Q1, Q3) were fixed to unit resolution, and the modifications were quantified by predetermined molecular transitions. Q1 was set to transmit the parent ribonucleoside ions, and Q3 was set to monitor the deglycosylated product ions, except for pseudouridine (Y), which was detected by setting Q1 to transmit the parent ion and Q3 set to monitor the m/z = 125 product ion resulting from pyrimidine ring fragmentation.[7,21] The dwell time for each ribonucleoside was 200 ms. The retention time, m/z of the transmitted parent ion, m/z of the monitored product ion, fragmentor voltage, and collision energy of each modified nucleoside and 15N-labeled internal standard are as follows: D, 2.2 min, m/z 247 → 115, 80 V, 5 V; Y, 2.3 min, m/z 245 → 125, 80 V, 10 V; m5C, 5.4 min, m/z 258 → 126, 80 V, 8 V; Cm, 6.4 min, m/z 258 → 112, 80 V, 8 V; m5U, 7.9 min, m/z 259 → 127, 90 V, 7 V; ncm5U, 8.7 min, m/z 302 → 170, 90 V, 7 V; ac4C, 19.0 min, m/z 286 → 154, 80 V, 6 V; m3C, 5.0 min, m/z 258 → 126, 80 V, 8 V; Um, 10.7 min, m/z 259 → 113, 90 V, 7 V; m7G, 8.5 min, m/z 298 → 166, 90 V, 10 V; m1A, 6.9 min, m/z 282 → 150, 100 V, 16 V; mcm5U, 14.6 min, m/z 317 → 185, 90 V, 7 V; m1I, 16.0 min, m/z 283 → 151, 80 V, 10 V; Gm, 17.2 min, m/z 298 → 152, 80 V, 7 V; m1G, 18.8 min, m/z 298 → 166, 90 V, 10 V; m2G, 22.2 min, m/z 298 → 166, 90 V, 10 V; I, 7.8 min, m/z 269 → 137, 80 V, 10 V; mcm5s2U, 31.3 min, m/z 333 → 201, 90 V, 7 V; [15N]5-dA, 30.0 min, m/z 257 → 141, 90 V, 10 V; m22G, 31.7 min, m/z 312 → 180, 100 V, 8 V; t6A, 32.8 min, m/z 413 → 281, 100 V, 8 V; Am, 33.1 min, m/z 282 → 136, 100 V, 15 V; yW, 34.1 min, m/z 509 → 377, 120 V, 10 V, and i6A, 34.5 min, m/z 336 → 204, 120 V, 17 V. The signal from each modified nucleoside was normalized by dividing by the signal from [15N]5-dA for the purpose of comparison between samples. Complete analysis was performed as technical duplicates of each of five biological replicates, with samples from each biological replicate analyzed as a single batch to minimize the influence of day-to-day variation in instrument performance.

Data Analysis

Exposure-induced changes in the quantities of modified ribonucleosides were calculated as fold-change values using the normalized MS signal intensities for the treated and unexposed samples. This was accomplished by averaging the normalized MS signal intensity for each ribonucleoside in the three control samples and then dividing the corresponding normalized MS signal intensity for each ribonucleoside in each treated sample by this average control value. These fold-change values are presented in Table S1, and the mean (±standard deviation) fold-change values for the five replicates are presented in Table S2. The fold-change values were then analyzed by hierarchical clustering using the centroid linkage algorithm in Cluster 3.0 following log2 transformation of the fold-change data; heat map representations were produced using Java Treeview, as described elsewhere.[7] The heat map for the entire data set is shown in Figure S3, and that for the averaged data, in Figure 3.

Figure 3

Hierarchical clustering analysis of average fold-change values for tRNA modifications in total tRNA from cells exposed to equitoxic (LD80) doses of different alkylating agents and oxidizing agents. On the basis of the dose–response curves shown in Figure S1, budding yeast cells were exposed to LD80 doses of the various agents, and tRNA modifications were quantified by LC-MS/MS. The fold-change values (Table S2, with significance testing) were derived from the average of normalized MS signal intensity data from five biological replicates (Figure S2 and Table S1), and hierarchical clustering analysis was performed in log space (log2).

Data-Driven Modeling

To assess the predictive power of the stress-altered tRNA modification patterns, the normalized MS signal intensities for each biological replicate of treated and control cells was used to develop a data-driven model (Figure 4A). The data for individual samples were labeled as one of three classes of toxicant exposure: unexposed control (CT), alkylating agent-exposed (AA), and oxidizing agent-exposed (OX). A classification model was then developed using K-nearest neighbor classification. From each exposure class, normalized ribonucleoside signal intensities were compared to the other classes using multiple t-tests with Bonferroni correction, and those with p values < 0.01 were assigned as unique features of that exposure group. All of these unique features were then set as parameters to construct a data-driven model using the programming software R, in which all data were randomly assigned into two groups: a training set to build the model and a testing set to evaluate the prediction accuracy of the model. For this evaluation, the confusion matrix method (Figure 4B) was used to determine prediction sensitivity and prediction specificity.

Figure 4

A data-driven model to define the ability of tRNA modification spectra to distinguish different stresses. K-nearest neighbor classification was used to develop a data-driven model (A) for the normalized MS signal intensities for each biological replicate of treated and control cells (Table S1; hierarchical clustering analysis of data sets shown in Figure S3). A confusion matrix (B) shows the number of false positives, false negatives, true positives, and true negatives that were used to determine prediction sensitivity (C) and prediction specificity (D) of the model, based on a total of 20 training and test cycles. A low standard error (<2%) demonstrates the stability of the data-driven model.

SILAC Proteomics

MMS-induced changes in protein levels were determined by SILAC proteomics.[8,22] To prepare isotopically labeled proteins, Lys1Δ yeast cells were grown in yeast nitrogen base (YNB) liquid medium containing 30 mg/L of l-lysine-U-[13C]6, [15N]2 (Isotec-SIGMA, Miamisburg, OH) for at least 10 generations, until they reached log-phase (OD600 ∼ 0.7).[23] Wild-type yeast cells were grown to log phase (OD600 ∼ 0.7) in YNB medium containing 30 mg/L of l-lysine. Cells were harvested by centrifugation at 1500g for 10 min at 4 °C and washed twice with ice-cold water. Cells were then lysed by suspension in alkaline buffer (2 M NaOH, 8% 2-mercaptoethanol v/v), and proteins were precipitated by adding 50% TCA and incubating on ice for 10 min, followed by centrifugation at 15 000g for 15 min at 4 °C. The pellet was resuspended in lysis buffer (8 M urea, 75 mM NaCl, 50 mM Tris, pH 8.2, 50 mM NaF, 50 mM β-glycerophosphate, 1 mM sodium orthovanadate, 10 mM sodium pyrophosphate, 1 mM phenylmethylsulfonyl fluoride).[24] Protein concentration was determined by the Bradford assay.[25] Heavy SILAC-labeled lys1Δ yeast proteins were used as a global internal standard.[26] For the MMS exposure studies, cells (OD600 ∼ 0.7) were treated for 1 h with 0.0125% MMS in YNB media at 30 °C, which caused less than 5% cell death based on a cytotoxicity assay. Proteins were then harvested as noted earlier.

For the proteomic analyses, internal standard was added to all treated and untreated protein samples (1:1), and the protein mixture was subjected to disulfide reduction by incubation for 2.5 h at 37 °C in 1 mM dithiothreitol, followed by thiol alkylation by 5.5 mM iodoacetamide for 40 min at ambient temperature in the dark. Proteins were then digested with 50:1 (w/w) trypsin overnight at 37 °C. Peptide mixtures were loaded onto a Vydac C18 trap column (150 μm × 10 mm; 5 μm diameter, 300 Å pore-size particle; Grace, Deerfield, IL) at flow rate of 5 μL/min and eluted onto a Vydac C18 analytical column (75 μm × 150 mm, 5 μm diameter, 300 Å pore-size particle) at 200 nL/min with a gradient of 2–98% acetonitrile with 0.1% formic acid over 180 min. Eluted peptides were analyzed by mass spectrometric analysis on a QSTAR-XL system (Applied Biosystems/MDS Sciex, Foster City, CA). Data were obtained from three biological replicates. Acquired MS/MS spectra were parsed by Spectrum Mill (Agilent Technologies, Foster City, CA) and searched against Saccharomyces Genome Database (SGD). SILAC peptide and protein quantitation was performed with differential expression quantitation, and SILAC protein ratios were determined as the average of all peptide ratios assigned to the protein. Differential protein expression was determined by Student’s t test between samples.

A total of 2381 high-confidence proteins with a false-discovery rate of 0.5% were identified in control and MMS-treated samples.[22] To control for protein expression regulated by transcription, rather than by translational control mechanisms, the protein expression data were corrected using microarray data from our previous studies with the same yeast strain and MMS treatment conditions.[6,27,28] After removing genes that show the same expression change (up, down, or no change) for both mRNA and protein, we identified 222 upregulated and 438 downregulated genes whose expression is significantly regulated by translational machinery (Table S3).

Gene Ontology Annotation

Gene functional categorization and pathway analysis were performed with DAVID Bioinformatics Resources 2011.[29] The annotated proteins are clustered according to the biological process branch of the Gene Ontology (GO) annotation. The statistical significance of over- or under-representation of proteins in each GO category was assessed using a hypergeometric distribution, and the significance indicated by the p values for each GO category.[29]

Results

Establishing Equivalent Stress Conditions for Oxidizing and Alkylating Agents

A quantitatively meaningful comparative analysis of stress responses requires that cells experience equivalent levels of stress. To this end, we chose cytotoxicity as a common end point, with doses producing 80% lethality for a 1 h exposure to the four oxidizing agents (ONOO–, H2O2, γ-rad, TBHP) and five alkylating agents (MMS, EMS, IMS, MMNG, MNU), the structures of which are shown in Figure 2. The dose–response curves for the nine agents are shown in Figure S1, from which we determined LD80 doses of 5 mM H2O2, 25 mM TBHP, 0.8 mM ONOO–, 510 G γ-radiation (116 Gy/min), 24 mM MMS, 190 mM EMS, 50 mM IMS, 82 mM MNNG, and 3 mM NMU. To assess the role of Trm140, which methylates C3 of cytosine at position 32 in several tRNAs (see below), in the cellular response to alkylating agent exposures, cell survival was determined for wild-type and Δtrm140 mutant S. cerevisiae exposed to the alkylating agents at LD80 doses determined for wild-type cells. As shown in Figure S3, loss of Trm140 modestly but significantly increased the sensitivity of the cells to SN2 alkylating agents MMS and EMS but not the other three SN1 alkylating agents.

Figure 2

Structures of the oxidizing and alkylating agents used in the studies.

Quantifying Exposure-Induced Changes in the Full Set of tRNA Modifications

Having defined equivalent exposure conditions, the effect of exposure to nine different agents on the relative quantities of 23 tRNA modifications was determined using LC-MS/MS. The analysis consisted of five biological replicate sets of samples each comprised of three unexposed controls and one exposure for each of the nine agents. Normalized mass spectrometric (MS) signal intensities for each modification in the three unexposed controls were averaged, and changes in the MS signal intensities for the exposed samples were determined relative to this control value. Table S1 shows the complete set of 1380 data elements for the five biological replicates of this analysis, with the mean fold-change for each ribonucleoside and each exposure shown in Table S2. The results reveal that the levels of 22 of the 23 modified ribonucleosides were changed significantly by the various exposures (p < 0.05 by Student’s t test); only m5U did not change significantly (Table S2).

Stress Induces Reprogramming of tRNA Modifications with Unique Class- and Agent-Specific Patterns

To identify exposure-dependent patterns in the relative quantities of tRNA modifications in cells exposed to the nine toxicants, the fold-change data in Tables S1 and S2 were subjected to hierarchical clustering analysis, with the heat map representations shown in Figures 3 and S3. Analysis of the averaged data sets in Figure 3 shows clear class-specific distinctions for the groups of oxidizing agents and alkylating agents. All alkylating agents caused significant (p < 0.05) increases in Um, m2G, mcm5s2U, mcm5U, and m1A, but these modifications were not significantly affected by the oxidizing agents. Furthermore, levels of yW and m1I decreased significantly in all alkylator-exposed cells, but the changes in oxidant-exposed cells were not significant. In contrast, the levels of ncm5U, m5C, and i6A did not change in cells exposed to any of the alkylating agents, but they increased significantly in response to all oxidizing agents. In addition to these class-specific signature patterns of modified ribonucleosides, there were also “sub-signatures” apparent for both classes of toxicant. This is most apparent for the alkylating agents, for which SN2 alkylating agents[30] (EMS, MMS) cause increases in m3C and m7G, whereas SN1 alkylating agents[30] (IMS, MNNG, MNU) uniquely increase Am (Figure 3).

A Data-Driven Model for Distinguishing Alkylation and Oxidation Stresses

The class- and agent-specific distinctions revealed by clustering analysis suggested that the stress-specific patterns of tRNA modifications could be used to predict exposure chemistry. To test this hypothesis, we developed a data-driven model (Figure 4A) using the K-nearest neighbor classification method to classify and predict the agent classes in term of tRNA modification spectra. Following 20 training and test cycles of supervised learning, a stable model (standard error <2%) was established. As shown in the confusion matrix in Figure 4B, which reports the average number of false positives, false negatives, true positives, and true negatives, the model proved to have sensitivities of 95% for the alkylating agent-exposed group (AA), 94% for oxidant exposures (OX), and 78% for unexposed cells (CT). The predictive specificities are 95% for AA, 76% for OX, and 98% for CT (Figure 4C). On the basis of this model, a set of 14 modified ribonucleosides was identified as contributing the most to distinguishing the exposures, with the group for predicting alkylating agents composed of m3C, m7G, yW, mcm5U, Am, Gm, m5C, mcm5s2U, and m1I, and those for identifying oxidizing agent-exposures consisting of m5C, m3C, m7G, Gm, ncm5U, m22G, i6A, yW, and Cm (Table 1). These features match those revealed in the hierarchical clustering analysis (Figure 3).

Table 1

tRNA Modifications That Contribute Most Significantly to the Data-Driven Model for Predicting Exposure Chemistrya

rN1	tRNA	position	modifying enzyme	rN1	tRNA	position	modifying enzyme
Modifications Increased by Alkylating Agents				Modifications Increased by Oxidizing Agents
m3C	multiple	32, e2	Trm140	m5C	tRNALeuCAA	34	Trm4
m7G	multiple	46	Trm8, Trm82		tRNAPheGAA	40	Trm4
mcm5U	tRNAArgUCU	34	Trm9, Elp1-6, Kti11-13		multiple	48	Trm4
mcm5s2U	tRNAGluUUC	34	Trm9, Nfs1, Elp1-6, KTI11-13		multiple	49	Trm4
Am	tRNAHisGUG	4	unknown	ncm5U	tRNAValUAC	34	Elp1-6, Kti11-13
Gm	multiple	18	Trm3	m22G	multiple	26	Trm1
	tRNAPheGAA	34	Trm7	i6A	multiple	37	Mod5
m2G	multiple	10	Trm11	Cm	multiple	32	Trm7
	tRNAValCAC	26	unknown		tRNATrpCCA	34	Trm7
Modifications Decreased by Alkylating Agents					multiple	4	unknown
yW	tRNAPheGAA	37	Trm5
m1I	tRNAAlaIGC	37	Tad1, Trm5

The modified ribonucleosides (rN) are organized into three groups, depending on whether their levels increased or decreased in response to the two classes of chemical stress (alkylating agents; oxidizing agents). Specific tRNAs are noted if the modification is located in a single tRNA species. Information was from Modomics Database,[10] Saccharomyces Genome Database,[43] and D’Silva et al.[57] The link between m3C and cell response to SN2 alkylating agents MMS and EMS is shown in Figure S3, in which it is shown that loss of Trm140 confers sensitivity these agents.

Linking tRNA Modification Patterns to Selective Translation: Analysis of MMS-Induced Changes in the Proteome

The observation that SN2 alkylating agents caused differential upregulation of m3C in tRNA is reminiscent of H2O2-induced increases in m5C (Figure 3) that led to selective translation of mRNAs enriched in the Leu codon TTG.[8] In addition, the increase in mcm5U following MMS exposure has been linked to Trm9 and to the selective translation of mRNAs enriched with the Arg codon AGA.[6,19] Following this line of logic with MMS-treated cells, m3C is present at position 32 of the anticodon loop of tRNAThr(IGU), tRNASer(UGA), and tRNASer(CGA) in S. cerevisiae, with these tRNAs reading codons ACT, TCA and TCG, respectively.[31] Position 32 is adjacent to the anticodon (positions 34–36) and has the potential to modulate codon–anticodon interactions. It was thus hypothesized that MMS treatment would cause upregulation of proteins containing high levels of Thr and Ser. To test this hypothesis, we performed a comparative analysis of MMS-altered protein expression in S. cerevisiae using SILAC proteomics to quantify soluble proteins.[22] This analysis yielded 2381 high-confidence proteins with a false-discovery rate of 0.5%.[22] To control for protein expression regulated by transcription, rather than by translational control mechanisms, we corrected the protein expression data using microarray data from our previous studies with the same yeast strain and MMS treatment conditions.[6,27,28] After removing genes that show the same expression change (up, down, or no change) for both mRNA and protein, we identified 222 upregulated and 438 downregulated genes whose expression is significantly regulated by translational machinery (Table S3). This corrected index of protein expression was then evaluated for codon usage to test the hypothesis that MMS alters translational efficiency of genes with biased use of codons related to m3C. To minimize the possibility that the influence of m3C is counterbalanced by other functional constraints, we considered only proteins with high usage of optimal codons (each being >3% of the total codons in the gene for the protein). Among the Thr codons regulated by tRNAThr(IGU) (ThrACA, ThrACC, ThrACG, ThrACT) and Ser codons regulated by tRNASer(UGA) and tRNASer(CGA) (SerTCA, SerTCG), only ThrACC and ThrACT are optimal codons, which are more efficiently translated and more frequently used by yeast genome in comparison with other synonymous codons (ThrACA, ThrACG).[32] As a result, other codons are rarely used, and their influence is therefore limited. Indeed, less than 20 of 660 MMS-altered proteins identified in this study have significantly high usage of these nonoptimal codons (hypergeometric distribution, p < 0.01; yeast genome as background). However, as shown in Figure 5A,B, genes upregulated by MMS treatment are significantly skewed to higher usage of ThrACC and ThrACT codons in comparison with downregulated genes (Student’s t test; ThrACC, p = 3.8 × 10–09; ThrACT, p = 8.0 × 10–7). This is also apparent as a shift in the distribution of proteins with different ACC and ACT codon content (Figures 5C,D; Kolmogorov–Smirnov test: ThrACC, p = 3.7 × 10–06, ThrACT, 1.9 × 10–5). The ThrACC codon is significantly enriched in 13% of upregulated genes (29 of 222; hypergeometric distribution, p < 0.01, with yeast genome as background), whereas only 7.5% of downregulated genes (33 of 438) have high usage of this codon (chi-squared test: χ2 = 5.29, p = 0.021) (Figure 5E). Similarly, the ThrACT codon is overrepresented in 9.0% of upregulated genes (20 of 222; hypergeometric distribution, p < 0.01) but only in 4.1% of downregulated genes (18 of 438) (chi-squared test: χ2 = 6.52, p = 0.011) (Figure 5F). Furthermore, if we consider the summed usage of the four Thr codons and two Ser codons from tRNAs possessing m3C, then no significant difference is observed between up- and downregulated genes, indicating the over-representation of ThrACC and ThrACT codons in upregulated genes is not due to the overall high usage of Thr and Ser in the proteins. These results support the idea that the four Thr codons are differentially recognized by m3C-modified tRNAThr and that m3C-modified tRNAThr predominantly interacts with ACU and ACC to enhance translation. An analysis of statistically significant use of ACC and ACT codons in genes in Gene Ontology categories reveals highly significant enrichment of translation functions in genes upregulated by MMS and in intermediary metabolism categories in genes downregulated by MMS, as shown in Figure S4.

Figure 5

Proteomic analysis reveals that exposure to the SN2 alkylating agent, MMS, causes selective translation of mRNAs possessing biased use of the threonine codons ACC and ACT that are decoded by of tRNAThr(IGU), one of three tRNAs possessing MMS-regulated m3C at position 34 (Figure 3). SILAC-based proteomic analysis revealed that MMS treatment induced up- or downregulation of 694 of 2381 proteins independent of changes in mRNA levels (Table S3). The upper box plots show that the threonine codons ACC (A) and ACT (B) are used significantly more frequently in genes for proteins upregulated by MMS exposure than in downregulated proteins, with significance demonstrated by a Student’s t test (p = 3.8 × 10–9 and p = 8.0 × 10–7, respectively). Similarly, the distribution plots in middle panels C and D demonstrate a significant MMS-induced shift in the population of proteins enriched in codons ACC (C) and ACT (D), with significance demonstrated by a Kolmogorov–Smirnov test (p = 3.7 × 10–6 and 1.9 × 10–5, respectively). The lower tables show that the threonine codons ACC (E) and ACT (F) are significantly enriched in genes upregulated by MMS exposure (29 of 222; hypergeometric distribution, p < 0.01), whereas the codons are under-represented relative to genome averages in genes downregulated by MMS (chi-squared test: ACC, χ2 = 5.29, p = 0.021; ACT, χ2 = 6.52, p = 0.011). Figure S4 shows an analysis of statistically significant use of ACC and ACT codons in MMS-regulated genes (Table S3) in Gene Ontology categories, and Table S4 lists proteins with threonine content >10%.

Discussion

We recently discovered a system of translational control of the cellular stress response, in which stress-induced reprogramming of the two-dozen modified ribonucleosides in tRNA causes the selective translation of codon-biased mRNAs for critical response proteins.[6−8,33] The studies revealed that mechanistically different chemical stressors produced signature patterns of changes in the relative quantities of the various tRNA modifications. For example, H2O2 treatment increased the levels of Cm, m5C and m22G, whereas HOCl, NaAsO2 and MMS either did not affect these modifications or lowered their levels.[7] These results raised the question of the predictive power of tRNA modification patterns to identify a specific stressor and to distinguish among chemically similar stressors. To address these questions, we developed a convergent platform of bioanalytical and mathematical tools to perform a systems-level analysis of patterns of change in the tRNA modifications in cells exposed to four oxidizing agents and five alkylating agents. This highly precise analysis (12% variance among biological replicates; see ref (7)) revealed significant class-specific features that distinguished oxidizing agents from alkylating agents (Figure 3), with 14 modifications forming the basis for a data-driven model that predicted toxicant chemistry with >80% sensitivity and specificity (Figure 4). Furthermore, tRNA modification spectra distinguished SN1 from SN2 alkylating agents, with SN2 agents causing a coordinated increase in m3C and selective translation of threonine-rich membrane proteins derived from ACT-rich genes. These results establish tRNA modifications as predictive biomarkers of exposure and illustrate a novel regulatory mechanism for translational control of cell stress response.

tRNA Modification Spectra as Biomarkers of Exposure

The class-specific patterns in the spectrum of tRNA modifications revealed by hierarchical clustering analysis (Figure 3) motivated the development of a data-driven model to assess the predictive power of the system of modified ribonucleosides. Such a model would provide insights into biological function of the system and the utility of tRNA modifications as biomarkers of exposure. The resulting model identified 14 of 25 tRNA modifications in S. cerevisiae as contributing most significantly to distinguishing the two classes of agents. As discussed shortly, several of these modified ribonucleosides have been observed to play roles in the cellular response to H2O2 and MMS exposures,[6−8] including m5C, mcm5U, and mcm5s2U, which supports the biological relevance of the model. The predictive power of the model, as revealed by confusion matrix analysis, was remarkably strong, with overall sensitivity and specificity over 80%. Examination of the individual specificity and sensitivity for oxidants and alkylating agents, however, revealed that the model had greater predictive power for the latter (>95% sensitivity and specificity versus 76–78% for oxidants; Figure 4). Although the basis for this difference is unclear, it is possible that alkylating agents produce a more complicated or extensive response in terms of activating more gene expression or requiring a broader set of responses than oxidizing agents for the same level of cytotoxicity. By any mechanism, it is striking that changes in the levels of 14 modified ribonucleosides are highly predictive of exposure chemistry in a manner similar to that with transcriptional profiling, proteomics, and metabolomics.[34−37] The predictive power of the tRNA modifications likely lies in their link to the system of biased codon use in families of stress response genes, for which we have observed that different stresses uniquely upregulate specific sets of survival and damage mitigation genes.[6−8,20,22,33] While further study is needed to compare the predictive power of tRNA modification patterns relative to other large ‘omic data sets following exposure to these toxicants, the apparent mechanistic interdependence of tRNA modifications, selective translation of survival proteins, and the resulting changes in cell phenotype suggests the potential for strongly parallel and complementary biomarker signatures.

Assigning Biological Function to Stress-Predictive Ribonucleosides

While some of the stress-regulated modified ribonucleosides in Table 1 are located at more than one position in a tRNA and in multiple tRNA species, which complicates interpretation of their biological function, 10 of them are found in the anticodon loop, which suggests possible involvement in codon recognition and interactions with aminoacyl-tRNA synthetases.[38] A mechanistic link between stress-induced tRNA reprogramming and cell stress response is illustrated by oxidation-induced increases in the level of m5C (Figure 3). Among other positions in several tRNAs, m5C is located at the wobble position of tRNALeu(CAA) that reads the UUG in mRNA, and it confers resistance to H2O2 by regulating translation of critical stress response proteins from TTG-enriched genes.[8] So, it is not surprising that m5C appears in the basis set of 14 ribonucleosides in the model. Similar arguments can be made for other ribonucleosides in the basis set (Table 1), several of which are located in a single species of tRNA (mcm5U, mcm5s2U, m3C, Am, yW, m1I). For example, alkylating agent-induced increases in the relative levels of mcm5U and mcm5s2U (Figure 3), which are part of the 14 ribonucleoside basis set for the exposure model, are consistent with our observation that modification of the wobble position of tRNAArg(UCU) with mcm5U confers resistance to MMS by promoting the translation of proteins from genes enriched with the cognate AGA codon, which represent a specific group of DNA damage-response genes.[6] So, again, it is not surprising to see elevations of mcm5U and mcm5s2U in MMS-exposed cells. While changes in the number of copies of specific tRNAs could contribute to the stress-induced changes in levels of tRNA modifications, along with increased modifying enzyme activity, our previous tRNA-seq analysis revealed that H2O2 and MMS induced nearly identical patterns of up- and downregulation for 58 tRNAs, with 18 tRNAs showing opposing changes for the stresses.[39]

Analogous arguments may provide insights into the role of other stress-regulated ribonucleosides in Table 1. For example, yW and m1I are Trm5-dependent modifications that are downregulated by exposure to the alkylating agents (Figure 3), with yW located at position 37 of tRNAPhe(GAA) and m1I located at the wobble position 34 of tRNAAla(IGC), respectively. Niu and co-workers have shown that trm5 is an essential gene involved in cell cycle regulation,[40] which raises the possibility that yW and m1I play a role in regulating the translation of cell cycle proteins during the alkylation stress response. Further evidence for mechanistic linkage between stress-regulated tRNA modifications and the stress response arises in the observation of subclass signatures for both oxidizing (Am changes with ONOO– exposure) and alkylating agents (m3C changes with SN2 alkyllators MMS and EMS), as discussed next.

tRNA Modification Patterns Distinguish Different Oxidizing Agents

While the behavior is more striking for the alkylating agents, the four oxidizing agents each showed unique patterns of change in the 23 tRNA modifications (Figure 3 and Table S2). For example, ONOO– uniquely and significantly caused an increase in Am, whereas no significant change in this ribonucleoside was observed with the other three oxidizing agent samples (Table S2). The basis for this behavior may lie in the different chemistries of reactive oxygen species (ROS) and reactive nitrogen species (RNS) in terms of molecular damage.[41,42] Am is found only at position 4 of tRNAHis(GUG) in S. cerevisiae.[10,43] Analysis of the yeast genome using the Codon Counting Database[44] reveals that the cognate codon for this tRNA, CAC, is found to be significantly enriched in 410 genes in the yeast genome (hypergeometric distribution, p value < 0.01; Table S4). Of these CAC-enriched genes, yeast flavohemoglobin is a nitric oxide oxidoreductase, which directly interacts with ONOO– and is involved in NO detoxification. It is well-known to play a role in the oxidative and nitrosative stress responses.[45]

Interestingly, among spectra of the other three oxidizing agents, H2O2 exposure is more similar to γ-rad than to another peroxide, TBHP. This may be due to the fact that H2O2 and γ-rad generate hydroxyl radicals that lead to chemically distinct types of molecular damage than the peroxyl radicals arising from TBHP. Thus, there is likely to be a discrete biological mechanism underlying the tRNA modification patterns for different chemical stresses, with the subclass signatures for SN1 and SN2 alkylating agents providing a striking illustration of the phenomenon.

Correlations between Subclass Signatures and Mechanisms of Action among Alkylating Agents

The hierarchical clustering analysis in Figure 3 reveals a significant differential alteration of modified ribonucleosides by the two chemically distinct types of alkylating agents. MMS and EMS represent one class of alkylating agent that appears to react with nucleophilic sites in nucleic acids, proteins, and other cellular molecules[46,47] by a bimolecular nucleophilic substitution reaction (SN2), whereas IMS, MNNG and NMU appear to react by a unimolecular mechanism (SN1).[48] The relevance of this mechanistic classification for genetic toxicology has been challenged[49] but the two groups of alkylating agents produce distinctly different types of molecular damage,[48] and we use the SN1/SN2 nomenclature for the sake of clarity. The analyses in Figure 3 and Table S2 clearly and significantly distinguish MMS and EMS from the other three agents in terms of tRNA modification spectra. For example, the relative levels of m7G and m3C in MMS- and EMS-exposed cells increased more significantly than in IMS-, MNNG-, and NMU-exposed cells. In contrast, the levels of Am and Um did not change in response to MMS and EMS, but they were elevated when cells were treated with the other three alkylating agents. These results suggest that the cells are responding to the damage caused by two classes of alkylating agents with different translational responses that reflect activation of different cellular survival pathways. This point was demonstrated in 1983 by Smith and Grisham[50] in studies of cytotoxicity in yeast. They observed that the toxicity of SN2 alkylating agents, such as MMS, arises mainly by damage to membrane proteins, whereas the toxicity of SN1 alkylating agents, such as MNNG, results from interference with DNA replication.[50] Our data provide a mechanistic model for these observations. The SN2 agents, MMS and EMS, increased the level of m3C, which is known to occur at position 32 of tRNAThr(IGU), tRNASer(UGA), and tRNASer(CGA) and in the variable stem of a Ser tRNA.[51−53] The link between these tRNAs and the proposed membrane protein target of SN2 agents may lie in the fact that Thr and Ser are significantly enriched in membrane proteins in eukaryotes.[54−56] Analysis of the S. cerevisiae genome reveals that, on average, 5.9% of each yeast protein is composed of Thr, whereas only 83 of ∼6000 genes in the yeast genome code for proteins composed of >10% Thr[29] (Tables S4 and S5). Functional analysis of the Thr-enriched proteins by annotation with Gene Ontology functional attributes using the DAVID bioinformatics resources[29] revealed that 61 of these 83 genes code for cell wall and cell membrane proteins (Tables S5 and S6). While membrane proteins were not well represented in our proteomic analysis, which is to be expected given the predominance of soluble proteins in our isolation method, there was a significant upregulation of proteins from genes enriched with ThrACT and the optimal codon ThrACC in cells exposed to MMS (Figure 5). Along with the observation that loss of the ability synthesize m3C confers sensitivity to EMS and MMS in Δtrm140 cells (Figure S2), these results are consistent with the idea that the EMS- and MMS-induced increase in m3C reflects the need for more efficient translation of Thr-enriched membrane proteins in the response of S. cerevisiae to exposure to the SN2 alkylating agents. There are certainly factors other than reprogramming of m3C-cointaining tRNAs and selective translation of ACT- and ACC-enriched genes in the translational control of the MMS stress response, as suggested by the downregulation of some proteins from ACT- and ACC-enriched genes (Figure 5), such as contributions from tRNA modifications and biased codon usage, as well as factors regulating tRNA charging, ribosome loading, and translational elongation rates.

Conclusions

We employed novel bioanalytical and computational tools to demonstrate that cells respond to exposure to different classes of chemical stressors by reprogramming a system of modified ribonucleosides in tRNA, with unique patterns distinguishing the class of chemicals, as well as the subclass. Multivariate statistical analysis and data-driven modeling of the tRNA modification patterns proved to be highly predictive of toxicant chemistry, with bioinformatic and proteomic analyses establishing a mechanistic link among exposure, tRNA modifications, and codon-biased translation of response proteins. With implications for other organisms, including humans, these results establish tRNA modifications as predictive biomarkers of exposure and illustrate a novel regulatory mechanism for translational control of cell stress response.