Quantitative mapping of pseudouridines in bacterial RNA

Domain: chemistry
Published: 2026-02-26
Editorial outcome: minor_revision

Abstract
RNA pseudouridylation is one of the most prevalent post-transcriptional modifications, occurring universally across all organisms. Although pseudouridines have been extensively studied in bacterial tRNAs and rRNAs, their presence and role in bacterial mRNA remain poorly characterized. Here, we used a bisulfite-based deep sequencing approach to provide a comprehensive and quantitative measurement of bacterial pseudouridines using E. coli, to provide proof of concept. We identified 1,954 high-confidence sites in 1,331 transcripts, which is 27 times above previous estimates and representing almost 30% of the transcriptome. Furthermore, pseudouridines were significantly associated with mRNA stability and enriched in transcripts associated with secondary metabolite production and adaptation to diverse environments. Finally, we mapped pseudouridines in oral microbiome samples of human subjects, demonstrating the broad applicability of our approach in complex microbiomes. This way, we observe that, although uridines are required for modification, mRNAs from GC-rich bacteria harbored more pseudouridine sites than AT-rich genomes in our dataset. Altogether, our work highlights the advantages of mapping bacterial pseudouridines and provides a tool to study posttranscription regulation in microbial communities.

Introduction
Pseudouridylation is one of the most abundant nucleotide modifications present in all domains of life1,2. Pseudouridines (Ψs) are found in rRNA, tRNA, and other non-coding RNA where they enhance base-pairing, RNA stability, and influence translation fidelity3,4. Recent transcriptome-wide studies in humans have further identified Ψs in eukaryotic mRNAs where pseudouridylation is altered in response to stress which suggests a regulatory role in eukaryotes5.

In many bacteria species such as E. coli, modification is carried out by eleven different pseudouridine synthase (PUS) enzymes ─ five of which specifically modify rRNA targets (RluB, RluC, RluD, RluE, and RsuA)6,7. On the other hand, TruA, TruB, TruC, and TruD explicitly modify tRNAs, while RluA and RluF pseudouridylate both rRNA and tRNA6. Although, lack of Ψ in eukaryotic rRNA severely impacts ribosomal function, a mutant E. coli strain devoid of any Ψ in the ribosomal subunit did not show any major effect on growth, decoding and ribosome biogenesis7. Among the tRNA PUS enzymes, TruA appears to be the most important. Mutations in TruA, which is responsible for pseudouridylation in anticodon loop, affects the translation machinery8. In particular, tRNA lacking Ψs in this loop, are unable to proceed at aminoacyl-tRNA transfer step due to lack of stability in the mRNA-tRNA complex8.

Most in vivo studies investigating the impact of pseudouridylation in bacteria have been focused on rRNA and tRNA. This is partly due to the unavailability of experimental toolset to investigate the distribution and function of Ψ in mRNA. Previous in vitro studies have demonstrated that replacement of uridine with Ψ in bacterial mRNA can impede amino acid addition and increase the occurrence of amino acid substitutions9. Similarly, Ψ at stop codons (UAA, UAG, or UGA) can enable ribosomal readthrough of the modified stop codon10. Based on these observations, Ψs in mRNAs are believed to influence protein translation in bacteria, but yet to be investigated in vivo.

There has been an attempt to identify pseudouridylation sites in bacterial mRNA at base-level resolution11. The authors relied on Pseudo-seq, a transcription-wide approach for Ψ profiling12. However, Pseudo-seq is known to suffer from low sensitivity and typically identifies very few pseudouridylation sites. In fact, newer methods that uses bisulfite (BS) to label Ψ (e.g. BID-seq13,14 and PRAISE15) have recently been shown to identify ~23 times more Ψ sites than Pseudo-seq15. Therefore, new approaches are required to evaluate the assumed widespread distribution of Ψ in bacterial mRNAs.

Furthermore, studying mRNA pseudouridylation in microbiomes may help elucidate the role of a set of genes under a given condition. Typically, analysis of gene expression of the microbiota is performed using RNA-seq (i.e. metatranscriptomics). In this case, mRNA level in sequence data is used as an indirect proxy for protein abundance. Given that mRNA pseudouridylation can potentially influence translation13, conventional RNA-seq analysis may not accurately recapitulate the protein translation landscape of the microbial community. Consequently, new methods capable of incorporating this widespread posttranscription modification will provide a more accurate representation of protein dynamics. To date, there is no study that has examined mRNA pseudouridylation or posttranscriptional regulation in microbiomes.

In this work, we adapted a BS-based approach to identify Ψ sites in E. coli RNA at base resolution and identified Ψ sites in one-third of mRNAs. We also establish that Ψ in protein-coding RNAs enhances transcript stability. Lastly, we demonstrate the applicability of our method in complex microbiome samples. Altogether, our work provides a tool to study posttranscription regulation in isolate bacteria and complex microbial communities.

Results
A BS-based approach to identify Ψ positions in bacterial RNA
BS-based approaches were recently developed to identify the positions of Ψs in human mRNA13-15. However, due to the instability of prokaryotic mRNAs and absence of poly-A tails, these methods are not readily applicable to bacteria. In BS methods, treatment of mRNA with BS generates BS-Ψ adduct, which subsequently induces deletion at Ψ sites during reverse transcription (RT). By mapping these single base deletions with untreated samples (which have no deletions), the exact location of Ψs can be identified (Fig. 1). We adapted the BS-based approach with modifications and applied it to E. coli samples.

We began by extracting total RNA and enriching for non-rRNAs through ribodepletion. Next, we performed fragmentation and split samples into two halves, with one half treated with BS (Fig. 1). After RT and sequencing, we mapped reads to the E. coli genome, realigned reads to increase sensitivity for reads that contain gaps16 and retrieved the coverage of each nucleotide. Based on thresholds from previous studies13,15, we considered a uridine position as pseudouridylated if (i) the coverage depth in both BS-treated and untreated samples is ≥ 20; (ii) deletion rate is ≥ 5% in BS-treated samples but less than 1% in untreated samples; and (iii) deletion count is above 5 in BS-treated samples. Therefore, the deletion ratio is proportional to Ψ levels. For instance, a deletion ratio of 0.5 suggests 50% of transcripts harbor Ψ in a site.

Validation of BS-based profiling of bacterial Ψs
To validate our approach, we used a wild-type (WT) E. coli and a mutant strain with knockout deletions in all seven rRNA PUS enzymes7 (rluA, -B, -C, -D, -E, -F, and rsuA). This mutant is incapable of pseudouridylating 16S and 23S rRNA and subsequently referred to as ΨΔrRNA strain. We cultured cells under optimal (37ºC) and various stress conditions (28ºC, ampicillin, gentamicin, and high NaCl) to maximize the detection of Ψ sites, especially in transcripts which are mostly expressed under sub-optimal conditions (Supplementary Fig. 1A). Using the pipeline described above, we initially sought to determine the distribution of all base deletions in our dataset. Because the GC content of E. coli is approximately 50% and assuming deletions occur randomly, the expected deletion counts for a given base, relative to other bases, is expected to be similar under non-treated conditions. In other words, the average deletion ratio for one base (e.g. U), relative to others should be ~ 1. (i.e. mean(U:A, U:G, U:C)). Indeed, in untreated samples, we observed the expected deletion frequencies for each base (Fig. 2A). In contrast, U-sites, but not other bases, were deleted > 5 times above the expected frequency in BS-treated samples. This is consistent with the hypothesis that BS-treatment results in deletions at Ψ positions13-15.

