Temperature-dependent tRNA Modifications in Bacillales

1. Introduction

More than 100 chemical modifications of the four canonical nucleotides, adenosine (A), cytosine (C), guanosine (G), and uracil (U) have been identified in RNA molecules. All modifications are posttranscriptionally incorporated by a diverse set of chemical reactions including deaminations, isomerizations, glycosylations, thiolations, methylations, and transglycosylations [1] and affect not only A, C, G, and U, but also create hypermodifications such as queuosine or wybutosine from already modified nucleotides [1,2]. The most diverse collection of modifications has been found in rRNAs and tRNAs [3].

With the advent of cost-effective high throughput sequencing, it has become possible to assess modified nucleotides on a genome- and transcriptome-wide scale. The common basis of a broad class of techniques is the use of reverse transcription (RT) that terminates just before certain modifications (e.g., for 1-methyladenosine (m¹A) and 1-methylguanosine (m¹G)) and/or misincorporates nucleotides at very high rates when encountering other modifications (e.g., inosine (I), m¹A, and 6-methyladenosine (m⁶A)). The resulting cDNA is then subjected to high-throughput sequencing, mapped to a reference genome, and computationally analyzed for strong arrests of primer extension (RT-stops) and position-specific mismatch profiles [4,5,6,7,8]. These events constitute a RT signature which varies depending on the modification type and the used reverse transcriptase [4,5,6,7,8]. While some modifications directly affect reverse transcriptase activity or fidelity, other modifications require a preceding chemical treatment to become visible. The treatments can lead to a conversion of the modified nucleotide or induce a block of the primer extension by the reverse transcriptase. As a result, the read-out of the chemically treated nucleotide is different compared to the untreated control. Thus, these modifications can directly be detected by comparing treated versus untreated samples. For example, pseudouridine (Ψ) becomes detectable after treatment with carbodiimides which acylate one of the nitrogen positions of the base, resulting in Ψ-1-cyclohexyl-(2-morpholinoethyl)carbodiimide (CMC) adducts that block RT [9,10,11]. The saturated pyrimidine ring of dihydrouridine (D) can be opened by sodium borohydride (NaBH₄) or under alkaline conditions. As a result, base pairing is impeded, leading to straightforward detection by RT signatures. In addition, 7-methyl-guanosine (m⁷G) reduction by NaBH₄ leads to a depurination in RNA followed by the cleavage of the RNA chain by β-elimination [12].

Mainly, tRNA modifications act as checkpoints for tRNA identity as well as integrity, regulate translation, modulate the structural stability, and ensure the correct folding of the tRNA molecule [13]. They are further involved in structural fine-tuning of local elements and rapidly adapt tRNAs to environmental changes, such as stress and temperature [14,15]. Most chemical modifications modulate the tRNA structure in terms of stability and flexibility [16]. m¹A at position 9 in the human mitochondrial tRNA^Lys, for example, disrupts a Watson-Crick base pair A9-U64 that stabilizes a non-functional rod-like structure and forces the tRNA to adopt the typical cloverleaf secondary structure and the L-shaped three-dimensional functional conformation [17]. Other modifications in the D- and TΨC- loop, carrying D and Ψ, respectively, are involved in tertiary base pairs and enhance the folding of the transcript into the L-shape [18,19]. Pseudouridine, in particular, contributes to the overall stability of tRNA structures with the help of water-mediated bridging interactions between modified bases and the RNA backbone [20,21]. Dihydrouridine, on the other hand, frequently found in the D-arm, is the only non-planar base and thus reduces stacking interactions, conferring a certain flexibility to tRNAs [22].

Initial studies indicate that the control of flexibility versus rigidity indeed plays a role in the temperature adaptation of microorganisms. Accordingly, high levels of D are found in psychrophilic organisms, where flexibility in macromolecules is interpreted as a frequent strategy for cold adaptation [23,24]. In psychrophilic Archaea, tRNAs contain a higher amount of D modifications than their thermophilic counterparts [1,25]. Pseudouridine levels in bacteria can also change as a response to growth temperature [26,27]. However, these analyses investigated the overall change of certain modifications, but did not analyze changes at the resolution of individual tRNA sequences. Hence, it is an interesting question as to how modification levels at individual positions in the tRNAs are adjusted in response to temperature changes to keep these molecules functional.

To investigate a possible correlation of temperature and posttranscriptional modifications in bacterial tRNAs, four closely related Bacillales with different growth temperature profiles (Bacillus subtilis (mesophilic), Exiguobacterium sibiricum, Planococcus halocryophilus (both psychrophilic), and Geobacillus stearothermophilus (thermophilic)) were used as model organisms. These strains were cultured at temperature ranges from 10-70°C, including the individual optimal growth temperatures as well as the minimum and maximum temperatures feasible in the laboratory. Together, this allowed a screen for temperature-dependent, short- term changes in the tRNA modification pattern within the organism and long-term evolutionary adaptation to their particular habitat.

2. Results

2.1. Specific Library Design for tRNA Sequencing and Subsequent Analysis of tRNA Modifications

For a systematic detection of D, Ψ, and m⁷G modifications in the tRNA mapping profiles, individual tRNA samples were chemically treated to convert these modifications to yield a specific read-out in the subsequent sequencing step compared to untreated control samples. From each bacterial culture, small RNA fractions were prepared by high salt precipitation of the high molecular weight RNA. A quality assessment of the small RNA preparations was performed on a BioAnalyzer 2100 device (Agilent) as well as on high-resolution denaturing polyacrylamide gels. Subsequently, individual small RNA preparations were treated with NaBH₄ (for D and m⁷G profiling) and CMTC (for Ψ identification). Treated as well as untreated samples were then used for library preparation and RNA seq on an Illumina MiSeq device (experimental protocol is illustrated in Figure 1). To this end, a 5’-phosphorylated DNA adapter was ligated to the 3’OH ends of the small RNAs. The 3’-end of this adapter was amino-C7-blocked to avoid concatemer formation. To ensure efficient ligation, the adapter was used in a two-fold excess compared to the small RNA preparation.

