Global Epitranscriptomics Profiling of RNA Post-Transcriptional Modifications as an Effective Tool for Investigating the Epitranscriptomics of Stress Response*

The simultaneous detection of all the post-transcriptional modifications (PTMs) that decorate cellular RNA can provide comprehensive information on the effects of changing environmental conditions on the entire epitranscriptome. To capture this type of information, we performed the analysis of ribonucleotide mixtures produced by hydrolysis of total RNA extracts from S. cerevisiae that was grown under hyperosmotic and heat shock conditions. Their global PTM profiles clearly indicated that the cellular responses to these types of stresses involved profound changes in the production of specific PTMs. The observed changes involved not only up-/down-regulation of typical PTMs, but also the outright induction of new ones that were absent under normal conditions, or the elimination of others that were normally present. Pointing toward the broad involvement of different classes of RNAs, many of the newly observed PTMs differed from those engaged in the known tRNA-based mechanism of translational recoding, which is induced by oxidative stress. Some of the expression effects were stress-specific, whereas others were not, thus suggesting that RNA PTMs may perform multifaceted activities in stress response, which are subjected to distinctive regulatory pathways. To explore their signaling networks, we implemented a strategy based on the systematic deletion of genes that connect established response genes with PTM biogenetic enzymes in a putative interactomic map. The results clearly identified PTMs that were under direct HOG control, a well-known protein kinase pathway involved in stress response in eukaryotes. Activation of this signaling pathway has been shown to result in the stabilization of numerous mRNAs and the induction of selected lncRNAs involved in chromatin remodeling. The fact that PTMs are capable of altering the activity of the parent RNAs suggest their possible participation in feedback mechanisms aimed at modulating the regulatory functions of such RNAs. This tantalizing hypothesis will be the object of future studies.

The viability of a cellular system is largely predicated on its ability to recognize and respond to a variety of environmental conditions (1,2). External stimuli elicit counteracting compensatory responses intended to maintain growth and ensure survival. Cell adaptation in response to stressors, such as sudden variations of temperature, osmotic conditions, and availability of nutrients, produces significant variations of enzyme activity and metabolism, which are orchestrated by well calibrated regulatory controls (3). Owing to its multifaceted activities in protein synthesis and newly discovered functions in gene regulation, RNA is uniquely positioned among cellular components to act as a communication node between metabolic and regulatory mechanisms of stress response (4,5). This principle is realized in the activity of mitogen-activated protein kinase (MAPK) 1 pathways induced in many eukaryotic systems by heat and osmotic stress (6,7). In S. cerevisiae, activation of the hyperosmolarity glycerol (HOG) MAPK pathway produces conspicuous variations of the activity of ion channels, glycerol export, general protein machinery, and cell cycle progression (8 -11). The transcriptional response to osmostress includes the induction and stabilization of hundreds of mRNAs (12)(13)(14) and the up-regulation of specific families of long noncoding RNAs (lncRNAs) (15), as well as a general down-regulation of gene expression (15)(16)(17)(18). Significantly, the latter class of RNA has been shown to induce chromatin remodeling and variations of nucleosome occupancy that have lasting effects on cellular memory and epigenetic makeup (21)(22)(23).
In recent years, covalent modification and nucleotide editing have increasingly been recognized as effective mechanisms for modulating the functions of RNA (24). Natural RNA contains over 100 variants of the four canonical bases, which are cataloged in the RNA Modification (25,26) and MODOMICS (27,28) databases. These post-transcriptional modifications (PTMs) are generated by specific biogenetic enzymes and possess the ability to stabilize single base pairs and alternative hydrogen bonding patterns, which contribute to the astonishing diversity of RNA structure (29,30). At the same time, the observation that unmodified tRNAs were not loaded efficiently with the respective aminoacyl groups, whereas specific modifications prevented incorrect loading altogether (31), revealed that PTMs play essential roles in protein-RNA recognition. Analogous effects were observed in RNA-RNA recognition for hypermodified nucleotides in the tRNA anticodon loop, which were shown to induce translational recoding by affecting the accuracy of codon-anticodon interactions (24). The detection of an increased incidence of anticodon loop modification upon treatment with hydrogen peroxide has led to the assertion that translational recoding is an essential mechanism of response to oxidative stress (32,33). Owing to the almost exclusive emphasis on tRNA and rRNA as abundant sources of PTMs, little else is known about their biological functions and distribution in other classes of RNAs. However, their enzymatic biogenesis supports the tantalizing hypothesis that PTMs may constitute key elements of common feedback mechanisms for modulating the activity of regulatory RNAs.
To support the investigation of entire epitranscriptomes, we have developed a novel approach for the analysis of all the PTMs contained in a target cell (34). This approach relies on the extraction of total RNA from the sample, followed by digestion into individual mononucleotides and analysis by direct infusion electrospray ionization (ESI) (35). The selected workflow eschews high-resolution chromatography to minimize possible separation bias, reduce losses during sample transfer/manipulation, and expedite the analysis by reducing the number of necessary operations. Alternative combinations of high-resolution and multistep tandem mass spectrometry (MS n ) (36,37), or ion mobility spectrometry-mass spectrometry (IMS-MS) (38 -40) and gas-phase fragmentation techniques can be independently employed to complete the characterization of the mononucleotide mixtures with excellent sensitivity (34,41). Either combination can afford both identity and abundance of any PTM present in the total RNA extract, thus providing comprehensive and unbiased profiles at the entire transcriptome level. This approach has been shown capable of revealing not only significant differences, but also subtle variations of PTM expression between samples grown in different media, or obtained altogether from different organisms (34).
In this report, we employed our global profiling approach to investigate the possible existence of a broad, general link between PTM expression and stress response, which goes beyond the known mechanism of tRNA-based translational recoding. Preliminary experiments were aimed at evaluating whether different types of stresses, such as heat shock and hyperosmotic stress, produced distinctive PTM profiles that could be ascribed to different response mechanisms. We subsequently focused on the HOG pathway because of the wealth of available information on its role in the hyperosmotic stress response (9,10), including the well-characterized induction and stabilization of mRNAs and lncRNAs (12)(13)(14)(15). A key component of this pathway is the highly conserved MAPK protein Hog1, a stress-activated protein kinase homologous to mammalian p38 and c-Jun N-terminal kinase, which is involved in multiple phases of the cell cycle and has the ability to reprogram gene expression in the presence of stress (42). For this reason, we first compared PTM profiles obtained from wildtype S. cerevisiae (i.e. WT) in the absence and presence of osmostress to identify possible PTM variations that correlated with the known induction of stress-response coding and noncoding genes (43). Then, we performed the profiling of a mutant strain in which the HOG1 gene was deleted (i.e. hog1⌬) to discriminate hog1-dependent PTMs from those that were induced by osmostress without falling under direct control of HOG1. Finally, we implemented an interactomewalk strategy in which the genes connecting HOG1 with PTM biogenetic enzymes in putative interactomics maps were sequentially deleted to retrace the signaling cascade and show a direct relationship between PTM expression and this stress response pathway. For the purpose of obtaining proof-ofprinciple and evaluating the significance of this approach, we explored the putative pathway leading from HOG1 to RIT1, the enzyme responsible for one of the most prominent stressinduced modifications observed in the study. The results clearly showed the merits of monitoring the broader impact of individual gene deletions on PTM biogenesis at the entire transcriptome level.