For downstream analysis, we included an additional filter to minimize the detection of false positives. We required a putative Ψ site to be identified in at least two samples ─ whether from biological replicates or different stress conditions. This resulted in 1,954 high-confidence sites in 1,331 transcripts, representing almost 30% of the transcriptome. Importantly, the deletion rates were highly reproducible across biological replicates, indicating a quantitative ability of BS-based methods to detect Ψ in bacterial mRNAs (Supplementary Fig. 1B, Supplementary Data 1). To further validate our approach, we searched for known Ψ sites in rRNA. Although all samples were subjected to ribodepletion during sample prep (Fig. 1), there is usually a small percentage (~ 0.5%) remaining which is sufficient for many analyses17. Our pipeline accurately identified the sole Ψ site in 16S rRNA, and 8 of 10 sites in 23S rRNA (Fig. 2B and C). In contrast, we did not detect any Ψ in rRNAs from ΨΔrRNA samples (Fig. 2B and C). Therefore, our approach accurately identified known Ψ positions in rRNA with no false positives detected.

Next, we sought to identify established Ψ sites in tRNAs. Pseudouridines are known to be deposited in position 13 (D-stem loop), positions 32, 38, 39, and 40 in the anticodon-stem loop, and positions 55 and 65 in the TΨC-loop6,18,19. Similarly to the rRNA data above, we accurately detected all known tRNA sites (Fig. 2D), despite differences in length due to the variable arm region. For instance, with our pipeline, we found position 72 in the TΨC loop of tRNA-Ser is pseudouridylated, which is equivalent to position 55 in the TΨC loop of tRNA-Val (Fig. 2E). Furthermore, we uncovered a previously unreported Ψ site in position 8 of tRNA-Tyr (Fig. 2E). This is unlikely to be a false positive since the sequence motif (GUUC) is similar to that found in TΨC loops13,15 (Fig. 2E).

Lastly, we obtained four different PUS knockout strains that are incapable of pseudouridylating specific positions in tRNAs and grown under optimal conditions (37ºC). We found ΔtruA strains were unable to deposit Ψ in the anticodon-stem loop (Fig. 2F). Similarly, we did not detect Ψs in TΨC-loop and D-stem loop of tRNAs in ΔtruB and ΔtruD strains, respectively (Fig. 2F, Supplementary Data 2). These observations are consistent with established target locations of TruA, TruB, and TruD PUS enzymes18,19. However, although TruC modifies position 65 in the TΨC-loop, we still observed Ψ at this location in ΔtruC mutant (Fig. 2F). One possible explanation could be that TruB may deposit Ψ in TruC sites since they recognize similar motif as described below. Thus, TruC and TruB may have redundant roles in this context. Altogether, our approach accurately identifies known Ψ positions in rRNA and tRNAs of E. coli.

A comprehensive Ψ landscape in E. coli transcriptome
In a previous study to map the repertoire of Ψ in E. coli RNAs, Schaening-Burgos et al identified Ψ in only 42 mRNAs11 using Pseudo-seq. However, using the highly sensitive BS-based approach, we detected high-confidence Ψ sites in 1,217 mRNAs, which is 27 times above previous estimates (Fig. 3A, Supplementary Data 1). Most mRNAs only harbor a single Ψ site (Supplementary Fig. 2A). In addition, we found few mRNAs that are always pseudouridylated irrespective of growth conditions. For example, the multidrug efflux transporter, mdtG, contained Ψs in both normal and all tested stress conditions (Fig. 3B). In addition, almost 20% of detected pseudouridylated mRNAs are either involved in the biosynthesis of secondary metabolites or adaptation to diverse environments (Fig. 3C). This suggests a possible regulatory role of Ψ in response to stress, similarly to eukaryotes5.

Next, we analyzed the motif frequency and distribution of all Ψ sites. The most frequent pseudouridylated motifs contained mostly G or C bases upstream of the Ψ site (Fig. 3D, Supplementary Data 3). Interestingly, the most abundant motif (GGUAU), which we detected in over 100 sites, is also reported as being a frequent Ψ motif in the human transcriptome15. To further understand the sequence preference for each PUS enzyme, we compared the Ψ proportion (i.e. deletion ratio) for every motif in WT and mutant strains. In other words, a motif with a significantly lower deletion rate in a PUS mutant is most likely a sequence preference. We established that target sites for rRNA PUS enzymes occur in a consensus sequence comprised of ‘UUGC’ (Fig. 3E) which agrees with previously reported RluA recognition motif11. Similarly, we observed ‘GUUC’ as the main recognition sequence for TruB, in strong agreement with the homologous TRUB1 targets in human and yeast13,20 (Fig. 3E). TruC, on the other hand, deposits Ψ mostly in ‘UCC’ sequences and possibly recognizes the ‘GUUC’ sequence similarly to TruB (Fig. 3E). However, we did not identify a consensus sequence for TruA or TruD. This may suggest a lack of sequence preference or pseudouridylation by other PUS enzymes.

We also examined the predicted secondary structure of all target sites using RNAfold21 and observed pseudouridylation occurs predominantly in unpaired uridine sites (Supplementary Fig. 2B). However, deletion rates were similar in both paired and unpaired sites (Supplementary Fig. 2C). Although, using synthetic probes, bisulfite treatment was previously shown to have a modest structural dependency due to both chemical accessibility of bisulfite15 and transcriptase preference13,15; however, all Ψ sites were identified albeit with varying deletion ratios suggesting secondary structures do not qualitatively impede Ψ detection. Overall, Ψs are widespread across E. coli mRNAs and are mostly deposited in loop or hairpin structures by rRNA and tRNA PUS enzymes.

Ψ stabilizes bacterial mRNA
Because Ψs are strongly linked to stability of non-coding RNA22,23, we examined their possible role in mRNA stability and gene expression. To begin, we retrieved transcripts containing Ψ that can be assigned unambiguously to a PUS enzyme. In this case, sites with significantly lower Ψ (i.e. deletion ratio) in a PUS mutant, relative to WT, were considered a recognition site for a given PUS. Since multiple PUS can deposit Ψ at different locations in a single transcript, we assigned those candidates to multiple PUS. Next, we determined the abundance levels of candidate mRNAs from read counts. We observed that mRNAs unable to be pseudouridylated in PUS mutants were significantly less abundant (Fig. 4A). This was true for all PUS targets, with the sole exception of TruA. In addition, we did not observe significant differences in mRNA abundance between WT and mutants when random control sets of equal size were compared which suggests a link between Ψ and abundance (Supplementary Fig. 3A).