Next, the ³²P-5’-labeled RT primer ) was hybridized to the sRNA/adapter ligation product. At its 5’-end, this primer carried the cDNA-3’-adapter region for subsequent circularization of the cDNA products by TS2126 RNA ligase (CircLigase). Furthermore, it contained the sequence corresponding to the Illumina P5 adapter (Figure 1). For cDNA synthesis, SuperScript IV reverse transcriptase was used, and the resulting radioactively labeled reaction products were separated and visualized on a denaturing 15% polyacrylamide gel (a representative example is shown in Figure 2A). Bands representing cDNA products were cut out, isolated, circularized and directly used for PCR amplification. A single band migrating at a higher position than the tRNA signals was identified by Sanger sequencing after cloning and omitted from further analysis. To avoid saturation effects and an amplification bias, the PCR reaction was limited to 10 – 12 cycles (Figure 2B, representative example). One PCR primer corresponded to the Illumina P5 flow cell linker sequence, the second one hybridized to the RNA-ligated 3’-adapter and introduced the Illumina sequence (Supplementary Table S1) as well as an index region for identification of the individual library sequences. The indexing allowed for a mixing of the individual amplification products and a combined sequence analysis. A total of 36 libraries per organism were prepared, consisting of three biological replicates, each grown at the three indicated temperatures. For each temperature sample, the following four treatments were performed: untreated sRNA preparation, NaBH4 treatment (detection of D, m⁷G), CMCT treatment and CMCT-negative control (detection of Ψ). Subsequently, the pooled libraries were sequenced on an Illumina MiSeq device (Azenta).

After sample quality assessment and normalization of the resulting sequences, we scanned for RT-stop sites that are significantly enriched (FDR adjusted P value < 0.01) in the mapping profile of the treated libraries compared to the untreated control sample according to the Poisson distribution. Additionally, we distinguished true modification sites from spurious noise. Thus, significant RT-stop sites with a log₂ FC ≤ 1 were filtered out, since we assumed that the modified sites of D, m⁷G, and Ψ showing a higher effect size due to the treatment than that from noise. Further, low read coverage at specific tRNA positions may result in overestimation of the ratio of RT-stops to the total number of reads. According to our observations, treatment-based RT-stop sites have a higher ratio of reads with RT-stops to reads without RT-stops the corresponding position compared to noise. Therefore, only RT-stop sites with a total number of RT-stops ≥ 20 and a percentage of RT-stops ≥ 2 were considered. The thresholds were chosen to be close to the maximum value of potential true positive sites (TPs) at which known sites are consistently predicted, while keeping the number of false positive sites (FPs) as low as possible. Figure 3 shows example data of B. subtilis at 20°C, illustrating both treatments and the controls for each replicate. The replicates show a consistent pattern for each condition, supporting the quality and reproducibility of the data.