MATERIALS AND METHODS
Preparation of Cellular Extracts-Saccharomyces cerevisiae strain BY4741 (defined here as wildtype, WT), as well as hog1⌬::kanMX, ssb1⌬::kanMXi, rit1⌬::kanMX, and hog1⌬::kanMX derivatives obtained according to established procedures (44,45), were purchased directly from Open Biosystems (Huntsville, AL). Yeast strains were grown in yeast extract, peptone, dextrose (YPD) medium. Each strain was streaked onto YPD agar and incubated at 30°C. For all studies, an individual colony was selected from a plate and placed into an individual tube containing 20 ml of YPD. The culture was incubated at either 30 or 37°C with 200 rpm gyration. Optical density at 600 nm (OD 600 ) was monitored on a ThermoFisher Scientific (Waltham, MA) Nanodrop 2000c spectrophotometer until a value greater than 0.3 units was achieved. Liquid cultures left untreated were diluted to a final 0.3 OD 600 and centrifuged at 6000 ϫ g for 5 min to obtain pellets that contained approximately the same number of cells. Each sample was grown in parallel in at least five separate cultures (biological replicates) and the material obtained from each of them was submitted to five independent analyses (technical replicates). To assess the possible effects of osmostress on RNA modification profiles, S. cerevisiae was grown to mid-log phase and induced with 0.4 M NaCl. Liquid cultures were then incubated for an additional 15 min, as these conditions have been found to specifically overexpress ϳ343 coding genes and ϳ173 lncRNAs (43). In all cases, the reported results were the average of at least five biological replicates, each of which submitted to five technical repeats.
Sample Preparation and Mass Spectrometric Analysis-Each pellet was disrupted by using Denaturation Solution (Life Technologies, Grand Island, NY). Total RNA was extracted by using the ToTALLY RNA Extraction Kit (Life Technologies), which is based on a typical phenol/chloroform procedure. The RNA was precipitated by using cold isopropanol and then treated with DNase 1 (New England Biolabs, Ipswich, MA) in 1X DNase buffer to remove any remaining DNA. The recovered RNA was subsequently desalted by ethanol precipitation overnight and reconstituted in 50 l of RNase-free water (Sigma-Aldrich, St. Louis, MO). The concentration of intact total RNA from each sample was measured by UV absorbance at 260 nm. Nuclease P1 and phosphodiesterase 1 from snake venom (Sigma-Aldrich) were employed to complete the digestion of RNA into individual mononucleotides, as previously described (46). Immediately before analysis, final samples were diluted 1:10 in 150 mM ammonium acetate and 10% isopropanol.
Mass spectrometric analyses were carried out by direct infusion ESI on either a Waters (Milford, MA) Synapt G2 HDMS IMS mass spectrometer, or a ThermoFisher Scientific LTQ-Orbitrap Velos mass spectrometer. All analyses were performed in nanoflow ESI mode by using quartz emitters produced in house. Up to 5 l samples were typically loaded into each emitter by using a gel-loader pipette tip. A stainless steel wire was inserted in the back-end of the emitter to supply an ionizing voltage that ranged between 0.9 and 1.2 kV. Source temperature and desolvation conditions were adjusted by closely monitoring the incidence of ammonium adducts and water clusters (47).
Analyses performed on the IMS-MS instrument provided heat maps that offered comprehensive representations of the entire complement of PTMs present in the sample. For these experiments, the instrument was calibrated by using a 2 mg/ml solution of cesium iodide in 50:50 water/methanol. This external calibration provided typical ϳ9 ppm accuracy (34,41). For comprehensive mixture analysis, the Tri-WAVE region was held at a pressure of ϳ4.40 mbar (uncalibrated gauge reading) by a 90 ml/min flow of N 2 and 180 ml/min of He. It was operated with an ϳ650 m/s IMS wave velocity, a 40 V wave height, a 109 m/s transfer wave velocity, and a 2.0 V transfer wave height. Time aligned parallel (TAP) dissociation and mass-selected time-aligned (MaSTeR) fragmentation of isobaric/isomeric species were performed as previously described (34,41). Briefly, TAP experiments involved transmitting all ions produced by the nanospray source through the quadrupole (Q) element of the instrument, and then injecting them into the Tri-WAVE region to enable their discrimination according to their arrival time (t D ). Activation was performed in parallel for all ions that reached the transfer region of the instrument located immediately after the Tri-WAVE element. The fragments were readily detected in the time-of-flight (TOF) analyzer, thus preserving the relationship between each set of fragments and the t D of the respective precursor(s) (41). In contrast, MaSTeR experiments involved transmitting only a specific population of ions sharing the same m/z value through the mass-selective Q. The mass-selected precursors, which could comprise multiple isobaric/ isomeric species, were then allowed to disperse on the time domain in the Tri-WAVE element and then activated in the transfer region (34). Therefore, the difference between these techniques is that the former can provide fragmentation data simultaneously for all possible ions contained in the sample, whereas the latter mimics a true tandem MS experiment in which precursor ions are sequentially selected according to both m/z and t D characteristics in the Q and Tri-WAVE regions, respectively.
Analyses performed on the LTQ-Orbitrap provided accurate mass and fragmentation data that were used to confirm the initial assignments afforded by the IMS-MS experiments. For high resolution/ accuracy determinations, the instrument was calibrated by using an anion mixture that contained sodium dodecyl-sulfate, sodium taurocholate, and Ultramark. We have previously shown that this external calibration can afford typical ϳ200 ppb accuracy (34,41). For analytes the size of ribonucleotides, this level of accuracy reduces considerably the number of possible elemental compositions that may fit the observed masses, thus increasing the level of confidence in the data interpretation. Tandem mass spectrometry (MS/MS) (48), multistep activation experiments (MS n ) (36,37), and consecutive reaction monitoring (CRM) (49) were performed as previously described (34,41). MS n and CRM were specifically employed to corroborate the identity of isobaric/isomeric species, which was initially obtained from IMS-MS determinations.
Data Treatment-The experimental masses obtained from either the IMS-MS or LTQ-Orbitrap instrument were transformed into a readable file format by using MSConvert (50), and then employed to search the RNA Modifications (51), MODOMICS (52), or a home-built database (34). The initial database "hits" were subsequently confirmed by TAP and MaSTeR determinations on the IMS-MS instrument. In the case of isomeric/isobaric species, the results were further corroborated by MS n and CRM analyses on the LTQ-Orbitrap. At the same time, the initial heat map data provided by IMS-MS determinations were also subjected to data reduction, alignment, and scaling to enable point-by-point data subtraction. The results were visualized in heat map format by using OriginPro 9.1 (Origin Lab, North Hampton, MA) to readily highlight subtle variations of PTM expression between samples.
The abundance of detected species was expressed in relative, rather than absolute terms, to enable an immediate appreciation of over-/under-expression under the selected experimental conditions. As discussed in detail in reference 34, the abundance of each species was expressed as percentage of the cumulative abundances of the canonical ribonucleotides detected in the same sample, defined as abundance versus proxy (AvP): in which ai x is the signal intensity of the analyte x, whereas cr i is the intensity of each of the four canonical ribonucleotides. Based on the observation that the canonical ribonucleotides are under genetic rather than metabolic control, their overall abundance is only minimally affected by the presence of modifications. For this reason, the canonical ribonucleotides can be pooled together to serve as a proxy internal standard for obtaining relative abundances expressed in percentage terms. In previous work, the standard addition method was employed to determine the absolute amounts of selected PTMs in yeast, which allowed the calculation of the corresponding AvP values to confirm the validity of this approach (34).
All AvP values included in this report were the average of the results acquired from a minimum of five different samples grown separately under the same experimental conditions (biological replicates). The data obtained from each strain/condition were compared with those afforded by the wildtype BY4741 strain grown in YPD at 30°C, as described above, which represented the baseline for these determinations. A paired two-tailed Student's t test was employed to assess the statistical significance of the comparison, which was accomplished at the 95% confidence level (i.e. corresponding to a p value of 0.05). The test was applied to each species detected in both the sample of interest and wildtype control. In this way, any variation with a p value smaller than 0.05 was considered significant and was marked with different shades of color in the tables according to whether the variation represented over-or underexpression in comparison to the wildtype.

