source: https://doi.org/10.48350/152107 | downloaded: 31.1.2022
Article
Transcriptome-Wide Analysis of Stationary Phase Small ncRNAs in E. coli
Nicole Raad1,2,†, Hannes Luidalepp1,†, Michel Fasnacht1,2and Norbert Polacek1,*
Citation: Raad, N.; Luidalepp, H.;
Fasnacht, M.; Polacek, N.
Transcriptome-Wide Analysis of Stationary Phase Small ncRNAs in E. coli.Int. J. Mol. Sci.2021,22, 1703.
https://doi.org/10.3390/ijms22041703
Academic Editor: Margaret L. Britz Received: 16 December 2020 Accepted: 5 February 2021 Published: 8 February 2021
Publisher’s Note:MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affil- iations.
Copyright: © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).
1 Department of Chemistry, Biochemistry and Pharmaceutical Sciences, University of Bern, Freiestrasse 3, 3012 Bern, Switzerland; nicole.raad@dcb.unibe.ch (N.R.); hannes.luidalepp@ut.ee (H.L.);
michel.fasnacht@dcb.unibe.ch (M.F.)
2 Graduate School for Cellular and Biomedical Sciences, University of Bern, 3012 Bern, Switzerland
* Correspondence: norbert.polacek@dcb.unibe.ch
† These authors contributed equally to this work.
Abstract: Almost two-thirds of the microbiome’s biomass has been predicted to be in a non- proliferating, and thus dormant, growth state. It is assumed that dormancy goes hand in hand with global downregulation of gene expression. However, it remains largely unknown how bacteria manage to establish this resting phenotype at the molecular level. Recently small non-protein-coding RNAs (sRNAs or ncRNAs) have been suggested to be involved in establishing the non-proliferating state in bacteria. Here, we have deep sequenced the small transcriptome ofEscherichia coliin the expo- nential and stationary phases and analyzed the resulting reads by a novel biocomputational pipeline STARPA (Stable RNA Processing Product Analyzer). Our analysis reveals over 12,000 small tran- scripts enriched during both growth stages. Differential expression analysis reveals distinct sRNAs enriched in the stationary phase that originate from various genomic regions, including transfer RNA (tRNA) fragments. Furthermore, expression profiling by Northern blot and RT-qPCR analyses confirms the growth phase-dependent expression of several enriched sRNAs. Our study adds to the existing repertoire of bacterial sRNAs and suggests a role for some of these small molecules in establishing and maintaining stationary phase as well as the bacterial stress response. Functional characterization of these detected sRNAs has the potential of unraveling novel regulatory networks central for stationary phase biology.
Keywords:stationary phase; sRNA; tRF; dormancy; deep sequencing; bioinformatics
1. Introduction
Most of the bacterial population resides in a non-proliferating growth state termed the stationary phase. One of the most common environmental signals that prompt bacteria into adopting a non-proliferative behavior is nutrient deprivation [1]. In fact, environ- mental conditions that would allow constant cell division and bacterial growth are rarely encountered in nature [2]. Stationary phase bacterial cultures are not homogenous but represent a mixture of actively dividing and dying cells [1]. It is also during this phase that a minor population of dormant and non-dividing bacteria, termed persisters, arises [3–5].
Even though persister cells are already detectable in exponential phase cultures, their abundance increases upon entry into the stationary phase [1]. Due to their metabolic inactivity, persister cells exhibit an antibiotic tolerance phenotype [6,7] that poses a severe health threat and has been linked to recurrent bacterial infections. Pathogenic bacteria within the mammalian host often resort to dormancy to bypass immune response and cope with nutritive stress [5]. Therefore, entry into and exit from the stationary phase are crucial for the long-term survival of a bacterial population and are, therefore, tightly regulated processes. A key player in this crucial transition in Gram-negative bacteria is the alternative sigma factor RpoS, whose expression peaks in the stationary phase [8]. RpoS orchestrates the global downregulation of gene expression in favor of the upregulation of a subset of
Int. J. Mol. Sci.2021,22, 1703. https://doi.org/10.3390/ijms22041703 https://www.mdpi.com/journal/ijms
genes critical for survival under these challenging conditions. Furthermore, RpoS has been reported to aid in persister cell formation [9,10] and appears to promote dormancy hand in hand with (p)ppGpp [11], the alarmone regulator of the stringent response in bacteria [12].
Interestingly, RpoS expression is under the control of several sRNAs, both positively and negatively [13]. As such, the bacterial stress response and stationary phase behavior appear to be regulated in a hierarchal fashion, often relying on small ncRNA regulators.
In bacteria, the discovery of sRNA regulators has preceded that of their counterparts in eukaryotes, and these ribo-regulators can be divided into several categories based on their nature, their various modes of action, and the physiological processes that they regulate [14–16]. Transcriptome profiles unraveled a whole depository of uncategorized sRNAs [17,18], whose expression levels are heavily affected by key environmental signals.
These profiles are significantly dominated by short, trans-encoded sRNAs that regulate gene expression at the post-transcriptional level. A single sRNA can regulate the fate of several target mRNAs via partial base pairing, thereby fine-tuning pathways involved in metabolism [19], quorum sensing [20], antibiotic resistance [21], stress response [22], or virulence [23]. For instance,ryhBappears to aid persister cell formation in uropathogenic E. coli[24], while this sRNA has been previously reported to regulate iron metabolism [25].
This further illustrates the intricate nature of RNA regulatory networks in bacteria. These regulatory sRNAs heavily rely on the presence of RNA-binding proteins that protect them against degradation and/or assist them during mRNA target recognition. Hfq, and recently ProQ and CsrA, have been extensively characterized as key RNA chaperones binding, assisting, and stabilizing sRNAs [26–28]. As such, the identification and characterization of novel sRNAs frequently utilized these chaperones as baits in pull-down studies [29–31].
Nonetheless, several sRNA could be acting independently of these so far characterized chaperones and would thus be undetected in these pull-down approaches.
A so far largely unnoticed class of regulatory RNAs that has recently emerged in transcriptome studies is rancRNAs (ribosome-associated noncoding RNAs). rancRNAs have been described in all domains of life, and a subset of experimentally characterized candidates bind the ribosome, are thereby protected from cellular degradation machiner- ies, and regulate translation under specific stress conditions [32]. In eukaryotic systems, rancRNAs were shown to be capable of either inhibiting or stimulating protein synthesis by various mechanisms [33–35]. rancRNAs have also been described in prokarya; in the archaeonHaloferax volcanii, two rancRNAs were shown to inhibit translation either globally by competing with mRNA binding [36] or locally by inhibiting the translation of specific mRNAs [37]. However, the contribution of sRNA in general, and rancRNAs in particular, to stationary phase biology has not been investigated in a systematic and comparable manner transcriptome-wide. Therefore, our knowledge about these ribo-regulators in the bacterial domain remains fractional.
As bacterial sRNAs are often expressed under specific physiological conditions in response to external stimuli and can modulate the levels of several transcripts at once, these small molecules can globally orchestrate the aftermath of the physiological condition under which they are produced. Our understanding of the magnitude of sRNA involvement in weaving regulatory networks in dormant bacteria is still elusive and far from being under- stood in molecular terms. In the present study, the small transcriptome (20–300 nucleotides) of the Gram-negative bacteriumEscherichia coli(E. coli) was extracted from the total, as well as ribosome-associated RNA pools, from bacteria growing in the exponential and sta- tionary phases and different growth media. Subsequently, these sRNomes were converted into cDNAs, deep sequenced, and finally bioinformatically analyzed and compared. We developed and reported a new pipeline called STARPA (Stable RNA Processing Product Analyzer) to identify, quantify, and characterize putative stable RNA processing products or novel sRNAs from multiple RNA-seq libraries. Dual analysis of libraries reveals various upregulated sRNAs in the stationary phase, offering a global repertoire of sRNAs expressed under this growth condition and thus, representing prime candidates for being involved in establishing the non-proliferating state in the stationary phase. To validate the STARPA
pipeline, we experimentally verified several predicted RNA processing products derived from various RNA types inE. coli.
2. Results
2.1. Library Preparation and Sample Collection
In this current study, we aimed at identifying stably expressed sRNAs specifically enriched during the stationary phase of growth, which would potentially regulate the entry and maintenance of this state inE. coli. To reveal and identify these potential candidates, and in the absence of adequate time-resolved sRNA expression profiles recorded over different growth phases, we deep sequenced the smallE. colitranscriptome in exponential and stationary phase deriving from total or ribosome-associated RNA pools (Figure1a, Supplementary Table S1).E. colicells were grown separately in two media: standard LB (lysogeny broth) and MOPS-Glc (morpholinepropanesulfonic acid supplemented with glucose), the latter serving as a minimal medium to sustain bacterial growth. Samples for RNA extraction were taken from two distinct time points for each growth medium and in two biological replicates. The sRNome was extracted from the total as well as from ribosome-associated RNA pools from every sample collection time point per medium of growth and reverse transcribed into cDNA, leading to eight distinct libraries per RNA source (total or ribosome-associated), 16 in total (Figure1a, Supplementary Table S1).
Figure 1. Sequencing library preparation and analysis workflow. (a) Schematic representation of sequencing library preparation:E. coliwere cultured in two different media, lysogeny broth (LB) and morpholinepropanesulfonic acid (MOPS), and in two biological replicates. Samples were taken at the indicated time points (OD600 = optical density measured at 600 nm, h = hour). From each growth condition, samples were used for total RNA extraction and crude ribosome extraction, followed by RNA extraction. RNA was size extracted from polyacrylamide gels with 20–300 nt cutoffs (nt = nucleotides).
Size selected RNA was Tobacco Acid Pyrophosphatase (TAP) treated then reverse-transcribed with 50- and 30- adapters to make cDNA libraries. cDNA libraries were sequenced. (b) Workflow of Stable RNA Processing Product Analyzer (STARPA) analysis: Reads were trimmed from adapter sequences then aligned to theE. colireference genome. Unaligned reads were discarded in a sorting step. Reads with paired-ends were transformed into pseudo-single-end (SE) reads. The remaining reads were identified, clustered by sequence or by overlap, and then quantified (RPM = reads per million).
2.2. STARPA Workflow
We developed a novel and universal pipeline, which we denoted STARPA, for the analysis of the sequencing data. STARPA fulfills requirements absent from existing tools, such as the capability to handle multiple cDNA libraries in a comparative way and the
ability to identify potential RNA processing products that are supported by full length reads. STARPA was designed in a modular setup, allowing the application for various types of RNA sequencing data. The pipeline consists of seven steps (Figure1b). Briefly, raw demultiplexed reads were first trimmed to remove adapters and low-quality nucleotides.
Subsequently, trimmed reads were aligned to theE. colireference genome. After alignment, reads were filtered to remove unaligned, mismatched, and nonspecific reads. Next, PE (paired-end) reads were converted to SE (single-end) reads to allow subsequent analysis by Flaimapper [38], which identifies processing products. To reduce redundancy, the identified RNA processing products were clustered by overlap and by genomic context.
Following clustering, quantification was carried out and reported in a normalized manner as read per million of mapped reads (RPM) for the identified RNA processing products (Supplementary Tables S1 and S2). Aside from libraries of ribosome-associated RNA under MOPS conditions, all biological replicates displayed positive correlations (Supplementary Figure S1a,b).
2.3. Overview of the Small E. coli Transcriptome
Our STARPA analysis identified 22,057 reads in total from all growth conditions and RNA pools. A read was kept if it possessed at least 10 reads per million in a single library, and as such, 12,933 distinct reads were retained after this initial filtering step (Supplementary Table S1). Libraries were then grouped by origin (total RNA and ribosome- associated RNA libraries), and reads within the new grouping were further discarded if they did not fulfill the previous filtering criterion. Interestingly, 1755 transcripts were captured by both total and ribosome-associated RNA libraries (Figure2a, Supplementary Table S2). Furthermore, the majority of the sequenced transcripts from both libraries originated from intergenic or coding regions (Figure2b), and the overall distribution of represented biotypes was very similar among the two conditions.
To predict sRNAs involved in stationary phase biology, we conducted a differential expression analysis based on the captured reads per million from STARPA. Of relevant interest were candidates upregulated in stationary phase libraries. Differential sRNA expression analysis revealed a proportion of candidates upregulated in the stationary phase (Table1, Supplementary Table S3). The majority of the sRNome captured in our sequencing experiments remained unchanged in both the exponential and stationary phase (Table1, Supplementary Figure S2a–d), and a larger proportion of the differentially expressed transcripts were downregulated in the stationary phase. A portion of the 1267 predicted upregulated transcripts in the stationary phase overlapped between libraries (Figure2c). As such, STARPA predicted 660 sRNAs to be upregulated in the stationary phase. Reads from similar libraries displayed segregated clustering (Supplementary Figure S2e,f). Interestingly, the biotypes of the upregulated reads had different proportions than those in the total sRNA transcriptome (Supplementary Figure S3), with the majority derived from rRNA and tRNA regions in total RNA libraries and from coding regions in ribosome associated RNA libraries.
Figure 2.STARPA revealed over 12,000 sRNAs. (a) Venn diagram showing the number of identified processing products in total and ribosome-associated RNA libraries, with a cutoff threshold of 10 reads per million or higher. (b) Pie charts showing the percentage of identified RNAs by biotype for total and ribosome-associated RNA libraries. (c) Venn diagram showing the number of upregulated processing products in the stationary phase, in total, and ribosome-associated RNA libraries deriving from cells grown in LB and MOPS.
Table 1.Differential expression analysis. Table showing the number of downregulated, unchanged (not significant), and upregulated sRNAs in total RNA libraries and ribosome-associated RNA libraries for both lysogeny broth (LB) and morpholinepropanesulfonic acid (MOPS) growth media.
All numbers were based on differential expression analysis in the stationary phase in comparison to the exponential phase.
Description LB Total MOPS Total LB Ribosome MOPS Ribosome
Downregulated 1024 1098 323 299
Not significant 1671 1574 1797 1909
Upregulated 393 416 273 185
2.4. Validation of Selected Candidates
Deep sequencing of the small transcriptome ofE. coliidentified several small ncRNA candidates that were differentially expressed during the stationary phase compared to exponentially growing cells. Highly expressed candidates from the stationary phase
samples (Supplementary Table S3) whose sequences mapped with defined read stacks on the reference genome were selected for subsequent experimental validation. Northern blot analyses confirmed the abundant expression of some candidates in the stationary phase (Figure3a), belonging to distinct classes of identified transcripts. Interestingly, some of the top enriched stationary phase-specific sRNAs were tRNA-derived RNA fragments (tRFs).
Northern blot analyses confirmed the enrichment of these transcripts in the stationary phase (Figures3a and4a).
Figure 3.Validation of identified sRNAs. (a) Northern blot analysis of three candidates:sRNA_35(FM_+NC_000913.2_
4541696_34),23S rRNAfragment (FM_+NC_000913.2_3942854_40), and 50-tRFtrpT(FM_+NC_000913.2_3944979_33). 5S rRNAwas used as a loading control. Zero point four and 2 were optical density measurements, h = hours, d = days.
tRF = tRNA derived RNA fragment. Predicted sizes are displayed on the right in nt (nucleotides). (b) North- ern blot analysis of candidates in A. P100 = pellet enriched for ribosome-associated RNA. S100 = supernatant con- taining non-ribosome associated RNA. (c) RT-qPCR mean and SEM for sRNA_dps (FM_-NC_000913.2_0848114_53), sRNA_bolA (FM_+NC_000913.2_0453657_58), sRNA_osmE (FM_-NC_000913.2_1820262_45), and sRNA_yhfg (FM_- NC_000913.2_3489623_47) in the stationary phase (two biological replicates each). log2fold change was based on comparison with exponential phase samples for both conditions: Total S = Total RNA samples in the stationary phase. P100 S = ribosome- associated RNA in the stationary phase.
Figure 4. Stationary phase-specific transfer RNA fragments. (a) Northern blot analysis probing for30-tRFLeuU (FM_- NC_000913.2_3320093_38),50-tRFGly(FM_+NC_000913.2_4390382_32), and50-tRFAlaV(FM_+NC_000913.2_0225499_28).5S rRNAwas used as a loading control. Zero point four and 2 were optical density measurements, h: hours, d: days. tRF = tRNA derived RNA fragment. (b) Northern blot analysis probing for30-tRFLeuU,50-tRFGly, and50-tRFAlaV. P100 = pellet enriched for ribosome-associated RNA. S100 = supernatant containing non-ribosome associated RNA.
To validate enrichment of some transcripts in ribosome-associated libraries, northern blot analyses were carried out on RNA samples extracted from ribosome enriched pellets (P100) after ultracentrifugation and corresponding supernatant fractions (S100). As ex- pected from the differential sRNA expression analyses,sRNA_35andtRFtrpTwere enriched in non-ribosomal fractions, whereas a23S rRNAfragment was abundantly retained in ribosome fractions in a growth phase-dependent manner. Overall, tRFs appeared to be predominately non-ribosome bound inE. coli(Figures3b and4b).
Our sequencing results showed that a large subset of differentially expressed sRNA derive from mRNA coding regions (Figure2b, Supplementary Figure S3). To test if these mRNA-derived fragments were real sRNA candidates and did not represent mere cDNA library artifacts, quantitative real-time PCR analyses were performed with size-extracted RNA. This analysis showed that fragments deriving fromdps, bolA, osmE,andyhfgcoding regions were enriched in stationary phase in total RNA samples (Figure3c). All four predicted sRNAs were differentially expressed in ribosome-associated libraries according to our analysis (Supplementary Table S3), but onlysRNA_osmEappeared to be slightly enriched in P100 samples in our RT-qPCR analysis (Figure3c). In summary, experimental expression profile analyses of selected sRNA candidates captured during exponential growth or in stationary phase largely confirm the differential sRNA expression seen in our sequencing library analyses, thus highlighting the applicability of the novel STARPA biocomputational tool.
We next wondered if these stationary phase-specific sRNAs had physiological roles.
We pickedsRNA_35because of its high abundance (Figure3a, Supplementary Table S3) for our preliminary analysis. To evaluate the potential importance ofsRNA_35in the stationary phase, we overexpressed the sRNA in the exponential phase of growth, when the sRNA were typically absent. For this, the sequence of the sRNA was introduced into pBbE6k plasmid (Figure5a), and the sRNA were overexpressed in the exponential phase along
a control spacer sequence. Overexpression ofsRNA_35in the exponential phase was far more potent than its canonical expression in the stationary phase (Figure5b), showing that the overexpression experiment was successful.sRNA_35overexpression in the exponential phase led to the upregulation of the RNA polymerase alternative sigma factor rpoS to similar levels typically only observed in the stationary phase (Figure5c). Overexpression of a spacer sequence as control did not influencerpoSmRNA levels, thus highlighting the specificity of the observed effects. These pilot experiments suggest that sRNAs expressed in the stationary phase and predicted by our biocomputational analysis likely contain meaningful candidates in the context of stationary phase biology and thus can possess roles in bacterial dormancy that are yet to be deciphered.
Figure 5.Overexpression ofsRNA_35upregulatesrpoSexpression. (a) Schematic representation of constructed plasmids forsRNA_35(top) and control (bottom) overexpression. Plac= lac operon promoter, inducible by IPTG. rrnB T1 and T7TE are terminator sequences. (b) Northern blot analysis ofsRNA_35(intergenic sRNA). E = exponential phase, S = stationary phase, OE = overexpression. Predicted sizes are displayed on the right in nt (nucleotides). (c) RT-qPCR mean and SEM for rpoSexpression under different conditions (two biological replicates each). E = exponential phase, S = stationary phase, OE
= overexpression. log2fold change is reported for all conditions in comparison to exponential phase condition.
3. Discussion
Previous whole transcriptome high-throughput sequencing studies gave rise to a plethora of new data and enabled mapping of several previously uncharacterized RNA species. Several of these newly identified transcripts have been shown to play regulatory roles governing gene expression, with versatile mechanisms affecting post-transcriptional and translational events in all three domains of life. In bacteria, the ever-increasing reper- toire of sRNAs positions these short transcripts in complex networks regulating virtually all bacterial processes, including stationary phase, nutritional stress, and virulence [39]. Our deep sequencing analysis of the small transcriptome ofE. coliin exponential and stationary phases added to this repertoire of bacterial sRNAs, identifying several small ncRNAs differentially expressed during the stationary phase. These stationary phase-specific RNA molecules identified herein are prime candidates for so far unknown sRNAs potentially governing stationary phase regulatory networks. Based on these findings, subsequent dedicated experimental work is required to unravel whether or not these sRNA candidates indeed represent functional ribo-regulators of stationary phaseE. coliand/or persister cell formation.
Prediction of processing products: Several existing biocomputational tools allow the prediction and identification of RNA processing products from sequencing data. Our previ- ously established APART pipeline [40] allows the identification of RNA processing products
through the detection of sharp increases and decreases in read coverages (Figure6a). This identification strategy might lead to the prediction of RNA species that are not supported by existing reads and overlook existing processing products. STARPA uses Flaimapper [38]
for processing product annotation, which overcomes the mentioned limitations. However, in Flaimapper, peak detection on 50and 30end densities allows the mapping of processing products (Figure6b). Nevertheless, by this approach, overlapping processing products can be overlooked. To overcome this limitation, STARPA segregates reads by size first and then employs Flaimapper for peak identification (Figure6c). Furthermore, because Flaimapper is only compatible with SE (single-end) reads, STARPA converts PE (paired-end) reads to pseudo-SE reads (Figure1b). This approach allows STARPA to be compatible with both SE and PE reads, reducing the complexity of the pipeline by employing the same downstream workflow. Furthermore, within the quantification step of STARPA, several statistics can be collected for predicted sRNAs, such as relative coverage, coverage at single position level, consensus sequence, quality of consensus sequence, and uniqueness (Supplementary Table S1).
Figure 6.Comparison of biocomputational methods for processing product identification. (a) Scheme showing processing product identification by APART [40]. The processing products were identified by sites of consequence, sharp, and reversed coverage shift (green circles). Other sites of sharp coverage shift were discarded (blue circles). Some processing products were accepted (green checkmarks), and others were discarded (red checkmarks). #Reads = number of reads. (b) Scheme showing processing product identification by Flaimapper [38]. The processing products were identified by peak detection on 50and 30end densities, followed by reconstruction of fragments. Red cross = discarded peak. #Ends = number of ends.
(c) Scheme showing detection of overlapping processing products by STARPA. Overlapping processing products (red checkmark) could be discarded. To avoid this, STARPA fractionates reads by size first, then employs Flaimapper to identify processing products from each fraction.
3.1. sRNAs in Stationary Phase Biology
Our work highlighted several sRNAs as potential regulators of stationary phase biology (Figure 2c). The number of differentially expressed sRNAs was higher than our initial prediction. In fact, several sRNAs candidates predicted to be deriving from coding regions were not detected in Nnorthern blot analyses (data not shown). Using a more stringent cutoff level for read counts may better recapitulate stationary phase enriched sRNAs. Within the top upregulated sRNAs, several candidates were validated by Nnorthern blot analyses (Figures3a and4a), suggesting that our biocomputational analysis captured the status quo of these transcripts in the stationary phase. By lowering the cutoff levels, we captured 1199 differentially expressed sRNAs in the stationary phase (Supplementary Table S4), of which the mRNA-derived sRNA_yhfkwas confirmed by
Northern blot analysis (Supplementary Figure S4). However,sRNA_yhfkwas lost from our analysis with the currently used more stringent cutoff. One potential solution to this problem is to add an extra step in STARPA, whereby reads are further segregated into new categories based on their biotype and then analyzed separately. For instance, reads deriving from rRNA regions are much more abundant than those derived from coding regions, and as such, a higher threshold for the former can be employed, and lower thresholds can be used for the latter.
Predicted sRNAs differed between libraries with some overlap (Figure2c), suggesting that these upregulated sRNAs belong to different classes of regulatory sRNA and may thus employ distinct mechanisms of action to establish or aid entry into the stationary phase. Our ribosome-associated sRNA libraries reported enrichment ofssrARNA that acted as tmRNA to release stalled ribosomes, a process more likely to occur upon nutritional starvation in the stationary phase [41]. As such, the upregulation ofssrAin the stationary phase is crucial for bacterial survival under challenging environmental conditions. These same sRNA libraries predicted the enrichment of a16S rRNAfragment in the stationary phase (FM_- NC_000913.2_3426703_81, Supplementary Table S3). We have previously demonstrated that this rRNA fragment is the result of16S rRNAcleavage specific to the stationary phase and that the 30S small ribosome subunits carrying the fragmented16S rRNA translate less efficiently than their exponential phase-derived counterparts [42]. Consequently, this reported processing product could represent a signal of ribosome shutdown in the stationary phase. Similarly, our analysis reported an abundant23S rRNAfragment in the stationary phase (Figure3a), whose function remains unknown, along with other rRNA- derived and ribosome-associated sRNAs predicted by STARPA (Figure2c, Supplementary Table S3). These sRNAs might be signals of processing events in broader physiological responses and not regulatory sRNAs in themselves, such as the characterized16S rRNA fragment. However, several of these processing products enriched in ribosome-associated sRNA libraries did not derive from abundant rRNA transcripts (Supplementary Figure S3). This suggests that these molecules were protected from degradation by binding the ribosome [32] and may thus act as rancRNAs. While rancRNAs are poorly characterized in bacteria, very recently, an antitoxin RNA has been identified as rancRNA inStaphylococcus aureusthat downregulates global protein synthesis by affecting tRNA binding, which, in turn, promotes persister cell formation [43]. As such, our enriched ribosome-associated sRNAs could play similar roles in dimming translation, a critical process for establishing and maintaining bacterial stationary phase.
Our analysis also reported previously characterized sRNAs with important roles in the stationary phase, particularly enriched in total RNA libraries. For instance, and not so surprisingly, our differential sRNA analysis positioned6S RNAas one of the most enriched sRNAs in the stationary phase (Supplementary Table S3). 6S RNA directly associates with RNA polymerase and its sigma factor 70, σ70, following its accumulation in the stationary phase and guides the downregulation of severalσ70-dependent promoters [44], a process coupled to the upregulation of RpoS-dependent transcripts. Furthermore, 6S RNA regulates the expression of several genes involved in metabolism [45], positioning this sRNA as a key regulator in establishing and maintaining bacterial dormancy.sibCis a cis-acting sRNA encoded antisense to its target mRNAibsCcoding for a toxin protein that is primarily involved in growth arrest [46]. The downregulation of theibsC-encoded toxin expression in the stationary phase is coupled to the upregulation of the sRNA [47]. This growth-phase dependent expression ofsibCremains poorly understood.ChiXis a trans- acting sRNA that strictly relies on Hfq to downregulate the expression of its targets in the stationary phase [48,49]. As our analysis captured several sRNAs with diverse regulatory mechanisms, we also predicted that the so-far uncharacterized portion of this sRNome could potentially be functional in the stationary phase. Out of the sRNA candidates that we validated experimentally, we also noted that the abundantsRNA_35(Figure3a) appeared to play a role in the stationary phase. Its overexpression is coupled to the upregulation of rpoS(Figure6c), a process innate to stationary phase [8] when it regulates transcription of
about 500 stationary phase genes [50]. RpoS is not only known as a master regulator of general stress response in bacteria but its upregulation has also been connected to elevated persister cell formation [9,10]. While it is not clear if our observed up-regulation ofrpoS was direct or indirect,sRNA_35appeared to be a functional stationary phase-specific sRNA with a role yet to be unraveled. More dedicated experimental work is required to gain functional insight into the roles of this and other sRNAs identified in this study.
3.2. tRF Enrichment in Stationary Phase
Interestingly, several tRNA-derived RNA fragments (tRFs) were enriched in the stationary phase (Figures3a and 4a). Recently, tRFs have gained vast recognition as regulatory molecules in all three domains of life [51]. tRFs appear to be involved in a wide range of biological roles, including the regulation of transcription, translation, stress granule formation, apoptosis, cell proliferation, tumor suppression, RNAi, vesicle-mediated intercellular communication, intergenerational inheritance, retrotransposons mobility, cell differentiation, and ribosome biogenesis [52–55]. What kind of functions can these RNA molecules serve in simpler forms of life, such as bacteria? Bacterial tRFs, while not yet characterized to an equally detailed level compared to other domains, appear to be involved in a more complex level of gene expression regulation. A 30-ETS (external transcribed spacer) of a pre-tRNA inE. coliassociates with the RNA chaperone Hfq, however, not to regulate gene expression as most Hfq-associated sRNAs do, but to sequester sRNAs away from their mRNA targets, thereby preventing transcriptional noise [56]. The same study suggests several other potential tRFs that are likely functional based on the conservation of tRNA ETS and ITS (internal transcribed spacer) sequences among enterobacteria. Our validated sequencing data for several tRFs inE. coliconfirmed the stationary phase-specific enrichment of these processing products. However, the functional scope of these tRF candidates likely surpassed that of serving as decoys as suggested earlier. Indeed, in the archeonHaloferax volcanii, a 50-tRF was produced under alkaline stress conditions and downregulates gene expression on a global level by binding the small ribosomal subunit, thus interfering with mRNA loading [36]. Most of these tRFs were enriched in total RNA libraries, and expression profiling shows that these molecules were likely not binding the ribosome directly inE. coli(Figures3b and4b). Unlike its peers,tRFAlaVwas enriched in ribosome-associated RNA libraries (Supplementary Table S3), but strikingly, Nnorthern blot analyses showed that this fragment was not interacting with ribosomes.
One explanation for this discrepancy between sequencing data and experiments may be that tRFAlaVwas interacting with other RNA molecules or protein partners closely associating with the ribosome, and as such, was temporally captured in our ribosome-associated sRNA libraries. In fact, the majority of the differentially expressed tRFs in our analysis may be interacting with RNA partners or RNA binding proteins; for instance,tRFGly(Figure4a) was enriched in ProQ RIL-seq (RNA interaction by ligation and sequencing) and deletion experiments [57]. This suggests that this tRF is interacting with the chaperone ProQ and other partners, likely RNA molecules. Overall, the stationary phase-specific expression behavior of these tRNA fragments argues for putative physiological roles, and further functional analysis needs to be conducted to understand if these molecules are serving as molecular decoys or trans-encoded sRNAs.
In this study, we surveyed the small transcriptomes ofE. coliin both the exponential and stationary phases by deep sequencing and showed that the composition of these sRNomes differed between both growth stages. Our refined biocomputational approach STARPA allowed the straightforward identification of several stable sRNAs that were likely aiding in the establishment of the resting phenotype or even orchestrating its establishment as part of regulatory networks. We suggest several sRNAs for further functional characteri- zation, including abundant uncharacterized tRF species whose roles remain enigmatic in the bacterial kingdom. To uncover the function of these tRFs (and other sRNA candidates identified in our stationary phase sRNome study), a combination of molecular, biochemical, and genetic approaches is required.
4. Materials and Methods 4.1. Strains and Media
E. colistrain MG1655 (F−λ−ilvG rfb-50 rph-1) [58] carrying pETgfp-mut2-AGGAGG(3) plasmid [59] was used. Bacteria were grown in LB-Miller media and in MOPS (mor- pholinepropanesulfonic acid) media supplemented with 0.1% glucose (MOPS-Glc), as described [60]. All growth media were supplemented with kanamycin to a final concen- tration of 25µg/mL. For overexpression experiments, MG1655 bacteria with and without pBbE6k plasmids [61] were used. Media was supplemented with IPTG (Isopropylβ-d-1- thiogalactopyranoside) (Roth. Germany) to a final concentration of 1 mM.
4.2. Cell Growth and Sample Collection
For sequencing and subsequent validation, 10µL of bacterial DMSO (dimethyl sulfox- ide) stock [60] was used to inoculate 2 mL of MOPS that was then incubated for 24 h at 37◦C. Two hundred microliters or 1 mL of this pre-culture were used to inoculate 200 mL media (LB and MOPS-Glc, respectively). The cultures were incubated at 37◦C, and cells were collected at desired time points (exponential phase: OD600 0.4, stationary phase:
20, 28, 40, and 48 h for LB, 2 or 5 days for MOPS) by pelleting for 5 min at 20◦C (for volumes≤200 mL) or for 10 min (for volumes > 200 mL) at 11,000×gat 4◦C. Cells were washed with 1-2 mL 1×PBS (Phosphate-buffered saline, 137 mM NaCl, 2.7 mM KCl, 10 mM Na2HPO4, 1.8 mM KH2PO4) and pelleted again by centrifugation for 5 min at 8000×gat 4◦C (Eppendorf 5804R Centrifuge, Hamburg, Germany). The pelleted cells were frozen in liquid nitrogen and stored at−80◦C for later RNA extraction.
For overexpression experiments, MG1655 bacteria carrying plasmids pBbE6k-35mer, pBbE6k-empty, or nothing, were grown overnight at 37◦C. Fifty-millimeter cultures were started from overnight cultures and grown to an OD600= 0.4. Induction with IPTG was then done. Samples for RNA extraction were collected before induction and 2 h after induction.
4.3. Construction of sRNA Expression Plasmids
sRNA_35overexpression plasmid, pBbE6k-sRNA_35, was constructed by insertion of thesRNA_35sequence at the transcriptional +1 site underPlacO control in pBbE6k- RFP [61] by MEGAWHOP cloning [62]. Briefly, primers NR0001 and NR0002 were used to create a mega primer spanning the genomicsRNA_35sequence by PCR using Phusion DNA Polymerase (NEB, Ipswich, MA, USA). The resulting PCR product was purified with Wizard® SV Gel and PCR Clean-Up (Promega, Madison, WI, USA). The purified mega primer was used in excess to insertsRNA_35into the pBbE6k-RFP plasmid and replace the RFP sequence with PCR. The PCR reaction was digested with DpnI (NEB, Ipswich, MA, USA) for 2 h at 37◦C and transformed into MG1655E. colicompetent cells. Plasmids were extracted from positive clones by Wizard® Plus SV Minipreps DNA Purification System (Promega, Madison, WI, USA) and sequenced with HL0199 primer (Microsynth, Balgach, Switzerland). pBbE6k-empty was constructed by one-step cloning PCR [63]
with primers NR0008 and NR0009 and pBbE6k-RFP plasmid. PCR reaction was DpnI treated as previously mentioned, and subsequent transformation, plasmid extraction, and sequencing were carried out, as described for pBbE6k-sRNA_35. The primers used for these manipulations are listed in Supplementary Table S5.
4.4. Preparation of Crude Ribosomes
The crude ribosomes were prepared according to standard protocols [34,42]. Briefly, total cell lysates were prepared by homogenization in FastPrep®-24 (MP Biomedicals, Illkirch Cedex, France) using 0.1 mm beads in lysis buffer (20 mM Tris-Cl pH 7.5, 100 mM NH4Cl, 10 mM MgCl2, 0.5 mM EDTA, and 6 mMβ-mercaptoethanol). Lysates were cleared from beads and cell debris by several rounds of centrifugation, and crude ribosomes were extracted following ultracentrifugation (Beckman Coulter Optima-XPN-80 Ultracentrifuge, Brea, CA, USA) at 100,000×g(P100), and supernatants were retained (S100).
4.5. cDNA Library Preparation and Deep Sequencing
cDNA libraries were prepared as described [34,42]. Briefly, 500 ng of size extracted RNA (18-300 nucleotides) from all listed experimental conditions and total RNA or ribosome-associated RNA pools were treated with TAP (Tobacco Acid Pyrophosphatase) (Epicenter/Lucigen, WI, USA) to remove pyrophosphate from the 50-end. The TruSeq Small RNA Library Prep kit (Illumina, San Diego, CA, USA) was used to prepare cDNA libraries, according to the manufacturer’s guidelines. A unique index primer was used for each RNA preparation, and two independent cDNA libraries were generated for each studied condition, all of which were subjected to paired-end deep sequencing analyses on an Illumina HiSeq platform (University of Bern, Bern, Switzerland).
4.6. RNA Extraction and Northern Blot Analysis
RNA extraction was performed as described [42]. Two to fifteen micrograms of total RNA or RNA isolated from crude ribosomes were separated on 8% polyacrylamide gels (7M Urea, 1×TBE), and gels were run for 2 h at 200 V. RNA was transferred to a nylon mem- brane (Amersham Hybond-N+, GE Healthcare, Chicago, IL, USA) using a semi-dry blotter (V20-SDB, Scie-Plas) and crosslinked to membranes using a microprocessor-controlled UV irradiation system (BLX-254, Vilber Lourmat, Witec AG, Switzerland). Hybridization was performed as described [64]. The sequences of end-labeled DNA oligonucleotides with [γ-32P]-ATP and used for hybridization are featured in Supplementary Table S5.
4.7. Quantitative Real-Time RT-PCR (qPCR)
One microgram of total RNA or 100 ng of size-extracted RNA (20-150 nt) from all experimental conditions were treated with DNase I (Thermo Scientific, Waltham, MA, USA) to digest any leftover DNA, according to the manufacturer’s protocol. Samples were then reverse transcribed into cDNA with SuperScript™ IV One-Step RT-PCR System (Invitrogen, Carlsbad, CA, USA) and random primer hexamers (Thermo Scientific, Waltham, MA, USA).
cDNA samples were treated with RNase H (NEB, Ipswich, MA, USA) to hydrolyze leftover RNA. qPCR was done using GoTaq®qPCR Master Mix (Promega, Madison, WI, USA), 50-fold diluted cDNA, and 500 nM to 1 uM of each primer. qPCR reactions were prepared by the CAS-1200 Corbett robot (Corbett Robotics, San Francisco, CA, USA) and were carried out using the Rotor Gene 6000, with suggested standard cycling conditions for gyrA and rpoS, and FAST cycling conditions for sRNA validation (Promega, Madison, WI, USA).5S andrecAwere used as an internal control for the normalization of gene expression. The samples were run in duplicates. The 2−∆∆CTmethod was used to calculate the fold-change relative to the control [65]. The mean log2fold-change and standard error of the mean were computed. Oligonucleotides used for RT-qPCR are provided in Supplementary Table S5.
4.8. Bioinformatics Analysis
EdgeR was used for differential gene expression analysis [66]. Default parameters of the STARPA pipeline were used. Instructions on how to obtain STARPA are located at https://github.com/luidale/starpa(accessed on 3 December 2020).
STARPA algorithm consisted of seven sequential tasks: trim, align, sam_sort, pseu- doSE, identify, cluster, and quantify. Each task can be run separately, or multiple tasks can be run sequentially. STARPA is compatible with both single-end (SE), or paired-end (PE) reads. Color space data is not supported. Each task requires specific input data, which is generally prepared by the preceding task. Task trim (read cleaning and first step of pipeline) supports data in FastQ format. FastQ files can be compressed as “.gz”, “.bz2” or
“.xz”. Reads must be previously demultiplexed.
4.8.1. Trim—Cleaning of Reads
Cutadapt [67] was used first to trim low-quality positions from the 30end (default cutoff phred quality score was set to 30) followed by adapter trimming, while also setting a minimum length (default: 18) for the processed reads. Default parameters of Cutadapt
were used unless otherwise stated. To determine the exact ends of the processing products or novel ncRNAs (to be predicted) and to increase the alignment level to the genome, the adapter removal must be maximized. In the case of PE (paired end) data, they were achieved by using 1 base minimum adapter overlap with the 30 end of the read. As a result, the tool sometimes removed some bases from the RNA insert sequence. In the case of overlapping reads, this loss will be covered by the second read from the pair when the PE reads are converted to SE read (step pseudoSE). When paired reads were not overlapping, only the length of the read was affected, while the ends of the RNA insert were not influenced. In the case of SE reads, the minimum adapter overlap with the 30end of the read was set to at least 3 bases. In addition, if an SE read was not adapter trimmed, it was rejected as it was not possible to determine the 30end of the RNA (RNA insert was longer than read). All rejected reads (too short for PE and SE, or untrimmed for SE) were saved in separate files to allow further assessment of quality issues of the libraries.
4.8.2. Align—Alignment to the Genome
Bowtie2 [68] was used to align reads to the reference genome with default parameters unless otherwise stated. To ease the downstream analysis of aligned reads, indels were avoided (--rdg 100,3, --rfg 100,3). As the pipeline was designed to also identify RNA processing products with potential modifications, Bowtie2 was adjusted to be more sen- sitive (allowing mismatches in seed sequence (-N 1)) and lowering minimum alignment score (--score-min L,0,−2). In addition, Bowtie2 was set to report all valid alignments (-a) reporting, thus all alignments for multi-mapped reads.
4.8.3. Sam_sort—Sorting of Aligned Reads
During sorting, alignments with the best alignment score (in case of PE alignments, alignment scores of paired alignments are summed) only were reported (best “stratum”) as SAM formatted output file. Unaligned reads were saved in separate files to allow further assessment of quality issues of the libraries. By default, steps align, and sam_sort were in sensitive mode. In this case, reads were initially aligned with default seed length (-L 22), then sorted and unaligned reads were aligned also using shorter seed length (-L 14), which was followed by sorting. In the end, aligned reads from both sorting steps were combined.
In sensitive mode, both align and sam_sort steps had to be run in combination and could not be run separately.
4.8.4. Pseudo SE—SE to PE Conversion
PE alignments were merged into pseudo SE alignments. Merging was conducted in overlapping regions by selecting the base from the alignment with higher sequencing quality, and combined sequencing quality was converted as described [69]. In case the paired alignments had no overlap, the gap between the aligned reads was filled with a genomic sequence, and a maximum sequencing quality score was given. Pseudo SE or SE (if input data for step pseudoSE is in a SE format) alignments with too many mismatches and reads with too many genomic alignments were discarded and saved in separate files to allow further assessment of quality issues of the libraries. All other reads received an NH tag (if not present) describing the number of reported alignments and were reported in SAM format. As poly- or oligoadenylation in bacteria can be relatively common on some RNA species [70,71], 30-oligo(A) tails were not considered as mismatches by default, and such reads were retained. This behavior can be switched off in the configuration file.
Reads with 30oligo(A) were also saved in separate files to allow further analysis of oligo or polyadenylation.
4.8.5. Identify—Identification of Processing Products
Flaimapper (version 3.0.0+) [38] was used to predict stable RNA processing prod- ucts. Flaimapper also predicts overlapping stable RNA processing products and is able to detect fragments sharing one end location if the length difference is more than 15
nucleotides. Next, the processing products were filtered by the minimum number of reads corresponding to any of the ends exceeding a set threshold (half the quantification threshold). The filtered predicted processing products were quantified more precisely via featureCounts [72] using parameters –s 1 (stranded), -M (allowing multi mapping reads), and -O (allowing multiple matches with predicted processing products). To take into consideration some positional inaccuracies in transcription initiation and RNA processing, parameters --nonOverlap and --nonOverlapFeature were used to set a number of allowed non-overlapping bases between reads and predicted processing products (default: 2).
4.8.6. Cluster—Clustering of Processing Products
Quantified input processing products were filtered by the read counts. Additional filtrating was done by relative coverage (average coverage of reads assigned to processing product divided by average coverage of all reads aligned to the positions of processing product) to remove candidates with a high background. Next, the processing products from all libraries analyzed were combined (identifying unique species) and clustered. Clustering was conducted as a two-step process.
Clustering by Overlap
As the prediction of processing products by Flaimapper is probabilistic, the predicted ends of the processing products in different libraries might vary slightly and might differ from the true ends of the fragments. As such, unique candidates with a set amount of non- overlapping positions were clustered by overlap, and representative processing products for individual clusters were selected. Representatives were selected in a repeating manner as a cluster member represents (set number of non-overlapping positions) the highest number of yet non-represented cluster members.
Clustering by Sequence
As most genomes contain repeating regions (repeat regions, rRNA operons, some tRNA genes, etc.), reads can be mapped to multiple positions. To reduce the number of identical candidates, clustering by sequence identity via CDI-HIT-EST [73] (100% identity, set by parameter -c 1) was employed. Because the genomic matches of these reads can be located in genomic regions with different surrounding sequences/contexts (e.g., different genes), clustering solely based on sequence identity can result in a loss of information. To avoid this, unique candidates which clustered by sequence identity had to be supported by the clustering of the contigs they overlap with (again via CDI-HIT-EST with identity threshold (-c) 90% and length difference cutoff (-s) 50%) and representative candidates for the clusters were selected.
4.8.7. Quantify—Quantification of Processing Products
Representative processing products were quantified and reported by parsing SAM files for every library. For multi mapped reads, each mapping gave 1 count. Additional characteristics (relative coverage, coverage at single position level, consensus sequence, quality of consensus sequence, genomic sequence, uniqueness) for each candidate were gathered and reported library wise to allow further assessment. Quantification was also reported in a normalized manner as read per million of mapped reads (RPM) and RPM of biotype (rRNA, tRNA, etc.) and RPM of biotype groups. In addition, general statistics from diverse steps of STARPA were collected and reported in a single file.
Supplementary Materials:The following are available online athttps://www.mdpi.com/1422-006 7/22/4/1703/s1.
Author Contributions: Conceptualization, N.P. and H.L.; experimental investigation, N.R. and H.L., software, H.L.; validation, N.R. and M.F.; writing-original draft preparation, N.R.; Funding acquisition and supervision, N.P. All authors have read and agreed to the published version of the manuscript.
Funding:This research was funded by the Swiss National Science Foundation grants 310030_188969 and 310030_197515 and by the Sciex grant 12.219 from the CRUS (Conférence des Recteurs des Universités Suisses).
Data Availability Statement:Sequencing data were deposited on the GEO database with the follow- ing accession number: GSE161907.
Acknowledgments:We thank Stefan Berger for experimental help during the initial phase of this project and Reidar Anderson and Maido Remm (University of Tartu) for their bioinformatics support and comments.
Conflicts of Interest:The authors declare no conflict of interest.
Abbreviations
cDNA complementary DNA CPM counts per million
d day
ETS external transcribed spacer FC fold change
h hour
ITS internal transcribed spacer mm millimeter
mM millimolar ncRNA noncoding RNA
nm nanometers
nt nucleotide
OD600 optical density at wavelength 600nm sRNA small RNA
tRF tRNA-derived RNA fragment
References
1. Roostalu, J.; Jõers, A.; Luidalepp, H.; Kaldalu, N.; Tenson, T. Cell division in Escherichia colicultures monitored at single cell resolution.BMC Microbiol.2008,8, 68. [CrossRef]
2. Jaishankar, J.; Srivastava, P. Molecular Basis of Stationary Phase Survival and Applications. Front. Microbiol. 2017, 8, 2000.
[CrossRef]
3. Jõers, A.; Liske, E.; Tenson, T. Dividing subpopulation of Escherichia coli in stationary phase.Res. Microbiol.2020,171, 153–157.
[CrossRef]
4. Levin, B.R.; Rozen, D.E. Non-inherited antibiotic resistance.Nat. Rev. Microbiol.2006,4, 556–562. [CrossRef]
5. Lewis, K. Persister cells, dormancy and infectious disease.Nat. Rev. Microbiol.2007,5, 48–56. [CrossRef]
6. Joers, A.; Tenson, T. Growth resumption from stationary phase reveals memory in Escherichia coli cultures.Sci. Rep.2016,6, 24055. [CrossRef] [PubMed]
7. Gray, D.A.; Dugar, G.; Gamba, P.; Strahl, H.; Jonker, M.J.; Hamoen, L.W. Extreme slow growth as alternative strategy to survive deep starvation in bacteria.Nat. Commun.2019,10, 890. [CrossRef] [PubMed]
8. Dong, T.; Joyce, C.; Schellhorn, H.E. The Role of RpoS in Bacterial Adaptation. InBacterial Physiology: A Molecular Approach;
El-Sharoud, W., Ed.; Springer: Berlin/Heidelberg, Germany, 2008; pp. 313–337.
9. Tkachenko, A.G.; Kashevarova, N.M.; Karavaeva, E.A.; Shumkov, M.S. Putrescine controls the formation of Escherichia coli persister cells tolerant to aminoglycoside netilmicin.FEMS Microbiol. Lett.2014,361, 25–33. [CrossRef] [PubMed]
10. Tkachenko, A.G.; Kashevarova, N.M.; Tyuleneva, E.A.; Shumkov, M.S. Stationary-phase genes upregulated by polyamines are responsible for the formation of Escherichia coli persister cells tolerant to netilmicin.FEMS Microbiol. Lett.2017,364. [CrossRef]
[PubMed]
11. Liu, S.; Wu, N.; Zhang, S.; Yuan, Y.; Zhang, W.; Zhang, Y. Variable Persister Gene Interactions with (p) ppGpp for Persister Formation in Escherichia coli.Front. Microbiol.2017,8, 1795. [CrossRef]
12. Traxler, M.F.; Summers, S.M.; Nguyen, H.T.; Zacharia, V.M.; Hightower, G.A.; Smith, J.T.; Conway, T. The global, ppGpp-mediated stringent response to amino acid starvation in Escherichia coli.Mol. Microbiol.2008,68, 1128–1148. [CrossRef] [PubMed]
13. Battesti, A.; Majdalani, N.; Gottesman, S. The RpoS-mediated general stress response in Escherichia coli.Annu. Rev. Microbiol.
2011,65, 189–213. [CrossRef] [PubMed]
14. Wagner, E.G.H.; Romby, P. Small RNAs in bacteria and archaea: Who they are, what they do, and how they do it.Adv. Genet.
2015,90, 133–208. [CrossRef] [PubMed]
15. Waters, L.S.; Storz, G. Regulatory RNAs in bacteria.Cell2009,136, 615–628. [CrossRef]
16. Fröhlich, K.S.; Gottesman, S. Small Regulatory RNAs in the Enterobacterial Response to Envelope Damage and Oxidative Stress.
InRegulating with RNA in Bacteria and Archaea; John Wiley & Sons: Hoboken, NJ, USA, 2018.
17. Raghavan, R.; Groisman, E.A.; Ochman, H. Genome-wide detection of novel regulatory RNAs in E. coli.Genome Res.2011,21, 1487–1497. [CrossRef] [PubMed]
18. Shinhara, A.; Matsui, M.; Hiraoka, K.; Nomura, W.; Hirano, R.; Nakahigashi, K.; Tomita, M.; Mori, H.; Kanai, A. Deep sequencing reveals as-yet-undiscovered small RNAs in Escherichia coli.BMC Genom.2011,12, 428. [CrossRef]
19. Miyakoshi, M.; Matera, G.; Maki, K.; Sone, Y.; Vogel, J. Functional expansion of a TCA cycle operon mRNA by a 30end-derived small RNA.Nucleic Acids Res.2019,47, 2075–2088. [CrossRef]
20. Feng, L.; Rutherford, S.T.; Papenfort, K.; Bagert, J.D.; van Kessel, J.C.; Tirrell, D.A.; Wingreen, N.S.; Bassler, B.L. A qrr noncoding RNA deploys four different regulatory mechanisms to optimize quorum-sensing dynamics.Cell2015,160, 228–240. [CrossRef]
21. Borgmann, J.; Schakermann, S.; Bandow, J.E.; Narberhaus, F. A Small Regulatory RNA Controls Cell Wall Biosynthesis and Antibiotic Resistance.MBio2018,9. [CrossRef]
22. Chao, Y.; Vogel, J. A 3’ UTR-Derived Small RNA Provides the Regulatory Noncoding Arm of the Inner Membrane Stress Response.
Mol. Cell2016,61, 352–363. [CrossRef] [PubMed]
23. Westermann, A.J.; Forstner, K.U.; Amman, F.; Barquist, L.; Chao, Y.; Schulte, L.N.; Muller, L.; Reinhardt, R.; Stadler, P.F.; Vogel, J.
Dual RNA-seq unveils noncoding RNA functions in host-pathogen interactions.Nature2016,529, 496–501. [CrossRef]
24. Zhang, S.; Liu, S.; Wu, N.; Yuan, Y.; Zhang, W.; Zhang, Y. Small Non-coding RNA RyhB Mediates Persistence to Multiple Antibiotics and Stresses in Uropathogenic Escherichia coli by Reducing Cellular Metabolism. Front. Microbiol. 2018,9, 136.
[CrossRef] [PubMed]
25. Massé, E.; Gottesman, S. A Small RNA Regulates the Expression of Genes Involved in Iron Metabolism in Escherichia coli.Proc.
Natl. Acad. Sci. USA2002,99, 4620–4625. [CrossRef]
26. Smirnov, A.; Forstner, K.U.; Holmqvist, E.; Otto, A.; Gunster, R.; Becher, D.; Reinhardt, R.; Vogel, J. Grad-seq guides the discovery of ProQ as a major small RNA-binding protein.Proc. Natl. Acad. Sci. USA2016,113, 11591–11596. [CrossRef] [PubMed]
27. Müller, P.; Gimpel, M.; Wildenhain, T.; Brantl, S. A new role for CsrA: Promotion of complex formation between an sRNA and its mRNA target in Bacillus subtilis.RNA Biol.2019,16, 972–987. [CrossRef]
28. Vogel, J.; Luisi, B.F. Hfq and its constellation of RNA.Nat. Rev. Microbiol.2011,9, 578–589. [CrossRef] [PubMed]
29. Melamed, S.; Peer, A.; Faigenbaum-Romm, R.; Gatt, Y.E.; Reiss, N.; Bar, A.; Altuvia, Y.; Argaman, L.; Margalit, H. Global Mapping of Small RNA-Target Interactions in Bacteria.Mol. Cell2016,63, 884–897. [CrossRef]
30. Iosub, I.A.; van Nues, R.W.; McKellar, S.W.; Nieken, K.J.; Marchioretto, M.; Sy, B.; Tree, J.J.; Viero, G.; Granneman, S. Hfq CLASH uncovers sRNA-target interaction networks linked to nutrient availability adaptation.Elife2020,9, e54655. [CrossRef]
31. Holmqvist, E.; Li, L.; Bischler, T.; Barquist, L.; Vogel, J. Global Maps of ProQ Binding In Vivo Reveal Target Recognition via RNA Structure and Stability Control at mRNA 3’ Ends.Mol. Cell2018,70, 971–982 e976. [CrossRef]
32. Pircher, A.; Gebetsberger, J.; Polacek, N. Ribosome-associated ncRNAs: An emerging class of translation regulators.RNA Biol.
2014,11, 1335–1339. [CrossRef] [PubMed]
33. Gonskikh, Y.; Gerstl, M.; Kos, M.; Borth, N.; Schosserer, M.; Grillari, J.; Polacek, N. Modulation of mammalian translation by a ribosome-associated tRNA half.RNA Biol.2020,17, 1125–1136. [CrossRef] [PubMed]
34. Fricker, R.; Brogli, R.; Luidalepp, H.; Wyss, L.; Fasnacht, M.; Joss, O.; Zywicki, M.; Helm, M.; Schneider, A.; Cristodero, M.; et al.
A tRNA half modulates translation as stress response in Trypanosoma brucei.Nat. Commun.2019,10, 118. [CrossRef] [PubMed]
35. Pircher, A.; Bakowska-Zywicka, K.; Schneider, L.; Zywicki, M.; Polacek, N. An mRNA-derived noncoding RNA targets and regulates the ribosome.Mol. Cell2014,54, 147–155. [CrossRef] [PubMed]
36. Gebetsberger, J.; Wyss, L.; Mleczko, A.M.; Reuther, J.; Polacek, N. A tRNA-derived fragment competes with mRNA for ribosome binding and regulates translation during stress.RNA Biol.2017,14, 1364–1373. [CrossRef] [PubMed]
37. Wyss, L.; Waser, M.; Gebetsberger, J.; Zywicki, M.; Polacek, N. mRNA-specific translation regulation by a ribosome-associated ncRNA in Haloferax volcanii.Sci. Rep.2018,8, 12502. [CrossRef]
38. Hoogstrate, Y.; Jenster, G.; Martens-Uzunova, E.S. FlaiMapper: Computational annotation of small ncRNA-derived fragments using RNA-seq high-throughput data.Bioinformatics2015,31, 665–673. [CrossRef]
39. Nitzan, M.; Rehani, R.; Margalit, H. Integration of Bacterial Small RNAs in Regulatory Networks.Annu. Rev. Biophys.2017,46, 131–148. [CrossRef]
40. Zywicki, M.; Bakowska-Zywicka, K.; Polacek, N. Revealing stable processing products from ribosome-associated small RNAs by deep-sequencing data analysis.Nucleic Acids Res.2012,40, 4013–4024. [CrossRef]
41. Janssen, B.D.; Hayes, C.S. The tmRNA ribosome-rescue system.Adv. Protein Chem. Struct. Biol.2012,86, 151–191. [CrossRef]
42. Luidalepp, H.; Berger, S.; Joss, O.; Tenson, T.; Polacek, N. Ribosome Shut-Down by 16S rRNA Fragmentation in Stationary-Phase Escherichia coli.J. Mol. Biol.2016,428, 2237–2247. [CrossRef] [PubMed]
43. Pinel-Marie, M.L.; Brielle, R.; Riffaud, C.; Germain-Amiot, N.; Polacek, N.; Felden, B. RNA antitoxin SprF1 binds ribosomes to attenuate translation and promote persister cell formation in Staphylococcus aureus.Nat. Microbiol.2021,6, 209–220. [CrossRef]
[PubMed]
44. Wassarman, K.M.; Storz, G. 6S RNA regulates E. coli RNA polymerase activity.Cell2000,101, 613–623. [CrossRef]
45. Steuten, B.; Hoch, P.G.; Damm, K.; Schneider, S.; Köhler, K.; Wagner, R.; Hartmann, R.K. Regulation of transcription by 6S RNAs:
Insights from the Escherichia coli and Bacillus subtilis model systems.RNA Biol.2014,11, 508–521. [CrossRef]
46. Han, K.; Kim, K.S.; Bak, G.; Park, H.; Lee, Y. Recognition and discrimination of target mRNAs by Sib RNAs, a cis-encoded sRNA family.Nucleic Acids Res.2010,38, 5851–5866. [CrossRef] [PubMed]
47. Jahanshahi, S.; Li, Y. An Effective Method for Quantifying RNA Expression of IbsC-SibC, a Type I Toxin-Antitoxin System in Escherichia coli.Chembiochem2020,21, 3120–3130. [CrossRef]
48. Rasmussen, A.A.; Johansen, J.; Nielsen, J.S.; Overgaard, M.; Kallipolitis, B.; Valentin-Hansen, P. A conserved small RNA promotes silencing of the outer membrane protein YbfM.Mol. Microbiol.2009,72, 566–577. [CrossRef] [PubMed]
49. Mandin, P.; Gottesman, S. A genetic approach for finding small RNAs regulators of genes of interest identifies RybC as regulating the DpiA/DpiB two-component system.Mol. Microbiol.2009,72, 551–565. [CrossRef] [PubMed]
50. Weichart, D.; Lange, R.; Henneberg, N.; Hengge-Aronis, R. Identification and characterization of stationary phase inducible genes in Escherichia coli.Mol. Microbiol.1993,10, 407–420. [CrossRef] [PubMed]
51. Polacek, N.; Ivanov, P. The regulatory world of tRNA fragments beyond canonical tRNA biology.RNA Biol.2020,17, 1057–1059.
[CrossRef] [PubMed]
52. Cristodero, M.; Polacek, N. The multifaceted regulatory potential of tRNA-derived fragments.Non-Coding RNA Investig.2017, 170, 61–71.e11. [CrossRef]
53. Goodarzi, H.; Liu, X.; Nguyen, H.C.; Zhang, S.; Fish, L.; Tavazoie, S.F. Endogenous tRNA-Derived Fragments Suppress Breast Cancer Progression via YBX1 Displacement.Cell2015,161, 790–802. [CrossRef] [PubMed]
54. Kim, H.K.; Fuchs, G.; Wang, S.; Wei, W.; Zhang, Y.; Park, H.; Roy-Chaudhuri, B.; Li, P.; Xu, J.; Chu, K.; et al. A transfer-RNA- derived small RNA regulates ribosome biogenesis.Nature2017,552, 57–62. [CrossRef]
55. Oberbauer, V.; Schaefer, M.R. tRNA-Derived Small RNAs: Biogenesis, Modification, Function and Potential Impact on Human Disease Development.Genes2018,9, 607. [CrossRef] [PubMed]
56. Lalaouna, D.; Carrier, M.C.; Semsey, S.; Brouard, J.S.; Wang, J.; Wade, J.T.; Masse, E. A 3’ external transcribed spacer in a tRNA transcript acts as a sponge for small RNAs to prevent transcriptional noise.Mol. Cell2015,58, 393–405. [CrossRef] [PubMed]
57. Melamed, S.; Adams, P.P.; Zhang, A.; Zhang, H.; Storz, G. RNA-RNA Interactomes of ProQ and Hfq Reveal Overlapping and Competing Roles.Mol. Cell2020,77, 411–425.e417. [CrossRef]
58. Blattner, F.R.; Plunkett, G.; Bloch, C.A.; Perna, N.T.; Burland, V.; Riley, M.; Collado-Vides, J.; Glasner, J.D.; Rode, C.K.; Mayhew, G.F.; et al. The Complete Genome Sequence ofEscherichia coliK-12.Science1997,277, 1453–1462. [CrossRef]
59. Vimberg, V.; Tats, A.; Remm, M.; Tenson, T. Translation initiation region sequence preferences in Escherichia coli.BMC Mol. Biol.
2007,8, 100. [CrossRef] [PubMed]
60. Luidalepp, H.; Joers, A.; Kaldalu, N.; Tenson, T. Age of inoculum strongly influences persister frequency and can mask effects of mutations implicated in altered persistence.J. Bacteriol.2011,193, 3598–3605. [CrossRef] [PubMed]
61. Lee, T.S.; Krupa, R.A.; Zhang, F.; Hajimorad, M.; Holtz, W.J.; Prasad, N.; Lee, S.K.; Keasling, J.D. BglBrick vectors and datasheets:
A synthetic biology platform for gene expression.J. Biol. Eng.2011,5, 12. [CrossRef] [PubMed]
62. Miyazaki, K. Chapter seventeen-MEGAWHOP Cloning: A Method of Creating Random Mutagenesis Libraries via Megaprimer PCR of Whole Plasmids. InMethods in Enzymology; Voigt, C., Ed.; Academic Press: Cambridge, MA, USA, 2011; Volume 498, pp.
399–406.
63. Qi, D.; Scholthof, K.B. A one-step PCR-based method for rapid and efficient site-directed fragment deletion, insertion, and substitution mutagenesis.J. Virol. Methods2008,149, 85–90. [CrossRef] [PubMed]
64. Gebetsberger, J.; Zywicki, M.; Kunzi, A.; Polacek, N. tRNA-derived fragments target the ribosome and function as regulatory non-coding RNA in Haloferax volcanii.Archaea2012,2012, 260909. [CrossRef] [PubMed]
65. Rao, X.; Huang, X.; Zhou, Z.; Lin, X. An improvement of the 2ˆ (-delta delta CT) method for quantitative real-time polymerase chain reaction data analysis.Biostat. Bioinforma. Biomath.2013,3, 71–85.
66. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. EdgeR: A Bioconductor package for differential expression analysis of digital gene expression data.Bioinformatics2010,26, 139–140. [CrossRef] [PubMed]
67. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads.Embnet J.2011,17, 3. [CrossRef]
68. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2.Nat. Methods2012,9, 357–359. [CrossRef] [PubMed]
69. Edgar, R.C.; Flyvbjerg, H. Error filtering, pair assembly and error correction for next-generation sequencing reads.Bioinformatics 2015,31, 3476–3482. [CrossRef] [PubMed]
70. Mohanty, B.K.; Kushner, S.R. The majority of Escherichia coli mRNAs undergo post-transcriptional modification in exponentially growing cells.Nucleic Acids Res.2006,34, 5695–5704. [CrossRef] [PubMed]
71. Fu, H.; Feng, J.; Liu, Q.; Sun, F.; Tie, Y.; Zhu, J.; Xing, R.; Sun, Z.; Zheng, X. Stress induces tRNA cleavage by angiogenin in mammalian cells.FEBS Lett.2009,583, 437–442. [CrossRef]
72. Liao, Y.; Smyth, G.K.; Shi, W. featureCounts: An efficient general purpose program for assigning sequence reads to genomic features.Bioinformatics2014,30, 923–930. [CrossRef]
73. Li, W.; Godzik, A. Cd-hit: A fast program for clustering and comparing large sets of protein or nucleotide sequences.Bioinformatics 2006,22, 1658–1659. [CrossRef] [PubMed]