To classify FP and TP sites of potential D, m⁷G, and Ψ modifications, knowledge from the RNA modification database MODOMICS [28] and scientific literature [3,5,29,30,31,32,33,34,35] was utilized regarding known bacterial tRNA positions where these modifications are present. The MODOMICS database (https://genesilico.pl/modomics/, accessed on January 3rd, 2024) contains 383 bacterial tRNA sequences, with only four tRNA sequences for G. stearothermophilus and 25 for B. subtilis, but none for P. halocryophilus and E. sibiricum. As only limited information is known about tRNA modifications in these bacteria, tRNA modifications from other bacterial species were also used as a reference. Thus, for the NaBH₄-treated RNA seq data, RT-stop sites for tRNA U positions 16, 17, 20, 20a_, 20b, 21, and 47 were counted as TP D modifications and signals at G46 as TP m⁷G modifications. Upon first inspection, we observed RT-enriched sites at position U8 for several tRNAs, which, due to their strong enrichment, did not appear to arise from noise. However, D modifications are not known to occur at U8 in bacteria; instead, s⁴U modifications are known [28,36,37,38]. To determine whether the modifications were D or s⁴U, we conducted tandem mass spectrometry (MS/MS) analysis in G. stearothermophilus (55°C growth temperature). We successfully identified s⁴U modifications in tRNA^Arg_ACG, tRNA^Arg_CCG, tRNA^Arg_UCU, tRNA^Glu_UUC, tRNA^Leu_GAG, tRNA^Phe_GAA, and tRNA^Trp_CCA (Figure 4), where we also detected significant RT-stop enrichments in their mapping profile (Supplementary Figure S6). Moreover, based on MS/MS analyses, we identified s⁴U modifications for tRNA^Cys_GCA, tRNA^Gly_CCC, tRNA^Leu_UAG, tRNA^Ser_UGA, and tRNA^Tyr_GAU for which we could not find significantly enriched RT-stops (Supplementary Figure S1). Based on this validation through an independent approach, we assume detecting s⁴U at position U8 using our NaBH₄-treated RNA seq data and have classified these sites as TPs.

For the CMCT-treated RNA seq data, Ψ modifications are known to occur in bacteria at U positions 13, 31, 32, 38, 39, 40, 55, and 65. Upon initial inspection, we observed RT-enriched signals at position U60 in several tRNAs (Supplementary Figures S7-S10), where Ψ modifications were previously unknown in bacteria. Through independent MS/MS analysis in G. stearothermophilus (grown at 55°C), we confirmed that these were Ψ60 modifications (Supplementary Figure S2). RT-stops at these positions were counted as TPs Ψ modifications. FPs were counted as type I errors, where all RT-stop sites for Ts and Gs (Gs only in the NaBH₄-treated data) occur at tRNA positions where no modifications are previously known. Type II FPs were counted as all RT-stop sites at C, A or G positions (Gs only in the CMCT-treated data). In summary, this systematic approach allowed us to accurately detect and classify D, Ψ, s⁴U and m⁷G modifications in tRNAs while minimizing false positive signals.

2.2. Variations in tRNA Modification Profiles Among Closely Related Bacillales

To investigate the general modification profiles of the four closely related Bacillales Planococcus halocryophilus, Exiguobacterium sibiricum, Bacillus subtilis, and Geobacillus stearothermophilus, the bacteria were cultivated at their optimal growth temperatures. These temperatures are 20°C for P. halocryophilus and E. sibiricum, 30°C for B. subtilis, and 55°C for G. stearothermophilus. A summary of all profiled modified tRNA positions in the four bacteria is illustrated in Figure 5.

Only the chemically untreated RNA seq data were utilized to identify tRNA modification patterns in the mapping profile that are clearly visible as accumulations of base misincorporations (Supplementary Table S2). In all organisms examined, more than 99% of the reads that could be assigned to tRNA^Arg_ACG clusters have a guanosine instead of the genomically encoded adenosine at tRNA position 34 in the anticodon loop. This is a strong indication of enzymatic deamination of A34 to I34 [5,7,8]. The altered base pairing properties of inosine lead to preferential base pairing with cytosine during cDNA synthesis, resulting in identification as guanosine in the mapping profile of RNA seq data. In the D-stem of several tRNAs of all four bacteria, a more extensive pattern of misincorporation was detected at position A22. In these modified tRNAs, the genomically encoded base is adenosine, which is typically identified as guanosine or thymine in the read pattern, indicating a m¹A22 modification [5,6,7,12]. Only in tRNA^Arg_CCG of P. halocryophilus, a mismatch pattern is found at position G37 of the tRNA anticodon. The underlying modification indicates a methylated guanosine m¹G37 [5,12].

Analyzing the chemically treated RNA seq data revealed numerous strong RT-stop sites (NaBH₄: Supplementary Table S3; CMCT: Supplementary Table S4). For a comprehensive overview of all identified s⁴U, D, m⁷G modifications, as well as type I FPs for each tRNA cluster and bacterium, please refer to Supplementary Figures S3-S6. Using the NaBH₄-treated RNA seq data, several tRNA genes show s⁴U8, D16, D17, D20, D20a, D20b, D47 and m⁷G46 modifications in each bacterium except for E. sibiricum lacking s⁴U and D20b modifications and B. subtilis lacking D47. Additionally, all four bacteria exhibit strong FPs RT-stop signals at G45. These RT-stop sites only occur if m⁷G46 is present.

In the CMCT-treated RNA-Seq data at the optimal growth temperature, all four bacteria show Ψ31, Ψ39 and Ψ55 modifications (refer to the overview of detected Ψ modifications and type I FPs in Supplementary Figures S7-S10). However, Ψ32 was not found in P. halocryophilus, Ψ38 is only present in the two psychrophilic bacteria, and Ψ60 is not present in E. sibiricum. Since Ψ32 and Ψ60 are present at the other growth temperatures considered, it can be assumed that these modifications occur in the respective bacteria.

Additionally, strong RT-stop sites at position U54 are present in many tRNAs across all bacteria. In bacteria, only 5-methyluridine (m⁵U) modifications are known at this position. Unfortunately, a CMCT treatment is not suited to detect such modifications, making it difficult to definitively determine whether it is a Ψ or another modification without further validation. Other FPs in both treatments seem to occur randomly and vary among bacteria.

Even though most of the profiled modifications occur in the tRNA pool of the four bacteria, they vary in the number of affected amino acid families (grouped by the corresponding encoded anticodon triplets) for each modified tRNA position. For example, while m¹A22 modifications were found in the same tRNAs (tRNA^Tyr_GUA, tRNA^His_GUG, tRNA^Glu_UUC, tRNA^Glu_UUG, tRNA^Cys_GCA, tRNA^Leu_UAA, tRNA^Leu_GAG, tRNA^Leu_CAA, and tRNA^Ser_GCU) in all four bacteria, m¹A22 modifications in other tRNA^Ser and tRNA^Leu molecules occur specifically for each bacterium, depending on the encoded anticodon (Supplementary Figure S11).

2.3. Specific tRNA Modifications Are Temperature-Dependent

To investigate tRNA-dependent temperature adaptation mechanisms of closely related psychrophilic (P. halocryophilus and E. sibiricum), mesophilic (B. subtilis), and thermophilic (G. stearothermophilus) Bacillales, different individual temperatures between 10 and 70°C were cultivated. This comprises in detail the optimal growth temperature, as well as (as far as feasible in the laboratory) minimum and maximum temperature: P. halocryophilus (10°C, 20°C, and 30°C), E. sibiricum (10°C, 20°C, and 30°C), B. subtilis (20°C, 30°C, and 37°C), and G. stearothermophilus (40°C, 55°C, and 70°C).

Regarding the modifications identified through base misincorporations in the chemically untreated RNA seq data, no strong difference was observed in the fraction of base mismatches for these specific modifications across the investigated bacteria at their respective specific growth temperatures (Supplementary Figure S12), suggesting no stronger temperature dependence of the m¹A22, I34, and m¹G37 modifications.

Analysis of D, m⁷G, and Ψ tRNA modifications revealed consistent patterns in P. halocryophilus, E. sibiricum, and B. subtilis within the different temperatures assessed, with similar levels of modified tRNA clusters observed for modifications such as D17, D20a, Ψ39 and m⁷G46 (Table 1). Interestingly, Ψ38 was uniquely identified in the two psychrophilic bacteria. Notably, a distinct trend was observed in G. stearothermophilus, where the number of modified tRNAs at specific positions increased with rising temperatures (Table 1). For instance, the abundance of modifications at D17 (40°C: 1; 55°C: 12; 70°C: 13), D20 (40°C: 6; 55°C: 18; 70°C: 19), and Ψ55 (40°C: 9; 55°C: 21; 70°C: 29) exhibited a clear temperature-dependent increase. In contrast, modifications like Ψ32 and m⁷G46 remained relatively stable across temperature variations in G. stearothermophilus. Furthermore, a striking observation was the significantly higher occurrence of s⁴U8 modifications in G. stearothermophilus (>24) compared to other psychrophilic and mesophilic bacteria (≤4).

2.4. Performance Validation

For analysis performance validation, we compared the tRNA modifications available in the MODOMICS database for each listed tRNA gene, including m¹A, m⁷G, s⁴U, D, and Ψ modifications, with our results for G. stearothermophilus and B. subtilis (Figure 6A, Supplementary Table S5). For the specific tRNAs available in the database for G. stearothermophilus (N = 4) and B. subtilis (N = 21, N = 4 not found with tRNAscan-SE), we were able to reproduce 89% of the modifications in our analyses, while only 10 out of 87 modifications could not be found. Most of these 10 modifications are present in our data as weak signals, so they were filtered out by our thresholds as we could not clearly distinguish them from noise. However, we discovered 35 additional modifications for these tRNA genes that are not listed in the database.

Additionally, to verify the validity of our findings using a different analysis method, we used the previously published mim-tRNAseq tool [39] for our B. subtilis (30°C) RNA seq data (Supplementary Table S6). Considering the affected amino acid families grouped by their corresponding encoded anticodon triplets and modified tRNA position (counted uniquely), a total of 61 modifications were identified by both mim-tRNAseq and our method (Figure 6B). In addition, 57 modifications were detected only by our analysis, while 14 modifications were found only by mim-tRNAseq. The number of type I FPs (mim-tRNAseq: 58; own: 63) is approximately the same for both methods (Supplementary Table S6).

In conclusion, our performance validation demonstrated the effectiveness of our method by identifying a large number of previously known modifications for specific tRNA genes. Our comparative study with mim-tRNAseq confirmed our results and showed that our method can reliably detect both known and novel modifications.

3. Discussion

This study aimed to investigate whether variations in tRNA modifications represent an adaptation of bacteria to different growth temperatures. To address this inquiry, we cultivated four closely related Bacillales species P. halocryophilus, E. sibiricum, B. subtilis, and G. stearothermophilus under three different temperature conditions (minimal, optimal, and maximal). We combined an RNA seq approach with chemical pre-treatment of the tRNA samples and developed an analysis strategy that allows for systematic profiling of D, s⁴U, and m⁷G (NaBH₄-treatment) as well as Ψ (CMCT-treatment) modifications in the transcriptome at single-nucleotide resolution. Based on our statistical quantification of the data, we detected a variety of RT-stop sites which are significantly and strongly enriched in the treated samples compared to the untreated control. A significant challenge was that the tRNA treatments were designed to detect up to three different modifications in the resulting reads. In contrast to gene expression analysis, where statistical robustness and performance can be increased by aggregating reads from an entire gene, measurements of tRNA modifications comprise only single nucleotides [40]. Another challenge was to distinguish true modification sites from noise. Such RT-stop sites occur since the reverse transcriptase is sensitive to certain structural features of the RNA template [41,42] or to other modifications which are not specifically enriched by the treatment. Complex processing of RNA, genomic misalignments of sequencing reads, and technical errors of the sequencing platform contribute to incorrect signals. We assumed that - due to the chemical tRNA treatments - the modified sites result in a considerably higher signal over noise ratio. The use of a specifically adapted FC cutoff is, therefore, suitable. In addition, low read coverage at specific tRNA positions may result in overestimation of the ratio of RT-stops to the total number of reads. In both treatments, these signals are particularly frequent in the TΨC-arm, D-arm, as well as the 5’-acceptor stem of tRNAs. In addition, we observed that the SuperScript IV reverse transcriptase exhibited an increasing tendency to incorporate incorrect bases in the 3’-range and terminate prematurely, leading to background noise in the data. Such signals were filtered out by only considering RT-stop sites that exhibited a minimum number and percentage of RT-stops at the respective position. We were not able to completely remove the background noise without accepting an increased false negative rate by more stringent parameter settings. As expected, an approximately equal numbers of type I FPs were observed using the mim-tRNAseq tool. However, such FP RT-stop sites showed variability across different bacterial species and growth temperatures, indicating a lack of clear patterns. As ignoring the impact of FP detection likely leads to data misinterpretation, a careful consideration and mitigation of background noise and its incorporation into threshold determination may enhance the accuracy and reliability of tRNA modification profiling studies. To determine the optimal cutoffs, we classified RT-stop sites as TP modifications or FPs based on established information from the RNA modification database MODOMICS [28] and relevant scientific literature [3,5,29,30,31,32,33,34,35]. It is important to note that this classification was not specific to individual tRNAs, as our focus was on less well-characterized bacteria. Due to the lack of comprehensive characterization of the investigated bacteria, this method is limited in its ability to detect unexpected modifications that may occur during a temperature shift, as the true positive positions of modified tRNA were identified based on existing knowledge.Irrespective of potential false negatives at unexpected sites, our data provide insights in temperature dependent changes at those sites where we observe modifications with high confidence.

Upon initial examination of the NaBH₄-treated RNA seq data, we identified RT-enriched sites at position U8 in several tRNAs, suggesting a specific modification rather than random noise. It is noteworthy that s⁴U modifications are typically associated with U8 in bacteria and archea [28,36,37,38]. Our mass spectrometry-based analysis (Figure 4) validated the presence of s⁴U8 modifications in G. stearothermophilus. Since we detected both MS/MS signals and enriched RT stops in the same set of 8 tRNAs, we hypothesize that we could also identify s⁴U8 modifications with the NaBH₄-treated data.

In addition to identifying the general tRNA modification patterns in the investigated bacteria, our study aimed to detect temperature-induced differences in tRNA modifications.The analysis revealed numerous strong RT-stop sites and base misincorporation accumulations in the mapping profiles, indicating the presence of various modifications such as s⁴U, D, m¹A, m⁷G, m¹G, I, and Ψ (Figure 5). Despite most of these modifications occurred in the tRNA pool of the four species, the number of each modified tRNA position varied in the individual isoacceptor families. This provided a unique modification profile for each bacterium, even though they are closely related (Supplementary Figure S11). It is important to note that the absence of certain modifications does not necessarily imply their non-existence but could be attributed to weaker expression levels or exclusion by our stringent cutoff criteria. For example, we did not observe strong RT-stops for Ψ38 in tRNAs of B. subtilis, despite previous reports of its presence [12,43]. Further, m¹G37, involved in frame shift suppression [44] was only found in P. halocryophilus, while it is also described in B. subtilis [35]. Further validation methods such as deletion of RT-stop signatures by enzymatic demethylation of m¹G by AlkB-type enzymes [45] or MS/MS analysis [46,47] could help confirm the presence of these modifications. Similarly, slight variations in the number of modified tRNA clusters within a species across different cultivation temperatures should not be overinterpreted as the variations may reflect differences in modification intensity rather than distinct modifications. Thus, our cutoff criteria allowed us to infer the presence of tRNA modifications with minimal intensity levels. In addition, strong RT-stop sites at position U54 are observed in numerous tRNAs across the bacterial species. In bacteria, only m⁵U54 (ribothymidine) modifications have been documented [3,28,35,48]. Unfortunately, it is not known whether a CMCT treatment can effectively detect such modifications, posing challenges in definitively classifying the RT-stop sites without further validation. However, the presence of robust RT-stop sites at position U54 in almost all tRNAs aligns with the listing of m⁵U modifications in the MODOMICS database, reinforcing the likelihood that this type of modification has been detected in this context.

The analysis of temperature-dependent tRNA modifications revealed intriguing patterns in psychrophilic and mesophilic bacteria, with consistent trends observed in these organisms. Interestingly, the thermophilic bacterium G. stearothermophilus exhibited an increase in the number of tRNAs modified at specific positions as temperatures increased. This observation suggests a temperature-dependent regulation of tRNA modifications. Of particular interest is the significantly higher occurrence of s⁴U8 modifications in G. stearothermophilus compared to non-thermophilic bacteria. Previous studies have highlighted the role of s⁴U modifications in stabilizing tRNA structures and influencing translation efficiency [15,49,50]. It could be further demonstrated that s⁴U8 is conserved in bacterial and archaeal tRNAs and plays a role as a UV protective [15,50]. The increased presence of s⁴U8 in G. stearothermophilus could be interpreted as an adaptive strategy to maintain tRNA stability and functionality under extreme high-temperature conditions, potentially serving as a factor for the thermotolerance of this organism. Furthermore, the presence of the Ψ55 modification has been identified as a control mechanism for thermal adaptation of tRNA functionality [15,51]. The presence of Ψ55 limits the incorporation of modifications that provide thermal stability, whereas its absence promotes the addition of such modifications, enhancing thermal stability. Originally identified as a modification that stabilizes tRNA structure, Ψ55 can also enhance flexibility and functionality of tRNAs at lower temperatures. While the investigated psychrophilic and mesophilic Bacillales did not show significant trends for Ψ55 modifications concerning temperature differences, G. stearothermophilus exhibited a substantial increase in Ψ55 from 9 modified tRNA clusters at minimal growth temperature to 29 modifications at maximal growth temperature. This suggests that Ψ55 may play a crucial role in the thermal adaptation of G. stearothermophilus by modulating stability and flexibility of tRNAs at varying temperatures.

It has been reported that some psychrophilic microorganisms have a higher dihydrouridine content in their tRNAs compared to mesophilic and thermophilic organisms environments [15,22,24,25,52]. This modification represents the only base that adopts a chair conformation. All other bases in nucleic acids are strictly planar and readily form base-stacking interactions that stabilize and rigidify the surrounding structure. Due to its non-planar structure, dihydrouridine interferes with base stacking and leads to an increased local flexibility in tRNAs, a feature that may represent an important adaptation to cold environments [15,22,24,25,52]. However, we cannot fully confirm this hypothesis, as we do not observe significant differences between psychrophilic and mesophilic bacteria. Nevertheless, the number of D modifications in tRNAs from psychrophilic and mesophilic bacteria is higher than in the thermophilic bacterium.

Only a few tools have already been published for the detection of tRNA modifications using treatment-based RT-stops in sequencing mapping profiles, for example see [39,53]. To compare our findings with those of an already established tool for such analyses, we utilized the mim-tRNAseq tool. The mim-tRNAseq methodology relies on prior knowledge to classify true positives, like the approach adopted in our analyses. Our method demonstrated high effectiveness, as we were able to identify over 81% of the modifications (grouped by amino acid families and their encoded anticodons) found by mim-tRNAseq in B. subtilis, along with an additional 57 modifications (Figure 6B). Furthermore, we successfully identified 89% of the specific modifications listed in the MODOMICS database for G. stearothermophilus and B. subtilis, as well as detecting an additional 35 modifications for specific tRNA genes (Figure 6A). This comprehensive coverage of known modifications highlights the robustness and reliability of the obtained data at single-base resolution. The ability to detect both previously known and novel modifications underscores the potential of our method for advancing the understanding of tRNA biology.

4. Materials and Methods

4.1. Cultivation of Bacteria

B. subtilis strain DSM 10, E. sibiricum strain DSM 17290, P. halocryophilus strain DSM 24742, and G. stearothermophilus strain DSM 22 were purchased at the German Collection of Microorganisms and Cell Cultures (DSMZ) and grown in the corresponding recommended liquid medium. Cultivation was carried out at various growth temperatures in a culture shaker at 180 rpm and a volume of 50-100 ml. Cell growth was photometrically monitored and cells were harvested by centrifugation for 10 min at 4°C and 3,500 g once they had reached the mean log phase. According to organism and growth temperature, OD600 was between 0.31 and 0.84.

4.2. Isolation of Total RNA and Small RNA Enrichment

Bacterial cell pellets were resuspended in TRIzol^TM. Cell disruption was performed using a FastPrep homogenizer (MP Biomedicals, Germany) for 40 s at 6 m/s, followed by incubation at room temperature for 5 min. 200 μl chloroform was added, mixed for 15 s and incubated for 3 min at room temperature. After centrifugation for 15 min at 12,000 g and 4°C, the supernatant was transferred to a sterile reaction tube, mixed with an identical volume of isopropanol, and incubated on ice for 15 min. RNA was pelleted by centrifugation and washed twice with 1 ml ethanol (70 %). After air drying at 30°C, the total RNA was dissolved in 200 μl sterile water. For small RNA enrichment, the total RNA preparation was mixed with 1/10 volume of 5 M NaCl and 50 % PEG 8,000, and incubated for 30 min at -20°C. After centrifugation for 30 min at 10,000 g and 4°C, the procedure was repeated with the supernatant. The resulting solution was mixed with 3 volumes of ethanol and incubated overnight at -20°C. Small RNAs were precipitated by centrifugation and resuspended in water.

4.3. NaBH₄-Treatment for the Detection of Dihydrouridine and 7-Methylguanosine

Treatment of the small RNA fraction was done as described [54]. 2-4 µg small RNA were adjusted with water to a volume of 50 µl, and 5 µl of a freshly prepared NaBH₄ solution (100 mg/ml in 0.01 N KOH) was added. The mixture was incubated on ice for 30 min. The reaction was stopped by the addition of 5 µl 0.01 N acetic acid followed by ethanol precipitation. The RNA was resuspended in sterile water and kept at -20°C until use.

4.4. CMCT-Treatment for the Detection of Pseudouridine

The procedure is based on the protocol as described [12] with slight modifications. Briefly, 2-4 μg small RNA were mixed with 30 μl of freshly prepared 1-cyclohexyl-(2-morpholinoethyl)carbodiimide metho-p-toluene sulfonate (CMCT) solution (165 mM CMCT, 50 mM Bicine, pH 8.0, 7 M urea, 4mM EDTA) or CMCT buffer (solution without CMCT) as control and incubated for 20 min at 37°C. The reaction was stopped by ethanol precipitation in the presence of 1 µg glycogen. The RNA pellet was resuspended in 40 µl 50 mM Na₂CO₃ solution (pH 10.3) and incubated for 2 h at 37°C. After another ethanol precipitation, the pellet was dissolved in sterile water.

4.5. Adenylation of Adapter Oligonucleotide and Ligation to Small RNA Samples

500 pmol of a 3’-terminally blocked adapter oligonucleotide (5’-pTGGAATTCTCGGGTGCCAAGG-amino-C7-3’) was mixed with 50 pmol of recombinantly expressed and purified TS2126 RNA ligase (CircLigase) in 20 μl adenylation buffer (50 mM MOPS, pH 7.5, 10 mM KCl, 5 mM MgCl₂) containing 2.5 mM MnCl₂, 1 mM DTT, and 50 µM ATP and incubated for 2 h at 60°C as previously described [55,56]. The reaction was stopped by heat inactivation at 80°C for 10 min. Adenylated adapter was purified by phenol/chloroform extraction and ethanol precipitation and resuspended in water at a concentration of 20 µM. For ligation to the 3’-end of the prepared small RNA, 20 pmol of the adenylated adapter were incubated with 10 pmol sRNA sample in 10 µl 50 mM Tris/HCl, pH 7.5, 2 mM MgCl₂, 1 mM DTT, 0.4 mM ATP, 20% PEG 8,000 and 14 pmol recombinantly prepared T4 RNA ligase 2, truncated KQ (T4Rnl2KQ). After incubation at 25°C for 20 h, the reaction was stopped by heat inactivation at 65°C for 20 min.

4.6. Reverse Transcription

10 pmol phosphorylated RT primer (spiked with 1 pmol of 5’-³²P-labeled RT primer for monitoring) were added to the adapter-ligated small RNA sample. Reverse transcription was performed with SuperScript IV according to the manufacturer's instructions (Thermo Fisher Scientific, USA). The reaction was stopped by the addition of 1 μl 5 N NaOH and incubation at 95°C for 3 min. The solution was neutralized by adding 1 μl 5 N HCl. 8 μl 3x RNA loading buffer were added to the mixtures. Reaction products were separated on a 15% denaturing polyacrylamide gel, visualized by autoradiography, cut out using a sterile blade and isolated according to Roth et al. [57]. As the RT primer (Supplemental Table S1) was used for circularization of the complementary DNA (cDNA), it contained two 18-atom hexa-ethylene glycol spacers (Sp18) to avoid rolling circle amplification during PCR [58,59,60]. In addition, the primer contained a binding site for the Illumina small RNA library primer for preparation of the sequencing library.

4.7. Circularization, Amplification and Indexing of cDNA

Gel-purified cDNA was incubated at 60°C for 2 h in the presence of 17 pmol TS2126 RNA ligase (CircLigase) in a volume of 20 µl adenylation buffer supplemented with 5 mM MnCl₂, 20 mM DTT, and 0.1 mM ATP. The ligation reaction was stopped by heat inactivation at 80°C for 10 min. The individual sequencing libraries were prepared by PCR according to the Illumina TruSeq smallRNA library instructions. 5 μl of circularized cDNA were added to a reaction volume of 50 μl containing 1x Phusion HF buffer, 0.2 mM dNTPs, 0.5 μM primer containing Illumina P5 sequence, 0.5 μM index primer, and 1 U of Phusion DNA polymerase (Thermo Fisher Scientific, USA). Each reaction was split and amplified in 10 to 12 cycles. PCR products were purified using the QIAquick PCR Purification Kit (Qiagen, Netherlands) according to the manufacturer’s instructions and the concentration was determined photometrically. P5-containing primer and primer containing Illumina P7 sequence are listed in Supplementary Table S1.

4.8. High Throughput tRNA Sequencing

Single-end RNA sequencing (RNA seq; 1x100 bp) of the cDNA libraries was performed on an Illumina MiSeq device using the MiSeq Reagent Kit v3. Three biological replicates were sequenced for the NaBH₄- and CMCT-treatment as well as for the corresponding untreated negative control for each bacterium and cultured temperature (in total 36 RNA seq libraries per species).

4.9. Annotation of tRNAs

Cytosolic tRNAs were annotated with tRNAscan-SE v2.0 [61] using the default model for bacteria. Obtained tRNA sequences were aligned to the tRNAdb [36] database using BLAST v2.4.0 [62] to fit individual tRNA secondary structure annotation to the corresponding standard tRNA model [63].

4.10. Mapping of tRNA Reads

To trim adapter sequences and low-quality portions of raw reads, Cutadapt v1.16 [64] was applied with a quality cutoff of 25 and a maximum error rate of 0.15. Only reads of 8 to 95 nts length (after trimming of adapter and low-quality bases) were captured. Read mapping and filtering were performed following the best-practice tRNA mapping strategy as described [65]. In brief, the reads were mapped against a modified target genome in which known tRNA loci were masked and instead intronless tRNA precursor sequences were appended as artificial “chromosomes”. In a first pass, reads displaying specific precursor hallmarks were filtered out. In the second pass, the mature tRNA reads were mapped against mature tRNA sequence clusters which assembled identical tRNAs helping to overcome the multi-copy nature of tRNA genes. Consequently, only uniquely mapped high confidence tRNA reads were used for downstream analyses. The genomes of the four bacteria of the corresponding strain were obtained from the database of the National Center for Biotechnology Information (NCBI) [66].

4.11. Library Normalization

For direct comparisons of control and treated libraries, mapped data were scaled library- and replica-wise. Library-wise normalization was performed by scaling the number of mapped reads for each tRNA cluster position to the total sample number of reads. For replica-wise normalization, the mean value of replicates was calculated for each genomic tRNA position.

4.12. Identification of Significant Base Misincorporation Sites

For detection of tRNA modifications directly visible as accumulation of nucleotide mismatches in mapped reads, bcftools v1.15.1 [67] was used. First, the mpileup command of bcftools was applied to generate genotype probabilities and read coverage at each genomic position. In a second step, the variants were read out using the call command with the -m option to determine rare variants to be included. Variants below the Phred-scaled confidence threshold of 20 and with a coverage of less than 10 reads were excluded.

4.13. Profiling Position-Wise RT-Stops

To identify whether a particular position in a tRNA was significantly altered by chemical treatment, the difference in the number of normalized RT-stop sites in the treated versus untreated libraries was determined. Drawing on the framework used in differential gene expression analysis, a fold change (FC) was calculated at first to differentially quantify effect sizes in RT-stop coverage between treatment and control samples and test for each tRNA position the null hypothesis that there is no significant change between both conditions. The FC for each tRNA and position n was calculated as

$$F{C}_n=S{top}_{n+1}^{+}/\left(\gamma {Stop}_{n+1}^{-}+\alpha \right)$$

(1)

where Stop⁺ and Stop^- are the number of RT-stops of the treated and untreated libraries, respectively. In case of modified tRNAs, the RT terminates mostly one position upstream of the modified nucleotide. Thus, values of ${Stop}_{n+1}^{+}$ and of ${Stop}_{n+1}^{-}$ pertain to the chemical modification at position n since the RNA template is read in 3’ to 5’ direction during RT. For between-sample normalization, the untreated library was scaled by the ratio $\gamma$ of total number of reads mapped per tRNAs cluster as scaling factor given by

$$\gamma =\mathsf{\Sigma }Read{s}_{clust}^{+}/ \mathsf{\Sigma }Read{s}_{clust}^{-}.$$

(2)

To handle complications arising from zero counts, a pseudocount of α = 1 was added.

We do not expect or observe substantial overdispersion and hence model the sampling of read ends (caused by RT-stops) as a Poisson process [68] as second step. Since we expect an enrichment of RT-stops only in the treated samples, the RT-stops derived from the untreated libraries were used to estimate the expected background signal as $\lambda =\gamma \backslash R{T}_{n+1}^{-}$. Since expression values are different for individual tRNAs (more precisely for the clusters of tRNAs with nearly Indistinguishable sequences), $\lambda$ is estimated separately for each tRNA cluster. This additionally reduced the influence of potential outliers. Subsequently, the P value for the observation of RT-stops k (including the pseudo number) was calculated for each tRNA position as follows:

$$P\left(x\ge k\right)=P_{\lambda }\left(k\right)=1-\left({\displaystyle \frac{{\lambda }^k}{k!}}{e}^{-\lambda }\right)$$

(3)

To account for multiple testing, the false discovery rate (FDR) was calculated by applying the Benjamini-Hochberg procedure [69]. Only tRNA positions having an FDR adjusted P value of less than 0.01 were considered as significantly modified.

4.14. Choice of Cutoff Values to Reduce Noise

However, regardless of filtering, a large number of presumably spurious noise remain. To further reduce the latter three different parameters were considered:

• Fold change (FC): A sufficiently large effect size (FC) is required, i.e., position-specific increase of the chemically introduced RT-stops in treated compared to the untreated library.
• Total number of RT-stops: To reduce signals with an overestimated log₂ FC, which can arise at positions with a low read coverage even at a few RT-stops.

Percentage of RT-stops: Treatment-based RT-stops have a higher ratio of reads with RT-stops to reads without RT-stops at the corresponding position compared to noise.

Therefore, an appropriate percentage of RT-stops helps to reduce noise. The thresholds were set as follows: FC > 1; total number of RT-stops ≥ 20; and percentage of RT-stops ≥ 2.

4.15. Outcome Comparisons with mim-tRNAseq

For the comparison study, mim-tRNAseq v1.1.6 [39] was utilized. A misinc-thresh threshold (total misincorporation rate at a position in a cluster used to call unannotated modification) of 0.005 was used to reduce noise. All other parameters were set according to default settings.

4.16. MS/MS Analyses

4.16.1. tRNA Cyanoethylation

tRNAs were chemically treated with acrylonitrile to target 4-thiouridine (s⁴U) modifications. The acrylonitrile derivatizations (Ces⁴U for s⁴U) result in a nucleotide mass increase of +53.0 Da, which can be readily detected using mass spectrometry (MS) [70,71].

4.16.2. tRNA Isolation by 2D-Gel Electrophoresis

tRNA isoacceptors were isolated using two-dimensional gel electrophoresis as described [72]. In short, tRNA was briefly heated to 90°C. Then, samples were separated in a first-dimension 12.5% polyacrylamide gel under denaturing conditions containing 1x TBE and 8 M urea, followed by a second dimension under semi-denaturing conditions using a 20% polyacrylamide gel, TBE 1x and 4 M urea at room temperature. Gels were stained with an ethidium bromide solution (10 µg/l) for 10 min. Spots containing isolated tRNA isoacceptors were visualized and excised under UV light (302 nm).

4.16.3. In-Gel RNase Digestion of tRNAs

Gel spots containing tRNAs were desalted by at least eight washes with 100 µl 200 mM NH₄AcO and dried under vacuum. For RNase T₁ hydrolysis, gel pieces were rehydrated by 20 µl of 0.1 U/µl of RNase T₁ or by 20 µL of 0.01 U/µL of RNase A (ThermoFisher Scientific) in 100 mM NH₄AcO (pH not adjusted). For RNase U₂ treatment, spots were digested by using 50 µl of recombinant RNase U₂ at 0.3 ng/µl in 50 mM NH₄AcO (pH 5.3) and incubated for 45 min at 65°C. After digestion, supernatants were dried under vacuum.

4.16.4. Nano Liquid Chromatography-MS/MS of RNA Oligonucleotides and Analysis

Pellets containing RNase digestion products were resuspended in 3 µl of milli-Q water and separated on an Acquity peptide BEH C18 column (130 Å, 1.7 µm, 75 µm × 200 mm) using a nanoAcquity system (Waters). The analyses were performed with an injection volume of 3 µl. The column was equilibrated in eluant A containing 7.5 mM TEAA (Triethylammonium acetate), 7.0 mM TEA (Triethylammonium) and 200 mM HFIP (hexafluoroisopropanol) at a flow rate of 300 nl/min. Oligonucleotides were eluted using a gradient from 15% to 35% of eluant B (100% LC/MS Grade methanol) for 2 min followed by an increase of buffer B to 50% in 20 min. MS and tandem mass spectrometry (MS/MS) analyses were performed using SYNAPT G2-S (hybrid quadrupole time-of-flight mass spectrometer) from Waters Corporation. All experiments were performed in negative mode with a capillary voltage set at 2.6 kV and a sample cone voltage set at 30 V. Source was heated to 130°C. The samples were analyzed over an m/z range from 550 to 1,600 for the full scan, followed by fast data direct acquisition scan (Fast DDA). Collision-induced dissociation (CID) experiments were achieved using Ar. All MS/MS spectra were deconvoluted using MassLynx software (Waters) and manually sequenced by following the y and/or c series (w ions were also used when sequencing was difficult or to confirm a sequence). tRNA identification was done bs comparisons with the genomic sequences.

5. Conclusions

In conclusion, this study provides insights into both short-term temperature-dependent changes and long-term evolutionary adaptations in tRNA modification patterns within bacterial organisms through specific modifications at single-nucleotide resolution. Our established computational pipeline for identifying chemical modifications in RNA sequencing data proved highly effective, requiring minimal parameters, and showing promise for accurate profiling studies. Importantly, we were able to construct a comprehensive tRNA modification profile for the four bacteria, particularly enhancing our understanding of E. sibiricum and P. halocryophilus where limited information was previously available. Furthermore, we successfully identified Ψ at position U60 in the four bacteria through both MS/MS and our RNA seq analyses, revealing the presence of Ψ at this position in bacteria where it was previously unknown. Our results indicate that the thermophilic G. stearothermophilus displayed a clear temperature-dependent increase in specific tRNA modifications (D17, D20, D20a, Ψ39, and Ψ55) at higher temperatures, suggesting a regulatory mechanism linked to temperature variations. The elevated presence of s⁴U8 and Ψ55 modifications in G. stearothermophilus compared to non-thermophilic bacteria may act as adaptive strategies for maintaining tRNA stability and functionality under extreme high-temperature conditions. Moreover, the higher abundance of D modifications in tRNAs from psychrophilic and mesophilic bacteria implies an essential adaptation to cold environments by enhancing local flexibility in tRNAs.