RESULTS AND DISCUSSION
Global Epitranscriptomics Analysis-The cellular response staged by an organism under stress involves the concurrent and coordinated activation of numerous metabolic pathways. Therefore, investigating the fate of an individual PTM can provide only a limited view of its role in the overall response mechanism. In contrast, the ability to monitor multiple PTMs at the same time can provide valuable insights into their synergistic activities and coordination at the system level. This type of information can be obtained through analytical strategies capable of multiplexing PTM detection and accurately assessing their mutual expression levels. In this direction, we have recently developed a strategy for obtaining complete PTM profiles, which involves the analysis of total RNA extracts obtained from cell lysates (34). The material is subsequently digested by selected exonucleases to produce mononucleotide mixtures that are analyzed without front-end separation procedures (see Materials and Methods). The complexity of these types of mixtures can be readily resolved by high-resolution MS analysis, followed by tandem MS for structure verification (34). Alternatively, IMS-MS analysis can accomplish the task by discriminating analytes according to both mass and ion mobility behavior. The latter is determined by the probability of ions to interact with background gas while traveling in a moderate electric field, which is a function of size and conformation (38 -40). The interactions affect the time of arrival (t D ) of each ion, which is plotted against the corresponding mass-to-charge ratio (m/z) to obtain comprehensive 3D-plots or heat maps. Fig. 1 shows representative heat maps in which the relative intensity of detected ionscorresponding to the third dimension of a 3D-plot-is represented by a color gradient. In this way, heat maps can provide comprehensive visual representations of all species in the samples, which enabled immediate comparisons between strains grown under different environmental conditions (e.g.  ) or B, 37°C (WT 37 , heat shock); and C, rit1⌬ strain grown at 30°C (rit1⌬). Panel D, shows the expanded regions containing the Ar(p) modification, which enable one to readily compare its expression levels in the respective samples. The intensity axis of the heat maps were scaled to the same level to enable direct comparisons by using the color gradient provided on the right. The expanded region in panel A shows the separation achieved on the time domain for the U/⌿ isobars. Fig. 1A and 1B), or between wildtype and mutant strains (e.g. Fig. 1A and 1C, vide infra).
The data depicted in Fig. 1A were obtained from S. cerevisiae strain BY4741, which was grown under normal conditions in YPD medium at 30°C (WT control, see Materials and Methods). The mass information provided on the y axis was readily employed to search a nonredundant database compiled in house from the entries available in the RNA Modifications (51) and MODOMICS (52) databases. The initial database "hits" included numerous isomeric/isobaric species that shared the same elemental composition and, thus, could not be readily discriminated on the basis of mass alone. However, as shown earlier for the uridine/pseudouridine couple (41) and a series of methyl-guanine products (34), the different structures associated with the isomeric species translated into distinctive ion mobility behaviors that were readily differentiated on the t D dimension (see for example Fig. 1A inset). Efforts to appropriately calibrate the t D axis in different hardware platforms are currently underway to enable the use of characteristic t D values as unique identifiers for individual PTMs. In the meantime, we capitalized on the ability to disperse isomeric/isobaric species on the time domain to enable their separate gas-phase activation and fragmentation. In what has been defined as mass-selected time-resolved (MaSTeR) dissociation determinations (34), each target precursor can be isolated in the mass-selective quadrupole of the instrument, allowed to disperse temporally in the ion mobility element, and then activated in the transfer region, before detection of the ensuing fragments in the time-of-flight analyzer (see Materials and Methods).
This type of analysis can be exemplified here for the species detected at m/z 365.024 in Fig. 1A, which was provisionally assigned to the isomeric/isobaric couple ac 4 C/f 5 Cm by the initial database search. The arrival time distribution (ATD) obtained by isolating the precursor ion in the mass-selective quadrupole is shown in Fig. 2A. Analyzed with the aid of an appropriate deconvolution algorithm (see Materials and Methods), the recorded signal revealed the presence of two discrete components that were discriminated on the t D dimension despite their widely different abundances. The identities of these components were confirmed by examining pertinent reconstructed ion chromatograms (RICs) obtained by extracting the masses of unique isomer-specific fragments from the recorded fragmentation dataset. Indeed, the t D of the signal present on the RIC of m/z 154 (Fig. 2B) matched that of the more abundant species in the overall ATD ( Fig. 2A), thus identifying this component as the ac 4 C isobar. In contrast, the signal on the RIC of m/z 136 (Fig. 2C) matched the less abundant species, which confirmed its identity as the f 5 Cm isobar. As corroborated separately by regular tandem MS, these isobar-specific fragments were produced by the loss of the corresponding nucleobase moieties that bore the covalent modifications (see supplemental Fig. S4). A manuscript in preparation will provide the unique precursor 3 fragment transitions necessary to confirm the identity of all PTMs observed in these yeast samples. This type of analysis was carried out for all initial database "hits" to complete their positive identification (see Supporting Material for additional representative data). The PTMs corroborated by fragmentation data were listed in Table IA.
The quantitative determination of PTMs has traditionally suffered from the dearth of suitable standards for the Ͼ100 known ribonucleotide variants. We have shown that purified, exogenous tRNA samples could be effectively used as controlled sources of PTM standards in a modified version of the standard-additions method (34). According to this strategy, accurately known amounts of exogenous tRNA were added to the total RNA extract, so that subsequent exonuclease treatment could release directly in situ the desired PTM standards present in its molecule. An alternative strategy was proposed for applications aimed at determining mutual relative variations of expression levels, rather than absolute PTM concentrations. According to this, the entire complement of canonical bases contained in the cell, which is under genetic rather than metabolic control, was employed as a proxy internal standard to express the relative abundances of detected PTMs (see Materials and Methods) (34). The abundances versus proxy (AvPs) produced by this treatment were compared with corresponding values calculated from absolute concentrations obtained by the standard-addition method. The excellent match between these