To further understand how Ψ may influence abundance, we performed rifampicin treatment and RNA sequencing (RIF-seq)24 to measure mRNA decay rates (Fig. 4B). We included spike-in controls during sample preparation for normalization25 and measured RNA abundance at 0, 3-, 6-, 12-, and 24-min post transcription arrest with rifampicin. We calculated mRNA half-lives in both WT and PUS mutants and observed RNAs unambiguously modified by a PUS enzyme were significantly less stable in mutants except ΔtruA (Fig. 4C). This suggests that the decrease in mRNA abundance is a consequence of transcript instability. We also compared the global RNA half-lives of all pseudouridylated and non-pseudouridylated transcripts and observed no significant differences in either the WT or mutant strains. (Supplementary Fig. 3B). Taken together, our results underscore a functional role of pseudouridylation in stabilizing bacterial mRNA.

Quantitative mapping of Ψ landscape in human microbiome samples
To examine the ability of our approach to quantify pseudouridylation of mRNA in complex microbial communities, we recruited 18 subjects (13 healthy and 5 periodontitis patients) and retrieved oral plaque samples from the subgingival region (Fig. 5A). After BS treatment and sequencing, we mapped reads to over 5,000 genomes assembled from publicly available oral metagenomes26 (Fig. 5A). Unlike our approach for E. coli isolates where we required a Ψ site to have a deletion rate < 1% in untreated samples, we considered sites with higher deletions to accommodate indels arising from strain heterogeneity in complex microbial communities. However, we required the deletion rate in BS-treated samples to be greater than twofold in untreated samples. To further minimize false positives, we used a statistical approach that takes into consideration total read depth, deletion rate, and relationship between these parameters in BS-treated and untreated samples14.

Next, we analyzed the microbiome composition and abundance using Kraken 2/Bracken27. We initially determined, using publicly available paired metagenomics and metatranscriptomics oral samples, that Kraken 2/Bracken accurately recapitulates microbial taxonomic information from metatranscriptome data (Supplementary Fig. 4). In our patient samples, the microbial community was dominated by diverse microbial species belonging to multiple phyla in agreement with the known microbial profiles of subgingival microbiomes28,29 (Fig. 5B). In total, we identified 3,534 Ψ sites from 3,135 protein-coding transcripts, distributed across 218 species (Fig. 5C and D, Supplementary Data 4). Similarly to the E. coli data above, many of the annotated mRNAs harboring Ψ sites are either involved in the biosynthesis of antibiotics or other secondary metabolites, or adaptation to diverse environments (Fig. 5E). This suggests that a regulatory role of Ψ in response to stress may be widespread across bacteria.

We also looked at the sequence preference for Ψ deposition and identified 254 different motifs ─ many of which were previously identified in E. coli. (Fig. 5F, Fig. 3D). Although we were unable to assign most motifs to a PUS enzyme, the top two motifs (CCUCC and GCUCC) resemble the TruC-associated motif in E. coli (Fig. 5F, Fig. 3E).

Because most of the enriched motifs are GC-rich, we hypothesized that pseudouridylation will be more predominant in mRNAs from GC-rich genomes. To exclude any potential sequence bias arising from chemical accessibility of bisulfite or transcriptase preference, we analyzed the number of identified Ψ sites per genome rather than the degree of pseudouridylation per site (i.e. deletion ratio). As mentioned above, the motif sequence does not qualitatively impede Ψ detection14. For each genome, we performed normalization by taking into consideration the coverage depth (i.e number of bases above the depth cutoff) and genome length. Interestingly, we observed significantly more Ψ sites in protein-coding transcripts from GC-rich genomes (Fig. 5G).

We also examined the relationship between Ψ and mRNA abundance in the community. Similarly, we focused on the number of Ψ sites observed in a transcript. Strikingly, the number of modification sites in an mRNA was positively correlated with transcript abundance (Fig. 5H). This data raises the possibility that multiple Ψ sites might be required to enhance the stability of an mRNA, though further investigation is necessary to confirm this requirement.

Finally, we sought to identify new Ψ sites in rRNA from our microbial community. The 16S of E. coli, for example, harbors a single Ψ site and the repertoire of rRNA modification sites across multiple bacteria have not been investigated. We began by retrieving 16S sequences from our genomes and only considered high-confidence sites to minimize false positives (i.e. sites identified in multiple genomes or samples) (Supplementary Data 5). We found new Ψ locations in 13 different species belonging to 5 phyla (Fig. 5I). Unlike E. coli, the 16S of some identified microbes had multiple modification sites. For instance, we identified 6 and 4 Ψ sites in the 16S RNA from Rothia aeria and Corynebacterium matruchotii, respectively. On the other hand, Aggregatibacter segnis, which belongs to the same Enterobacterales order like E. coli, harbors a high confidence modification site (position 1054) in addition to the established E. coli site (position 512) (Fig. 5I, Fig. 2B). Modifications in similar sites (positions 1054 – 1057) were observed in other Proteobacteria species such as Eikenella corrodens, Citrobacter freundii, and Haemophilus pittmaniae.

Altogether, Ψ modifications are prevalent in microbial communities and our approach provided an avenue to study posttranscription regulation in microbiomes.

Methods
Bacterial strains and growth conditions
Unless stated otherwise, wild-type (MC415), ΨΔrRNA (MC452), ΔtruA, ΔtruB, ΔtruC, and ΔtruD mutants were routinely grown on Luria-Bertani (LB) media at 37ºC (Supplementary Data 6). When required, 25 μg/ml of kanamycin was added to the media of Δtru mutants. To induce stress conditions, bacteria (MC415 and MC452) were grown in the presence of ampicillin (0.4 μg/μl), gentamicin (0.4 μg/μl), or high salinity stress conditions (4% NaCl). In the latter, cells were cultured in Tryptic soy broth (TSB).