cerevisiae grown at A) 30°C (WT, control); B) 37°C (WT 37 , heat shock); C) 30°C after treatment with 0.4 M NaCl (WT S , osmostress); D) 37°C after treatment with 0.4 M NaCl (WT S,37 ); hog1⌬ grown at E) 30°C (hog1⌬); F) 37°C (hog1⌬ 37 ); G) 30°C after treatment with 0.4 M NaCl (hog1⌬ S ); and H) 37°C after treatment with 0.4 M NaCl (see Materials and
Methods) values provided the necessary validation for this quantification strategy (34). The IMS-MS data depicted in Fig. 1A were processed according to this treatment to assess the relative expression levels of all PTMs in the sample (Table IA), regardless of the availability of individual standards (see Materials and Methods).

Effects of Heat Shock and Osmostress on PTM Expression-
The experimental strategies summarized above were employed to obtain the global profiles of samples that had been subjected to different types of stresses. For example, the heat map in Fig. 1B was obtained from a sample of wildtype S. cerevisiae BY4741 that was grown in YPD medium at 37°C (WT 37 , see Materials and Methods), under conditions known to induce a typical heat shock response (53). Visual comparison between these data and those provided by WT control (Fig. 1A) revealed remarkable differences, which were confirmed by completing the identification of the PTMs in the sample (Table IB). The baseline represented by the profile of WT control consisted of a total of 41 distinct PTMs (Table IA). In contrast, the WT 37 sample displayed 50 PTMs, out of which 32 were present also in WT control and the remaining 18 were newly expressed PTMs unique for this type of sample (marked in red in Table IB). The observed increase of PTM content was consistent with previously published reports that heat shock activation caused broad overexpression of RNA components within the cell (54,55). Among the 32 PTMs in common between the two types of samples, four were detected with comparable relative abundances in both, whereas 14 and 14 were, respectively, under-and over-expressed in WT 37 (marked respectively with darker and lighter shades of blue in Table IB). It must be emphasized that at least five biological replicates were analyzed for each type of sample to assess the uncertainty intrinsic in the quantitative determination (see Materials and Methods). Only PTMs with relative abundances that deviated from those of the WT control with a p value Ͻ 0.05 were deemed significant and marked as being under-/ overexpressed. Smaller deviations with p Ͼ 0.05 matched the typical biological reproducibility obtained from yeast samples submitted to this type of analysis and, thus, were considered as normal experimental fluctuations (see Materials and Methods).
In analogous fashion, S. cerevisiae strain BY4741 was grown in YPD medium at 30°C and then exposed for 15 min. to an environment containing 0.4 M NaCl (WT S , see Materials and Methods), which induces a characteristic osmostress response (15). Also in this case, the respective heat map (not shown) provided an outcome that was significantly different from that of WT control and, remarkably, from that of the WT 37 sample, as well. The PTM profile revealed by completing positive identification (Table IC) showed that WT S contained a total of 38 PTMs, with eight newly expressed (red) and 30 in common with WT control (blue), out of which 7 and 22 were, respectively, under-and overexpressed (shades of blue). In addition, seven of the eight newly expressed PTMs (absent in WT control) were also present in WT 37 , whereas the remaining one was unique for this sample.
These experiments clearly indicated that the cellular responses to these different types of stresses involved profound changes in the production of specific RNA PTMs. The changes were not limited to variations of the expression levels of common PTMs, but comprised also the activation of new pathways that were previously dormant in the absence of stressors. According to the information included in the RNA Modification (51) and MODOMICS (52) databases, 11 of the 18 new PTMs activated by heat shock and three of eight associated with osmostress have been previously described in eukaryotic, but not necessarily yeast tRNA. Representative tandem MS data obtained for this set of PTMs is included in the Supporting Material section to corroborate their detection in this sample. In this regard, it is important to note that the information available on the distribution of PTMs in different organisms and classes of RNAs is still rather fragmentary. In most cases, the provenance of a certain PTM cited only the initial study that reported its first detection/characterization, rather than the results of comprehensive surveys. The vast majority of studies have traditionally focused on tRNA and rRNA, not only because of their abundant PTM content, but also for the availability of convenient isolation protocols. Further, it is not clear whether any distribution information present in the common databases was obtained from samples harvested under stress conditions. Therefore, the provenance information included in Table I should be considered with discretion. These considerations aside, it was very significant that none of the new PTMs detected under stress corresponded to those implicated in the mechanism of oxidative stress response involving tRNA-based translational recoding (32,33). This observation suggested that a sizeable portion of the newly activated pathways were likely to rely on very different mechanisms.
A closer examination of the new PTMs expressed in WT S revealed that seven were in common with the WT 37 sample, but 1 and 11 were unique for either the former or the latter, respectively, thus pointing to the existence of distinct response pathways involving PTM biogenesis. It is well established that thermal stress activates a response based on the prominent expression of heat shock proteins (HSP), which exert their protective effects by limiting protein denaturation (30,56). However, it has been shown that the HOG pathway can be also activated at higher temperatures (57), thus raising the possibility that the respective signaling pathways may be closely intertwined. If proven correct, this hypothesis would provide an excellent explanation for the seven common PTMs that were induced by both heat shock and osmostress in our experiments.
In addition to the induction of new PTMs absent in WT control, stress conditions affected the expression levels of existing PTMs in different ways. Indeed, the relative abundance of 14 and 22 existing PTMs were significantly increased in the WT 37 and WT S samples, respectively, over the WT control (i.e. displayed p Ͻ 0.05). In contrast, 14 and 7 existing PTMs were significantly decreased in the same samples, whereas four and one ceased to be detectable altogether under the standard analytical conditions. Also in this case, the various PTMs could be readily categorized in distinctive groups that were either stress-specific or not. Of the 20 downregulated modifications in WT 37 , 14 displayed the same behavior in WT s . At the same time, of the 32 species up-regulated/unique in WT 37 , 20 displayed the same behavior in WT s . In some cases, however, the variation directions did not match in both samples, as shown for t 6 A that was up-regulated under osmostress as compared with the control, but disappeared under heat shock, or the methyl-A isomers that had the exact opposite fate (Table I). The disappearance of specific PTMs under stress conditions was particularly intriguing and suggested the possibility that they may not be essential to the cellular response and that their functions may be effectively taken over by other stress-specific PTMs. The picture emerging from these initial experiments clearly ruled out the notion that PTM biogenesis might be controlled by a single regulatory pathway, but suggested instead the presence of multiple individual controls.

Identification of HOG-Dependent Expression Pathways-To
unambiguously identify the stress-induced PTMs that were under HOG control, we analyzed a hog1⌬ strain, which lacks the stress-activated protein kinase Hog1, the master regulator of gene expression reprogramming in response to osmostress. The hog1⌬ mutant was first grown at 30°C to identify its baseline profile in the absence of stress, which comprised the induction of four new PTMs absent in WT control (Table  IE). The mutant strain was subsequently grown at 37°C (i.e. hog1⌬ 37 ) and treated with 0.4 M NaCl at both 30°C (hog1⌬ S ) and 37°C (i.e. hog1⌬ 37,S ). The samples displayed complex profiles with up/down variations of the majority of the PTMs observed in WT under corresponding environmental conditions, as well as production of distinctive sets of new ones absent in WT (marked respectively with different shades of blue, or solid red in Table IF-IH).
The newly activated PTMs observed exclusively under stress conditions were examined according to a process of elimination to identify their putative regulatory relationships (summarized in Table II). In particular, the mutation-induced elimination of any of the newly expressed PTMs (i.e. absent in WT control) was interpreted as evidence of HOG control. The reasoning behind the elimination process could be illustrated  1 Full names available at http://mods.rna.albany.edu/. 2 Putative classes of RNAs include: ribosomal RNA (r); transfer RNA (t); unknown class (?). 3 WT, hog1⌬, ssb1⌬, and rit1⌬ strains were grown at 30°C under normal conditions (no subscript); 37°C (37 subscript); 30°C after treatment with 0.4 M NaCl (S subscript); and 37°C after treatment with 0.4 M NaCl (S,37 subscript) (see Materials and Methods). 4 Red color marks newly expressed PTMs absent in the WT sample. * Modifications previously unreported for this microorganism.
here for the newly expressed ho 5 U modification. This species was absent in WT control, but present in both WT 37 , WT S , as well as the double-stressed WT 37,S sample (see Materials and Methods), thus cementing its status as a stress-induced PTM. At the same time, this species was prominently absent in all of the hog1⌬ samples regardless of conditions, which indicated beyond doubt that its expression was linked to the HOG pathway. In contrast, the i 6 A and ncm 5 Um modifications observed in both WT 37 and WT S were not eliminated by the mutation, thus suggesting that their stress-induced biogenesis was independent from the HOG pathway.
It is interesting to note that although the production of these PTMs was triggered by both heat shock and osmostress, others displayed prominent specificity for either type of stress. For example, X was detected only in the WT S sample (and, consistent with this finding, in the double-stressed WT 37,S ) to earn the status of osmostress-induced and hog1dependent. Its reproducibility and abundance in the analyzed samples dissipated any doubt about its significance. At the same time, 10 putative PTMs were found to be heat shockinduced and hog1-dependent, but only one was heat shockinduced and hog1-independent (Table II). Finally, f 5 U was observed only when the hog1⌬ mutant was subjected to higher temperature, thus suggesting the possible activation of an alternative mechanism of heat response, which may compensate for the elimination of the Hog1 kinase.
A Strategy for Integrating PTM Biogenesis within the Stress-Response Regulatory Network-The results provided by these experiments clearly showed that the effects of heat and osmotic stress on PTM biogenesis were not mediated by a single regulatory pathway, but resulted instead from the coordinated activities of a more complex regulatory network. In S. cerevisiae (budding yeast), extensive investigation of gene expression under stress has led to the identification of genes that are consistently up-or down-regulated in response to specific stimuli (58). These genes have been defined as specific environmental stress response genes (SESR, Scheme 1) (17). In S. pombe (fission yeast), instead, common sets of genes defined as core environmental stress response genes (CESR) are involved in nearly all types of stresses (17). In this case, a central hub integrates the stimuli provided by various stressors and coordinates a general cellular response. The detection of PTMs that were activated by both heat and osmotic stress suggests that a similar organization may be present in the stress-induced PTM biogenesis observed in S. cerevisiae.
To test this possibility, we devised a strategy that used gene deletion strains to retrace putative regulatory pathways. This strategy required the prerequisite formulation of valid hypotheses on the network's composition and architecture to guide the selection of the genes to be eliminated. For example, Scheme 2A depicts the organization of genes involved in the typical cellular response to heat shock and osmostress in S. cerevisiae, as deduced from a wealth of available information (20, 42, 59 -63). In contrast, Scheme 2B was drawn by using the Cytoscape server to examine databases of genetic and protein-protein interactions, as well as transcriptional and post-transcriptional regulation networks (64 -66). This specific subsection of the extensive stress-response interactome delineates a putative network connecting essential components of the established stress-response pathways, such as HOG and various heat shock genes, to the biogenetic enzymes responsible for the formation of target PTMs, which are listed in the MODOMICS database (52). It is interesting to note that, although the organization of the established response genes (Scheme 2A) fully conformed to the prevailing view of parallel, largely independent pathways for different stimuli, the downstream region covering the putative regulation of PTM expression (Scheme 2B and 2C) assumed a much more complex architecture characterized by numerous crossovers. The interactomic map clearly suggested that genes such as HOG1 and FUS3 could provide high-level integration between the heat shock and osmostress signaling pathways, whereas SSA1, SSB1, SSB2, SSE1, and HSP82 could provide additional low-level integration and fine tuning of PTM production. This set of genes codes for proteins with broad chaperone activities, which are known to prevent aggregation, promote folding of proteins to their native states, and solubilize and refold aggregated proteins (55,67). Often, these heat shock chaperones work closely with cochaperones that are involved in signal transduction, cell cycle regulation, cellular differentiation, and cell death. Under stress, eukaryotic systems rely heavily upon Hsp70 chaperones such as Hsp82 to stabilize death-inducing mutations, such as those propagated by oncogenes (55,67).
Based on this map, we have initiated the systematic investigation of selected gene deletion strains to substantiate the putative signaling pathways responsible for PTM regulation. The PTMs of most interest were the 11 induced in either heat shock or osmostress, which were found to be hog1-dependent (Table II). The results obtained from the ssb1⌬ mutant showed that 11 of the 19 PTMs found uniquely in stressed states were ssb1-dependent, indicating that this pathway may SCHEME 1. Regulation of stress genes in fission and budding yeast. In the former, different stresses share a common core of environmental stress response genes (CESR), whereas stress-specific environmental response genes (SESR) predominate in the latter (17,58).
indeed be acting as an integration hub responsible for the downstream initiation of multiple enzymatic processes involved in PTM biogenesis. The fact that the majority of these PTMs are controlled by genes located downstream in these established stress-response pathways appears to rule out their possible participation in sensing mechanisms, which are typically located upstream. Additionally, the two PTMs that appeared in WT, but not in any specific stress states, were also found to be ssb1-dependent (Table II). Further, a total of eight PTMs were found to be expressed even in the absence of SSB1 indicating that SSB1 cannot be the only hub responsible for signaling PTM biogenesis. The next step in retracing the proposed map consisted of analyzing strains in which the specific biogenetic enzymes were deleted. The case of the RIT1 gene (Scheme 2C), which appears to be connected directly to SSB1 (Scheme 2B), could serve here to exemplify possible outcomes. The Rit1 enzyme is known to catalyze the formation of Ar(p) by the addition of a 2Ј-O-ribosyl phosphate group to adenine (52). A typical substrate is represented by adenine 64 of the initiator tRNA-Met(i), which prevents its use during the elongation step of protein synthesis in yeast (68,69). Our data indicated that this PTM was uniquely expressed in the WT strain during heat shock and was both hog1-and ssb1-dependent. When the rit1⌬ mutant was examined (see for example the heat map obtained from rit1⌬ under normal growth conditions, Fig. 1C), no Ar(p) could be detected with or without stress, thus confirming that RIT1 is indeed essential for the production of this PTM. Further, in the context offered by the mutants in the study, this outcome corroborated the putative HOG1 3 SSB1 3 RIT1 pathway as the signaling control for Ar(p) (marked in light blue in Scheme 2B and 2C). Surprisingly, close examination of rit1⌬'s global profile revealed the absence of nine PTMs that were both hog1-and ssb1-dependent (Table II), which potentially places them under the same HOG1 3 SSB1 3 RIT1 regulatory pathway. In addition, three other PTMs that were not hog1-or ssb1-dependent were found to be rit1-dependent. These observations suggest that RIT1 may itself act as an additional hub in the control network. The fact that the majority of rit1-dependent PTMs do not SCHEME 2. A, Established relationships between genes involved in response to hyperosmotic stress and heat shock. B, Putative interactomic network connecting established stress-response genes with PTM biogenetic enzymes. C, Known biogenetic enzymes (red boxes) for selected PTMs of interest (blue). Arrows indicate proven functional relationships. Vertical lines indicate possible functional links. Horizontal lines indicate possible co-expression. Light blue arrows and boxes identify the HOG1 3 SSB1 3 RIT1 pathway corroborated by gene deletion.
involve adenine modification rules out the possibility that Ar(p) itself might serve as their direct biosynthetic precursor. An alternative hypothesis could be that activation of RIT1 may result in the co-expression of other biogenetic enzymes yet to be identified. In this direction, it must be noted that, of the 12 stress-induced PTMs that are rit1-dependent, Ar(p) is the only one with a known biogenetic enzyme. Considering the large number of genes without assigned function in the yeast genome, the functional link with RIT1 activation may provide the basis for further studies aimed at the identification of the missing biogenetic enzymes.