RNA extraction, bisulfite treatment and library preparation
All fresh cultures were grown until OD of 0.3 and subsequently preserved in RNA protect (Qiagen, #76506). Next, RNA was extracted from the cultures using RNeasy mini kit (Qiagen, #74104) and residual DNA was removed by DNA-freeTM DNA Removal Kit (Thermo Fisher Scientific, AM1906). After RNA quantification and quality analysis by Qubit Flex and Bioanalyzer, RNA samples were diluted with nuclease-free water to give approximately 400-700 ng in 22 μl of volume which was thereafter split in two halves (i.e ‘BS-treated’ and ‘untreated’ halves).

Ribosomal RNA was depleted in both halves using Ribo-Zero Plus rRNA Depletion Kit (Illumina, #20040525) followed by fragmentation of all samples by addition of 0.9 μl of fragmentation reagents (Thermo Fisher Scientific, #AM8740). The mix was incubated at 95°C for 20 s and fragmentation was stopped using 0.9 μl of stop reagent. Samples were immediately placed on ice.

Fresh bisulfite reagent (BSR) was prepared by adding 0.27 g of sodium sulfite (Sigma Aldrich, #901916) and 0.034 g of sodium bisulfite (Sigma Aldrich, 799394) to 900 μl of DEPC treated water. Next, 45 μl of freshly prepared BSR was added to the 11 μl of fragmented ‘BS-treated’ half and incubated at 70°C for 3 h. After incubation, 75 μl of nuclease-free water was added to the mix followed by 270 μl of RNA binding buffer (RNA Clean and Concentrator-5 column kit, Zymo Research, R1015) and 400 μl of 100% ethanol. The entire mixture (~800 μl) was loaded onto the column and centrifuged for 30 s. The column was then washed with 200 μl of RNA wash buffer and subjected to desulfonation by adding 200 μl of RNA desulphonation buffer (Zymo Research #R5001-3-40) to the column. This was subsequently incubated at room temperature for 75 min. The column was washed using RNA wash buffer and eluted in 10 μl of elution buffer.

For the other ‘untreated’ half, samples were purified after fragmentation and eluted with 10 μl of elution buffer. Both purified BS-treated and untreated RNA samples were then annealed to random hexamers (50 μM) and subjected to cDNA first strand synthesis using SuperScript™ IV reverse transcriptase (Thermo Fisher Scientific, #18090010). The reaction parameters were: 23°C for 10 min, 50°C for 1 h, and 80°C for 10 min. Next, the second strand of cDNA was synthesized by DNA polymerase 1 (Thermo Fisher, #EP0041) in a reaction mixture containing 0.8 μl of RNAse H (Thermo Fisher, #EN0202) followed by incubation at 15°C for 2 hr and thereafter 75°C for 10 min. This double stranded cDNA was purified by adding 90 μl of Ampure beads (Beckman Coulter, A63881) using exact instructions provided by Illumina (Illumina Stranded Total RNA Prep, Ligation with Ribo-Zero Plus reference guide). Subsequent steps for 3’ ends adenylation, anchor ligation and indexing were performed using the Illumina Stranded Total RNA Prep, Ligation with Ribo-Zero (Illumina, #20040525), anchor plate (Illumina, #20040899), and DNA/RNA UD indexes set (Illumina, #20091646), respectively. Generated libraries were evaluated using Bioanalyzer and sequenced on the NovaSeq platform.

Sample preparation for half-life estimation
Fresh bacterial cultures were grown in LB media to an OD of 0.3. Transcription initiation was halted using rifampicin (0.5 mg/ml) and samples were collected at 0, 3, 6, 12 and 24 min post rifampicin treatment. The reaction was immediately stopped using stop solution (95% ethanol + 5% Qiazol lysis solution (Qiagen, #79306)) followed by snap freezing in liquid nitrogen. The samples were allowed to thaw on ice, centrifuged, and bacterial pellets were resuspended in TE buffer containing 15 mg/ml lysozyme and 20 ul/ml proteinase K. To this mixture, 2 ul of a 1/10 of ERCC RNA spike-in sequences25 (Invitrogen, P#4456740) was added for downstream normalization. These spiked mixtures were further processed for RNA extraction and library preparation according to Illumina protocol.

Human oral sample collection
Collection of human samples was performed on an IRB clinical protocol approved at the National Institutes of Health (NIH) Clinical Center (ClinicalTrials.gov ID NCT01568697). This study included 18 participants: 13 healthy volunteers and 5 patients with chronic periodontitis. All study participants provided written informed consent for participation in this study. Participants were deemed systemically healthy based on detailed medical history and select laboratory work up. In addition to systemic screening, periodontal status was assessed through detailed clinical oral evaluation. Subgingival plaque samples (tooth adherent biofilm) were removed using a Gracey Curette (HuFriedy Group) from patients and samples were immediately placed in DNA/RNA Shield Stabilization Solution (Zymo Research, #R1100-50).

RNA extraction from human samples
Samples were vortexed to dissolve pellets and then added to 2ml ZR BashingBead Lysis Tubes (0.1 & 0.5 mm) (Zymo Research, #S6012-50) with 600 μl of Qiazol lysis solution (Qiagen, #79306). Bead beating was performed for 5 min in FastPrep 24 homogenizer at maximum speed. After centrifugation, 180 μl of chloroform was added to the supernatant and centrifuged for another 15 min at 4°C for phase separation. The aqueous phase was transferred to a new tube and 1.5:1 v/v 100% ethanol was added. The mixture was mixed by pipetting and RNA was extracted. Ribodepletion was performed on the extracted RNA using Illumina® Ribo-Zero Plus rRNA Microbiome depletion kit (Illumina, #20072062) and subsequent steps were processed as described in the previous section.

Bioinformatics pipeline to detect Ψ
Adapter removal and quality trimming of raw reads was performed using Trim Galore34. We mapped clean reads to E. coli BW25113 genome (NZ_CP009273.1) using bwa-mem (v0.7.17)35 and realigned reads using ABRA2 to improve detection of indels in downstream analysis16. Next, bam-readcount (v1.0.1) was used to retrieve nucleotide coverage and sequence variant information36 and the resulting output was subsequently parsed using brc-parser.py37. To be considered a Ψ site, a uridine position needs to meet the following criteria: (i) coverage depth ≥ 20 in both BS-treated and untreated samples; (ii) deletion count ≥ 5 in BS-treated samples; and (iii) deletion ratio ≥ 5% in BS-treated samples but less than 1% in untreated samples. Ψ was derived as the difference in deletion ratio between BS-treated and untreated sites and a Ψ site was considered as high confidence if it was identified in at least 2 samples in the main dataset (Supplementary Fig. 1A).

Motif analysis
To determine the sequence preference for each PUS, we retrieved flanking sequences of mRNA sites having deletion ratio > 6%, which is the median deletion ratio in our dataset. We further excluded low occurrence motifs which were identified less than 6 times. For every motif passing the threshold, deletion ratios were compared between WT and ΨΔrRNA strains, or between the four different tru mutants. A significantly decreased motif (P value < 0.05, two-sided Wilcoxon rank sum test) in a mutant was considered a sequence preference for that PUS. To identify consensus sequences, motifs were submitted to MEME and the parameter to search the reverse complement strand was disabled38.

Pathway analysis
The pathways of genes associated with pseudouridylation were analyzed using KEGG. The KEGG Mapper search tool was employed to map and assign these genes to the corresponding pathways39.

RNA structure analysis
tRNA structures were predicted with RNAfold web server using default parameters21. To investigate base-pairing preference for pseudouridylation, 12-mer sequences flanking each side of the Ψ residue were used as input in the command-line version of RNAfold (-p -T 37 parameters).

Estimation of mRNA abundance from RNA-seq and influence of Ψ on stability
Because untreated samples are equivalent to traditional RNA-seq, we quantified the mRNA abundance from read counts in untreated samples. To begin, the generated .bam files from previous mapping step were used as input in featureCounts (v2.1.1)40 to retrieve unique read counts mapped to each gene. We included a pseudocount of one to each gene and normalized using ‘transcript per million’ (TPM).

To evaluate the role of Ψ on mRNA abundance, we initially calculated the Ψ-strength for each transcript13. Ψ-strength is defined as the sum of deletion ratios at all Ψ sites within one RNA. The gene-level deletion ratios were compared between WT and PUS mutants using only samples obtained at 37ºC because tru mutants were solely grown at optimum temperature. For a gene to be considered pseudouridylated by a given PUS, the difference in gene-level Ψ-strength between WT and a mutant must be ≥ 0.05 in both biological replicates.

To measure mRNA half-lives, reads from samples collected at timepoints post rifampicin treatment were mapped to E. coli genome plus the ERCC spike-in RNA sequences25 and read counts per gene determined using featureCounts (v2.1.1)40. We used the 25 most abundant ERCC spike-in RNA for normalization and all timepoints were scaled relative to timepoint 0. Here, a normalization constant (k) was determined for each ERCC spike-in RNA using k = ERCCn / ERCC0 * Total0 / Totaln, where ERCCn and ERCC0 are RNA read counts at timepoints n and 0, respectively, while Totaln and Total0 are the total library size in samples Tn and T0 41. The geometric mean of all k values (at each timepoint), called normalization factor, was used to normalize the library. This was achieved by dividing raw read counts by the normalization factor. Transcripts per million (TPM) for each library were subsequently calculated. To account for the initial stability period observed post-rifampicin treatment, mRNA decay was modeled using a delayed first-order exponential function42. This two-phase model assumes a constant mRNA abundance during an initial lag phase, followed by exponential decay. Model fits for selected candidates were manually inspected to ensure data fidelity. Candidates with poor model fit were excluded.

To identify transcripts that are not PUS substrates (i.e. non-pseudouridylated), we filtered for those with, at least, moderate expression across all samples (TPM > 20) but no detectable Ψ in any sample. Requiring consistent expression ensured these high-confidence candidates were unlikely to be false negatives resulting from insufficient sequencing coverage in our pseudouridine detection pipeline.

Microbiome analysis
We downloaded 5,113 oral isolates and metagenome-assembled genomes (MAGs) from the Human Reference Oral Microbiome (HROM) database26. For each genome, we predicted and annotated genes using Prokka (v. 1.14.6)43. To reduce the computational load associated with analyses of a large dataset, we used a coverage cutoff to determine microbes for downstream analysis. For each subject, only genomes with > 20X coverage in both BS-treated and untreated samples were processed. Reads were first mapped to the human genome to exclude non-microbial reads. The filtered reads were subsequently mapped to HROM genomes using bwa-mem (v0.7.17)35, realigned using ABRA216, and per-base coverage was retrieved using bam-readcount (v1.0.1)36. Uridine positions passing the following filters were considered as pseudouridylated: (i) read depth ≥ 20 in both BS-treated and untreated samples; (ii) deletion count ≥ 5 in BS-treated samples; (iii) deletion ratio ≥ 0.02 (i.e. 2%) in BS-treated samples; (iv) deletion ratio in BS-treated samples is, at least, 2 fold above deletion ratio in untreated samples; and (v) the P-value from Fisher’s exact test should be < 0.01 when total number of deletions and reads depth in BS-treated and untreated samples are compared. For sites passing the above threshold, the proportion of Ψ was derived as the difference in deletion ratio between BS-treated and untreated sites (i.e. Δ deletion ratio).

Data regarding the GC content of each genome was derived from Cha et al. (2025)26. To estimate the number of pseudouridine sites per microbe, we initially determined the total number of base positions (i.e. A, T, G, or C) with depth > 20X in both BS-treated and untreated samples. This number was divided by 1,000 to derive bases per Kb. The total number of Ψ sites was subsequently divided by the derived bases/Kb factor.

For 16S analyses, we used strict threshold cutoff to minimize false positive discoveries. First, we only processed Ψ sites with Δ deletion ratio > 0.1 (i.e. 10%). Second, we required putative candidates to be present in at least 3 genomes or less than 3 genomes but identified in 3 or more samples. Finally, to evaluate the accuracy of Kraken 2/Bracken in metatranscriptomics data, we retrieved paired metagenomics and metatranscriptomics oral samples from Belstrøm et al44.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper presents a well-developed approach to studying pseudouridine (pseU) modifications in bacteria using a bisulfite (BS)-based method, providing valuable data on pseU sites in rRNA and mRNA. The authors effectively highlight the novelty of reporting the landscape of pseU sites in low-abundance RNA species, such as mRNA. The bisulfite reaction employed mirrors that used in PRAISE and BID-seq, and the application of this technology to bacterial samples, where polyA enrichment of mRNA is challenging, is commendable. This well-established protocol will undoubtedly benefit future research on bacterial pseU and other mRNA modifications.\nHowever, considering the recent publications in Nature Communications that highlight significant discoveries and comprehensive studies on the biological functions of RNA modifications in various systems, including mammals and plants, this paper is not suitable for the journal. It has not adequately addressed the actual function of pseudouridine (pseU) in regulating bacterial gene expression through a comprehensive study. Here are a few suggestions for improvement:\n1. In Fig. 2, the description of tRNA pseU could be more precise. The reverse transcriptase used in this study cannot read through heavily modified or methylated regions within bacterial tRNA. It would be beneficial to conduct an AlkB-based demethylation assay first to remove methylated bases, followed by library construction. This demethylation treatment could be applied immediately after RNA fragmentation, allowing the reverse transcriptase to read through all tRNA regions with accurate coverage, thereby ensuring reliable detection of pseU.\n2. While the results showing approximately 1,200 pseU sites on bacterial mRNA in Fig. 2 appear solid, additional evidence is needed to elucidate the functional roles of these pseU modifications. Conducting mRNA stability assays or translation assays (such as Ribo-seq) would help clarify how mRNA pseU regulates mRNA stability and translation. Additionally, gene knockdown or knockout of pseU writer proteins could serve as a valuable perturbation to assess changes in mRNA pseU levels. Site mutations at pseU-modified U sites could also be employed to validate the effects of pseU on mRNA stability and translation.\n3. In Fig. 6, considering the complexity and dynamic nature of mRNA expression in clinical samples, it would be prudent for the authors to monitor mRNA profiles before discussing pseU modifications. For example, if a pseU-modified mRNA from a specific bacterium uniquely decays in one patient while remaining intact in other patient samples, it would be essential to consider these mRNA expression dynamics in the discussion of pseU modification changes.\n4. There is a minor point regarding the frequently used term \"deletion fraction,\" which may cause confusion. Typically, the term \"deletion ratio\" is used to describe the observed reverse transcriptase deletion signature. Through calibration curves, these deletion ratios can then be converted into estimated pseU modification fractions. The pseU modification fraction is commonly referred to as pseU modification stoichiometry or pseU stoichiometry.", "reviewer_2": "The manuscript titled “Quantitative mapping of pseudouridines in bacterial RNA” by Sharma et al. describes the results of bisulfite-based deep sequencing in E.coli and oral microbiome samples to identify RNA pseudouridylation.\nThe manuscript presents a stronger section describing the methodology and overall results up to line 166 (Figures 1-3), but the subsequent discussion on potential Ψ-dependent RNA stabilization, novel transcripts, asRNAs, sRNAs, and the oral microbiome is substantially weaker.\nMany of the claims rely on assumptions, and alternative explanations have not been sufficiently explored. The manuscript would benefit from clearer, more rigorously supported arguments. For example, the assertion that “Pseudouridine mapping enabled the investigation of differential expression of very low abundance transcripts associated with the general stress response” lacks sufficient supporting data to be considered valid.\nIn conclusion, I do not recommend the manuscript for publication in its current form.\nMajor points:\n• Lines 86-91: Calculation of Ψ Proportion\no A primary concern is the lack of consideration for the effect of secondary structures. The authors should clarify whether bisulfite treatment has a structural dependency similar to SHAPE reagents. Later (lines 163-166), a weak structure dependency of pseudouridylation is reported. Do the authors have data to demonstrate whether this is due to an enzymatic preference of the PUSs or the chemical accessibility of bisulfite? This should be discussed in more detail.\no Even though the observed effect is relatively small, the background deletion rate should be subtracted before calculating the proportion. Additionally, data from bisulfite-treated PUS deletion strains should be considered as a background control to normalize for Ψ-independent effects of bisulfite treatment.\n• The statement “One possible explanation could be low efficiency of TruC in vivo since it is known to modify only two tRNAs” does not logically explain the observed phenomenon. It should be clarified whether the authors are proposing that TruC is not relevant in this context or that TruC and TruB may have redundant roles.\n• Figure 3E:\no Why is the wild-type deletion fraction not shown for GUUC and UCC? Based on the boxplot medians, it appears that the deletion rate in the ΔtruC strain is also lower, suggesting that the GUUC motif is used by both TruC and TruD. However, there is a higher variance. This point should be explicitly discussed.\n• Ψ Stabilization of Bacterial mRNA\no The observed reduction in normalized read counts of selected RNAs in PUS deletion strains needs stronger validation.\no A potential pleiotropic effect of PUS deletion should be ruled out. The authors should analyze random control sets of equal size for comparison.\no Replicates appear to be averaged; a proper differential expression analysis using DESeq or a similar tool would strengthen the claim.\no Even if the effect is statistically significant, a change in transcript stability is only one possible explanation. Indirect effects on RNA synthesis rates must be considered.\no To conclude that Ψ increases stability, global half-life measurements are necessary. Additionally, even if a change in half-life is observed in a specific RNA, it could be an indirect and Ψ-independent effect of PUS deletion. Depending on the individual RNA, stability changes could be also either positive or negative.\n• Ψ Profiling and Detection of Low Abundance Transcripts\no The claim that Ψ profiling identifies differentially expressed transcripts better than conventional RNA-seq relies on the assumption that pseudouridylation rates are dependent on cellular RNA concentrations. However, this assumption is not supported by data in the manuscript.\no Pseudouridylation may be influenced by co-factors, PUS concentrations, RNA-binding proteins, other RNAs, or other unknown factors.\no Among the 14 genes with differential pseudouridylation, only 7 showed differential expression in RT-qPCR, which challenges the underlying hypothesis.\no The manuscript should clearly distinguish between increased sensitivity for differential expression analysis and increased sensitivity for RNA detection in general.\no The detection of new asRNAs and sRNAs is not inherently dependent on the bisulfite method, rather it results from different experimental conditions. Also e.g. an increased sequencing depth alone could enhance detection sensitivity. This should be explicitly stated throughout the manuscript.\n• sRNA Functional Analysis\no The discussion of sRNA function is superficial and does not substantially contribute to the main findings.\no TargetRNA3 was trained on interactome data and may not be the most suitable tool for detecting functional sRNA interactions. Other tools, such as IntaRNA, might be more appropriate.\no To confirm that sucA is a genuine sRNA target, compensatory mutations in both the target and sRNA should be tested to rule out secondary effects.\n• Microbiome Analysis\no The microbiome section could serve as a proof of principle, but given that E. coli alone has 1,331 Ψ sites, the detection of only 174 sites across the microbiome suggests that a significant fraction remains undetected. This indicates that the sensitivity of the approach needs substantial improvement before strong conclusions can be drawn.\no The statement that TruB is the major PUS enzyme in the microbiome is weakly supported and should be rephrased. A more precise statement would be that a majority of the limited number of detected sites resemble the TruB-associated motif in E. coli.\no The authors should also consider that the motif might be recognized by other enzymes in different bacterial species and that TruB may have additional or alternative motifs in other bacteria. Additionally, the observed enrichment of TruB sites appears weak and should be discussed accordingly.\nMinor points:\n• The abstract prominently states that bisulfite sequencing eliminates the need for strand-specific RNA sequencing of asRNAs. However, this is not a major advantage. In contrast, the detection of low-expressed unmodified asRNAs will be hindered without stranded sequencing.\n• The sentence “Therefore, new approaches are required to evaluate the widespread distribution of Ψ in bacterial mRNAs” assumes a priori that Ψ is widely distributed. It should be revised to say “evaluate the assumed widespread distribution of Ψ” or “evaluate the widespread distribution of Ψ detected in this study.”\n• The sentence “We began by extracting total RNA and enriching for mRNA through ribodepletion” is imprecise. The method enriches for all RNA species except ribosomal RNA, not just mRNA.\n• Line 162: The structural analysis appears to be based on predictions. This and the used tool should be clearly stated in the text.\n• Line 199: The sentence “Altogether, Ψ profiling can facilitate the discovery of low-expressed genes associated with bacterial adaptation to varying environmental conditions which may be difficult to study using conventional RNA-seq analysis pipelines” should be refined. Even if the underlying hypothesis is correct, the claim should be limited to facilitating the discovery of differential expression, rather than gene discovery in general.\n• If a computational sRNA target prediction tool such as TargetRNA3 is used to generate data discussed in the main text, the tool should also be explicitly mentioned in the main text."}, "author_response": "We thank the reviewers for their insightful comments. This has significantly helped to improve our manuscript. Based on their comments, we restructured our manuscript to focus on the functional role of pseudouridines in bacterial mRNA. In this case, we excluded weaker sections such as those pertaining to sRNA analysis and detection of low abundance transcripts. In addition, we included 18 microbiome samples and new sequence data to measure mRNA global half-life in E. coli. Please find below our responses to the comments.\nReviewer 1\nThis paper presents a well-developed approach to studying pseudouridine (pseU) modifications in bacteria using a bisulfite (BS)-based method, providing valuable data on pseU sites in rRNA and mRNA. The authors effectively highlight the novelty of reporting the landscape of pseU sites in low-abundance RNA species, such as mRNA. The bisulfite reaction employed mirrors that used in PRAISE and BID-seq, and the application of this technology to bacterial samples, where polyA enrichment of mRNA is challenging, is commendable. This well-established protocol will undoubtedly benefit future research on bacterial pseU and other mRNA modifications.\nWe thank the reviewer for the positive feedback on our work\nHowever, considering the recent publications in Nature Communications that highlight significant discoveries and comprehensive studies on the biological functions of RNA modifications in various systems, including mammals and plants, this paper is not suitable for the journal. It has not adequately addressed the actual function of pseudouridine (pseU) in regulating bacterial gene expression through a comprehensive study. Here are a few suggestions for improvement:\n1. In Fig. 2, the description of tRNA pseU could be more precise. The reverse transcriptase used in this study cannot read through heavily modified or methylated regions within bacterial tRNA. It would be beneficial to conduct an AlkB-based demethylation assay first to remove methylated bases, followed by library construction. This demethylation treatment could be applied immediately after RNA fragmentation, allowing the reverse transcriptase to read through all tRNA regions with accurate coverage, thereby ensuring reliable detection of pseU.\nWe agree with the reviewer that a demethylation step may improve detection. However, there are already known pseudouridine sites in E. coli tRNAs. In our data, we identified all these sites. In addition, our tRNA data was simply for validation purposes rather than discovery of new tRNA sites.\n2. While the results showing approximately 1,200 pseU sites on bacterial mRNA in Fig. 2 appear solid, additional evidence is needed to elucidate the functional roles of these pseU modifications. Conducting mRNA stability assays or translation assays (such as Ribo-seq) would help clarify how mRNA pseU regulates mRNA stability and translation. Additionally, gene knockdown or knockout of pseU writer proteins could serve as a valuable perturbation to assess changes in mRNA pseU levels. Site mutations at pseU-modified U sites could also be employed to validate the effects of pseU on mRNA stability and translation.\nWe totally agree. We performed RIF-seq experiments to measure global E. coli mRNA half-lives in both wild-type and PUS mutants. We observed that pseudouridine significantly enhances mRNA stability (Fig. 4).\n3. In Fig. 6, considering the complexity and dynamic nature of mRNA expression in clinical samples, it would be prudent for the authors to monitor mRNA profiles before discussing pseU modifications. For example, if a pseU-modified mRNA from a specific bacterium uniquely decays in one patient while remaining intact in other patient samples, it would be essential to consider these mRNA expression dynamics in the discussion of pseU modification changes\nWe agree. Rather than using gene clusters (which are derived from multiple bacteria), we studied mRNA dynamics in individual microbes present in the community. Our new analysis is similar to that for E. coli isolates, but in this case, extended to multiple bacteria simultaneously (Fig. 5).\n4. There is a minor point regarding the frequently used term \"deletion fraction,\" which may cause confusion. Typically, the term \"deletion ratio\" is used to describe the observed reverse transcriptase deletion signature. Through calibration curves, these deletion ratios can then be converted into estimated pseU modification fractions. The pseU modification fraction is commonly referred to as pseU modification stoichiometry or pseU stoichiometry.\nWe agree and have implemented the change throughout the manuscript.\nReviewer 2\nThe manuscript presents a stronger section describing the methodology and overall results up to line 166 (Figures 1-3), but the subsequent discussion on potential Ψ-dependent RNA stabilization, novel transcripts, asRNAs, sRNAs, and the oral microbiome is substantially weaker.\nMany of the claims rely on assumptions, and alternative explanations have not been sufficiently explored. The manuscript would benefit from clearer, more rigorously supported arguments. For example, the assertion that “Pseudouridine mapping enabled the investigation of differential expression of very low abundance transcripts associated with the general stress response” lacks sufficient supporting data to be considered valid.\nWe agree with the reviewer. We have excluded weaker sections on sRNAs and low abundance transcripts and refocused our manuscript on the functional role of pseudouridines. In addition, we increased the number of microbiome samples analyzed from 2 to 18 which has substantially improved our analysis.\nMajor points:\n• Lines 86-91: Calculation of Ψ Proportion\nA primary concern is the lack of consideration for the effect of secondary structures. The authors should clarify whether bisulfite treatment has a structural dependency similar to SHAPE reagents. Later (lines 163-166), a weak structure dependency of pseudouridylation is reported. Do the authors have data to demonstrate whether this is due to an enzymatic preference of the PUSs or the chemical accessibility of bisulfite? This should be discussed in more detail.\nBisulfite-based protocols such as BID-seq1 and PRAISE2 (which we adopted) do not specifically report on structural dependency of bisulfite treatment. However, we analyzed their data and made the following observations. (i) In Zhang et al1., using synthetic probes containing different sequence motifs, bisulfite treatment was shown to have a very minor structural dependency. In one instance, deletion ratio for AGΨCU was 85% whereas AGΨUG was 50%. We suspect this is both due to chemical accessibility of bisulfite and transcriptase preference1,2. Nevertheless, pseudouridines were robustly identified albeit with varying deletion ratios suggesting secondary structures do not qualitatively impede Ψ detection. We have included these points in the main manuscript text (lines 158 – 165).\no Even though the observed effect is relatively small, the background deletion rate should be subtracted before calculating the proportion. Additionally, data from bisulfite-treated PUS deletion strains should be considered as a background control to normalize for Ψ-independent effects of bisulfite treatment.\nWe agree. We have now subtracted background deletion rates prior to calculating deletion ratios. In addition, it is difficult to use PUS deletion strains to normalize for Ψ-independent effects because some sites are modified by more than one PUS enzyme (Lines 129 – 130).\nThe statement “One possible explanation could be low efficiency of TruC in vivo since it is known to modify only two tRNAs” does not logically explain the observed phenomenon. It should be clarified whether the authors are proposing that TruC is not relevant in this context or that TruC and TruB may have redundant roles.\nWe apologise for the lack of clarity. We have modified the sentence to suggest that TruC and TruB may have redundant roles since they recognize similar motifs (Lines 130 – 133).\nFigure 3E:\no Why is the wild-type deletion fraction not shown for GUUC and UCC? Based on the boxplot medians, it appears that the deletion rate in the ΔtruC strain is also lower, suggesting that the GUUC motif is used by both TruC and TruD[B]. However, there is a higher variance. This point should be explicitly discussed.\nWe have included boxplots for the WT (Fig. 3E). We have also included the suggestion on the possibility of TruC and TruB recognizing the GUUC motif (Lines 154 -156).\nΨ Stabilization of Bacterial mRNA\no The observed reduction in normalized read counts of selected RNAs in PUS deletion strains needs stronger validation. A potential pleiotropic effect of PUS deletion should be ruled out. The authors should analyze random control sets of equal size for comparison.\nWe agree. We have analyzed random control sets of equal size which further corroborates a link between pseudouridines and mRNA abundance (Supplementary Fig. 3).\no Replicates appear to be averaged; a proper differential expression analysis using DESeq or a similar tool would strengthen the claim.\nWith DESeq or similar tools, it requires, at least, triplicate samples for robust analysis. In our case, we only had two biological replicates.\no Even if the effect is statistically significant, a change in transcript stability is only one possible explanation. Indirect effects on RNA synthesis rates must be considered.\nWe agree. We further performed RIF-seq to measure global RNA half-lives in WT and PUS mutants. Our data indicate pseudouridine stabilizes E. coli mRNA (Fig. 4).\no To conclude that Ψ increases stability, global half-life measurements are necessary. Additionally, even if a change in half-life is observed in a specific RNA, it could be an indirect and Ψ-independent effect of PUS deletion. Depending on the individual RNA, stability changes could be also either positive or negative.\nWe agree. Please see response above\n• Ψ Profiling and Detection of Low Abundance Transcripts\no The claim that Ψ profiling identifies differentially expressed transcripts better than conventional RNA-seq relies on the assumption that pseudouridylation rates are dependent on cellular RNA concentrations. However, this assumption is not supported by data in the manuscript.\no Pseudouridylation may be influenced by co-factors, PUS concentrations, RNA-binding proteins, other RNAs, or other unknown factors.\no Among the 14 genes with differential pseudouridylation, only 7 showed differential expression in RT-qPCR, which challenges the underlying hypothesis.\no The manuscript should clearly distinguish between increased sensitivity for differential expression analysis and increased sensitivity for RNA detection in general.\no The detection of new asRNAs and sRNAs is not inherently dependent on the bisulfite method, rather it results from different experimental conditions. Also e.g. an increased sequencing depth alone could enhance detection sensitivity. This should be explicitly stated throughout the manuscript.\nWe totally agree. As mentioned above, we have excluded the section on detection sensitivity of low abundance from our manuscript.\n• sRNA Functional Analysis\no The discussion of sRNA function is superficial and does not substantially contribute to the main findings.\no TargetRNA3 was trained on interactome data and may not be the most suitable tool for detecting functional sRNA interactions. Other tools, such as IntaRNA, might be more appropriate.\no To confirm that sucA is a genuine sRNA target, compensatory mutations in both the target and sRNA should be tested to rule out secondary effects.\nWe agree with the reviewer. We have excluded this section as mentioned above.\n• Microbiome Analysis\no The microbiome section could serve as a proof of principle, but given that E. coli alone has 1,331 Ψ sites, the detection of only 174 sites across the microbiome suggests that a significant fraction remains undetected. This indicates that the sensitivity of the approach needs substantial improvement before strong conclusions can be drawn.\nWe agree. To improve this section, we analyzed 18 microbiome samples (compared to 2 in our initial submission) and utilized a more efficient microbiome ribodepletion kit to enrich for non-rRNA. We identified 3,534 Ψ sites from 3,135 protein-coding transcripts, distributed across 218 species in the microbiome (Fig. 5).\no The statement that TruB is the major PUS enzyme in the microbiome is weakly supported and should be rephrased. A more precise statement would be that a majority of the limited number of detected sites resemble the TruB-associated motif in E. coli.\no The authors should also consider that the motif might be recognized by other enzymes in different bacterial species and that TruB may have additional or alternative motifs in other bacteria. Additionally, the observed enrichment of TruB sites appears weak and should be discussed accordingly.\nWe agree and have modified the text to exclude claims of TruB being the major PUS enzyme (Line 211).\nMinor points:\n• The abstract prominently states that bisulfite sequencing eliminates the need for strand-specific RNA sequencing of asRNAs. However, this is not a major advantage. In contrast, the detection of low-expressed unmodified asRNAs will be hindered without stranded sequencing.\nWe have excluded the section on asRNA from the manuscript.\n• The sentence “Therefore, new approaches are required to evaluate the widespread distribution of Ψ in bacterial mRNAs” assumes a priori that Ψ is widely distributed. It should be revised to say “evaluate the assumed widespread distribution of Ψ” or “evaluate the widespread distribution of Ψ detected in this study.”\nWe agree and have modified the text appropriately (Lines 54 -55).\n• The sentence “We began by extracting total RNA and enriching for mRNA through ribodepletion” is imprecise. The method enriches for all RNA species except ribosomal RNA, not just mRNA.\nWe agree and have modified the sentence (Line 79).\n• Line 162: The structural analysis appears to be based on predictions. This and the used tool should be clearly stated in the text.\nThis has been included in the main text (Line 158).\n• Line 199: The sentence “Altogether, Ψ profiling can facilitate the discovery of low-expressed genes associated with bacterial adaptation to varying environmental conditions which may be difficult to study using conventional RNA-seq analysis pipelines” should be refined. Even if the underlying hypothesis is correct, the claim should be limited to facilitating the discovery of differential expression, rather than gene discovery in general.\nAs mentioned above, we have excluded this section from the manuscript.\n• If a computational sRNA target prediction tool such as TargetRNA3 is used to generate data discussed in the main text, the tool should also be explicitly mentioned in the main text.\nSimilarly, we have excluded this section from the manuscript.\nReferences\n1. Zhang, L.S. et al. BID-seq for transcriptome-wide quantitative sequencing of mRNA pseudouridine at base resolution. Nat. Protoc. 19, 517-538 (2024).\n2. Zhang, M. et al.Quantitative profiling of pseudouridylation landscape in the human transcriptome. Nat. Chem. Biol. 19, 1185-1195 (2023)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The manuscript has been significantly improved, and I support its publication in Nature Communications. The authors should further review the main text, legends, and data figures to ensure accuracy and correctness, after AIP.", "reviewer_2": "The authors addressed the points mentioned in the reviewer comments and significantly improved the manuscript.\nHowever, the required and important inclusion of the global RNA stability data raised a new issue which needs to be addressed prior to acceptance.\nThe method to calculate the RNA stability based on reference 41 seems not to be the state of the art. It was shown that rifampicin decay curves show a characteristic “delay” of the exponential decline of the RNA dependent on the distance to the transcriptional start site. Not considering this delay leads to a false half-life which is often overestimated (PMID: 25583150, PMID: 37454177).\nPlease also compare the global RNA stabilities (all transcripts without pseudouridine) vs pseudouridine RNAs in WT and the mutant strains to look for pleiotropic patterns. How do the non-modified RNAs react in the mutants?\nAssuming steady state, what proportion of the foldchange in the transcript abundance is explained by the stability differences for your selected candidates?\nMinor observation:\nGroups in Fig 5h seem to have unequal variances and sample sizes which violates the assumptions of the Wilcoxon test. The Brunner-Munzel might be more appropriate. The largely different sample size also hampers the general interpretations about stability and abundance of multiple modified RNAs."}, "author_response": "Response to reviewers\nWe thank the reviewers once again for their positive comments. Please find below our responses to the remaining comments.\nReviewer 1\nThe manuscript has been significantly improved, and I support its publication in Nature Communications. The authors should further review the main text, legends, and data figures to ensure accuracy and correctness, after AIP.\nWe thank the reviewer for the positive feedback.\nReviewer 2\nThe method to calculate the RNA stability based on reference 41 seems not to be the state of the art. It was shown that rifampicin decay curves show a characteristic “delay” of the exponential decline of the RNA dependent on the distance to the transcriptional start site. Not considering this delay leads to a false half-life which is often overestimated (PMID: 25583150, PMID: 37454177).\nWe thank the reviewer for pointing us in this direction. We have now accounted for the initial stability period observed post-rifampicin treatment by modelling mRNA decay using a delayed first-order exponential function. This two-phase model assumes a constant mRNA abundance during an initial lag phase, followed by exponential decay. The new results are similar to previous ones, and the conclusions remain the same (Fig. 4C).\nPlease also compare the global RNA stabilities (all transcripts without pseudouridine) vs pseudouridine RNAs in WT and the mutant strains to look for pleiotropic patterns. How do the non-modified RNAs react in the mutants?\nWe compared the global RNA half-lives of all pseudouridylated and non-pseudouridylated transcripts and observed no significant differences in either the WT or mutant strains. (Supplementary Fig. 3B).\nAssuming steady state, what proportion of the foldchange in the transcript abundance is explained by the stability differences for your selected candidates?\nWe are careful not to assign a specific numerical value to the contribution of half-lives to abundance, as such an estimate would likely be inaccurate given the complexity of bacterial gene expression. For instance, accelerated degradation may feedback to alter the transcription rate for specific genes.\nMinor observation:\nGroups in Fig 5h seem to have unequal variances and sample sizes which violates the assumptions of the Wilcoxon test. The Brunner-Munzel might be more appropriate. The largely different sample size also hampers the general interpretations about stability and abundance of multiple modified RNAs.\nWe agree. We performed statistical analysis using Brunner-Munzel test and still observed significant differences between the groups (similarly to Wilcoxon test results) (Fig. 5H). Nevertheless, we have modified our interpretation to make it cautionary (Line 230)."}]

Decision letter