CONCLUSIONS
Monitoring the epitranscriptomic profile of S. cerevisiae as a function of different experimental variables has unambiguously shown that PTM biogenesis is profoundly influenced by environmental conditions. The fact that the production of numerous new PTMs was readily activated by increasing temperature or salt concentration provides very strong evidence for the broad implication of PTMs in stress response. This finding is consistent with the role of selected tRNA modifications that are known to participate in translational recoding triggered by oxidative stress (32,33). However, the majority of the activated modifications observed in our study are not known to engage in these types of mechanisms. The direct comparison of PTM profiles has clearly revealed the presence of stress-specific modifications, as well as nonspecific ones, which could be part of a more general cellular response. These observations support the presence of alternative mechanisms involving PTM activity, which are likely to include the broad participation of diverse classes of RNAs.
With the goal of investigating the relationships between established stress-response signaling and PTM activation, we devised a strategy that involves evaluating the simultaneous effects of gene deletion on the expression of all PTMs in the cell. This system-biology approach was inspired by the wealth of genomic and functional data available for S. cerevisiae. In what could be described as an interactome-walk, genes connecting the putative biogenetic enzymes with known response genes were systematically eliminated to retrace a putative regulatory network inferred from available information. Preliminary results have allowed us to clearly differentiate PTMs that fall under HOG control from those that do not. The HOG pathway consists of a cascade of MAP kinases, which is preeminently responsible for the response to osmostress, but has also tenuous communication with parallel heat shock controls (20, 42, 59 -63). Our results support the presence of numerous integration hubs between these stress-specific pathways, which coordinate the expression of selected PTMs. This type of regulatory architecture could help explain crosstolerance effects, in which preconditioning to a specific type of stress can increase the tolerance to a different one (70). These observations should be considered also in the context of the effects of HOG1 on chromatin remodeling (71,72) and activation of DNA-binding transcriptional activators (73)(74)(75)(76)(77). The facts that HOG1 increases the stability of a well-defined group of osmostress mRNAs (12,13) and induces specific lncRNAs (15) under osmostress leave ample space for the possible existence of feedback mechanisms based on PTM activity. In this direction, we could speculate that chromatin remodeling by stress-induced lncRNAs may affect the expression of genes involved in the general cellular response, as well as PTM regulatory pathways. In turn, the latter could integrate signals from other (likely protein-based) pathways and modify lncRNAs to reduced/enhance their chromatin remodeling activity. Therefore, determining the distribution of PTMs in these classes of RNAs and monitoring their levels under stress will provide the basis for exploring these possible mechanisms.
After establishing a direct link between stress and selected PTMs, the next challenge will be the identification of the actual roles played by such PTMs in cellular response. This important task is severely hampered by the fragmentary knowledge of their distribution outside tRNA and rRNA, and the absence of annotation on the extensive genomic/transcriptomic data available to the public. Toward this goal, we are currently developing technologies capable of determining the class of RNA containing a specific PTM, the identity of the parent species, and the sequence position of the modification. This information will help interrogate the available transcriptomic data, which will lead to testable hypotheses by suggesting possible correlations with known functions associated with the parent RNA. This information will also guide the implementation of strategies based on the mutation of nucleotides susceptible to modification, or silencing of the RNA bearing the modification, aimed at testing their effects on the metabolic and epigenetic state of the phenotype. The results will provide valuable insights into the functions of PTMs in stress response and, more broadly, will be expected to drive a re-evaluation of their overall biological significance. * This work was supported by The RNA Institute of the University at Albany and by the National Institutes of Health (GM064328 -14). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.