Dynamic insights on transcription initiation and RNA processing during bacterial adaptation.

Transcription initiation and RNA processing govern gene expression and enable bacterial adaptation by reshaping the RNA landscape. The aim of this study was to simultaneously observe these two fundamental processes in a transcriptome responding to an environmental signal. A controlled σE system in E. coli was coupled to our previously described tagRNA-seq method to yield process kinetics information. Changes in transcription initiation frequencies (TIF) and RNA processing frequencies (PF) were followed using 5' RNA tags. Changes in TIF showed a binary increased/decreased pattern that alternated between transcriptionally activated and repressed promoters, providing the bacterial population with transcriptional oscillation. PF variation fell into three categories of cleavage activity: (i) constant and independent of RNA levels, (ii) increased once RNA has accumulated, and (iii) positively correlated to changes in TIF. This work provides a comprehensive and dynamic view of major events leading to transcriptomic reshaping during bacterial adaptation. It unveils an interplay between transcription initiation and the activity of specific RNA cleavage sites. This study utilized a well-known genetic system to analyze fundamental processes and can serve as a blueprint for comprehensive studies that exploit the RNA metabolism to decipher and understand bacterial gene expression control.

INTRODUCTION

In response to an environmental cue, gene expression reprogramming triggers changes in RNA synthesis, RNA processing and/or degradation, resulting in a physiological response. Transcription initiation and RNA processing frequencies for genes not involved in these switches remain unchanged during the adaptation process and their primary and processed RNAs remain at constant levels. In contrast, for genes mediating the adaptation process, the transcription initiation frequency and/or the RNA processing frequency vary resulting in novel or transient changes in RNA levels (e.g., Reznikoff et al. 1985; Phadtare and Severinov 2010; Rochat et al. 2013; Bouloc and Repoila 2016). These events occur in minutes, or less, and in a fraction of the bacterial population doubling time (Anderson and Dunman 2009; Esquerré et al. 2014).

Transcriptomic approaches are suited to provide information on adaptation processes in a genome-wide manner as they enable comparisons of RNA landscape snapshots between different environmental conditions or genetic backgrounds (Croucher and Thomson 2010; Filiatrault 2011; Güell et al. 2011; Mader et al. 2011). The majority of these studies compare steady-state RNA levels in bacterial offspring against the initial bacterial population. As such, they lack information on the dynamics of RNA remodeling, the chronology of events, and the mechanisms responsible for changes in gene expression. Pioneering studies addressing RNA dynamics during adaptation processes have been reported. For instance, in Caulobacter crescentus, activations and repressions of gene expression during the cell cycle development was visualized using DNA microarrays (Laub et al. 2000). In Escherichia coli K12, the expression kinetics of regulons under control of the extracytoplasmic sigma factor σE and the noncoding RNA RhyB were observed (Masse et al. 2005; Rhodius et al. 2006; Bury-Mone et al. 2009; Gogol et al. 2011). However, these studies remain unable to distinguish transcriptional or post-transcriptional effects on RNA levels. Similarly, genomic-scale measurements of RNA stability uncouple transcription from RNA processing and degradation by using rifampicin, an antibiotic blocking transcription initiation (Mosteller and Yanofsky 1970; Redon et al. 2005; Kristoffersen et al. 2012; Esquerré et al. 2014; Dar and Sorek 2018). This technique cannot map the RNA processing sites that are key actors in RNA decay (Mohanty and Kushner 2016), and may mask interplay between transcription, RNA processing and degradation as demonstrated in eukaryotic organisms (Dahan and Choder 2013; Singh et al. 2015; Peck et al. 2019). RNA-seq methods have been gradually adapted in an effort to probe regulatory mechanisms of bacterial transcriptomes. Differential RNA-seq (dRNA-seq) compares total transcriptomes to those enriched for 5′ triphosphate RNA ends to aid in transcription start sites (TSSs) prediction (Sharma and Vogel 2014). Genome-wide single-nucleotide-resolution mapping of 5′ RNA termini can be performed by methods such as tagging RNA-seq (tagRNA-seq) or EMOTE (Fouquier d'Herouel et al. 2011; Linder et al. 2014; Innocenti et al. 2015). We demonstrated, in E. coli K12 and Enterococcus faecalis, that tagRNA-seq can map and distinguish TSSs and RNA processing sites (PSSs), and that tagging efficacy is proportional to the abundance of 5′ termini (Fouquier d'Herouel et al. 2011; Innocenti et al. 2015). These new techniques, while powerful, come with several caveats. dRNA-seq relies on the 5′-phosphate-dependent exonuclease (TEX) to degrade RNAs from 5′ monophosphate groups, yet TEX is sensitive to RNA folding and can generate artefactual internal TSSs (Szittya et al. 2010; Conway et al. 2014; Prados et al. 2016). Along the same lines, tagRNA-seq generates a significant number of RNA ends, coined “undetermined” (UND) that cannot be classified as TSS or PSS without additional and individual analysis (Innocenti et al. 2015). Nevertheless, tagRNA-seq provides the most promising means to investigate the dynamics of key events reshaping the transcriptome during the bacterial adaptation process.

The maintenance of envelop homeostasis in E. coli is ensured by regulatory networks including the σE regulon (Ades 2008; Silhavy et al. 2010; Grabowicz and Silhavy 2016). The extracytoplasmic σ factor (σE) encoded by rpoE is normally sequestered by RseA at the inner surface of the cytoplasmic membrane. Perturbations of surface proteins or envelop integrity cause σE to be released into the cytoplasm where it binds to RNA polymerase (RNAP). In turn, σE-RNAP binds to specific promoters directing transcription of ∼100 encoding sequences (CDSs) including surface proteins, enzymes, envelop compounds, transcription and translation components, and rpoE itself. Regulatory factors, including three small RNAs (sRNAs), MicA, RybB and MicL, are also expressed. These sRNAs provide negative feedback to the σE regulon, modulating the translation and/or the stability of surface proteins mRNAs (Rhodius et al. 2006, 2012; Mutalik et al. 2009; Gogol et al. 2011; Guo et al. 2014; Shimada et al. 2017).

The aim of this study was to evidence at the genomic scale changes in transcription initiation frequency (TIF) and RNA processing frequency (PF) during an adaptation process, the activities of two key molecular events generating 5′ RNA ends. The σE regulon exhibits the crucial aspects of RNA metabolism governing gene expression control: RNA synthesis, processing, and degradation. It consists of a manageable number of genes, making it an attractive model to simultaneously observe RNA transcription initiation and RNA processing for the first time at the genomic scale in bacteria. We used tagRNA-seq to study the kinetics of the RNA pool within a single bacterial generation in response to σE induction. The results presented, including the synchrony between changes in the cleavage efficacy measured at specific PSSs, and the activity of transcription initiation at the cognate TSSs, provide evidence to support the hypothesis that interplay exists between transcription initiation and RNA cleavage.

RESULTS AND DISCUSSION

Dynamics of the RNA levels in response to σE induction

The σE regulon was used as a model system in this study to observe an evolving bacterial RNA landscape. The σE encoding sequence, rpoE, was expressed under the control of an inducible promoter responding to anhydrotetracycline (aTc) in E. coli K12 (Bury-Mone et al. 2009). After the addition of aTc, the dynamics of the RNA pool was monitored by collecting samples every 5 min over a period of 20 min. Samples were then analyzed by tagRNA-seq (Fig. 1A; Innocenti et al. 2015). RNA-seq raw data were processed as presented in the “Materials and Methods” section. Supplemental Tables S1, S2 provide a comprehensive list of normalized data, statistical treatments, and temporal changes in RNA levels. As detailed and discussed in Supplemental Section S1, changes observed in RNA levels over the experimental time course demonstrated that the σE regulon was induced and functional, a prerequisite for this study. The bulk of gene reprogramming was observed within the 10 first min and was nearly complete at 20 min (Fig. 1B). The majority of changes in RNA levels showed a monotonic variation (increased or decreased) and a minority was “transient.” mRNA targets of sRNAs synthesized by σE-RNAP (MicA, RybB, and MicL) showed differential sensitivity to sRNA-mediated degradation. These observations confirmed that our σE experimental system produced the expected gene expression reprogramming patterns that had been previously reported (Supplemental Section S1; Supplemental Tables S2, S3). The results further demonstrated that the inducible σE system was functional and that resulting RNA landscape kinetics could be used to analyze changes in TIF and PF.

FIGURE 1.

Changes in transcription initiation frequency at σE-TSSs. (A) Principle of the tagRNA-seq method used in a time course experiment. Aliquots of growing bacteria were harvested every 5 min over 20 min following the addition of aTc. Total RNA was extracted and sequenced by tagRNA-seq. RNA samples at t0 were duplicated: one treated with tobacco alkaline phosphatase (+TAP); the other was not treated (TAP no). (B) Number of CDSs whose RNA level varies at least fourfold compared to t0. The selection process is nonrepetitive: A CDS selected at one time point was excluded from the selection process at subsequent time points. Total CDSs selected are shown by the blue histogram. Among those, the red histogram indicates the number of selected CDSs known or predicted to be transcribed by σE-RNAP. (C) Heatmap correlating patterns of tagging acceleration values (A) measured at σE-TSSs assigned over the course of the experiment (Table 2). (*) A values are below the statistical confidence (pA) or differences in A values are below the biological threshold of variation imposed (>|0.7| tag/min²). D–G refer to graphs on the left of the figure. Numbers attached to TSSs refer to chromosomal coordinates. (D) Biphasic pattern of changes in TIF (A values expressed in tag/min²) observed for about half of σE-TSSs mapped. (E–G) Variation from the biphasic patterns for certain σE-TSSs mapped. (X-axis, min) Experimental time points of the kinetics. (Y-axis) A values measured. Supplemental Figure S3 shows the patterns of changes in TIF for reported σE-TSSs mapped as UNDs and PSSs in this study.

TABLE 2.

Tagging accelerations measured at 5′ RNA ends mapped and previously reported as σE-TSSs

Frequencies of processes generating 5′ RNA ends

TagRNA-seq discriminates 5′ RNA ends generated by transcription initiation and RNA cleavage events by use of two short RNA sequences termed “TSS-” and “PSS-tag.” Criteria utilized to assign 5′ RNA ends are found in Supplemental Section S2 and the modulation of tags in the RNA pool is summarized in Supplemental Table S4. We identified 1147 5′ RNA ends as TSSs, and 594 as PSSs, while 1706 remained UNDs (Table 1). We previously demonstrated that tag-counts provide good estimates for the relative abundance of 5′ termini (Fouquier d'Herouel et al. 2011; Innocenti et al. 2015). Although current RNA-seq methodologies exhibit strong technical variability (McIntyre et al. 2011; Evans et al. 2018), the consistency analysis of tag-counts between kinetic series indicated that mean tag-counts reflect dynamics despite the variability between individual series (see Materials and Methods and Supplemental Table S4). Based on this result, we report mean tag-count values that relate to frequencies of events generating 5′ RNA ends. We defined the “tagging rate” (R) as the average number of tags appearing or disappearing over time, a magnitude coupled to the frequency of the process generating any given RNA terminus. An initial steady state of expression is assumed prior to the induction of σE where R0 = 0 tag/min. The tagging rate over an interval of time, Rinterval is Ri = [(ni–ni − 1)/(ti–ti − 1)]; where ni and ni − 1 are mean tag-counts at corresponding time points ti and ti − 1. From Ri, we focused on changes in tagging frequency as a metric of gene reprogramming. Changes in tagging frequency are evaluated by comparing Ri at different time points to yield a “tagging acceleration” (A). For consecutive experimental time intervals, Ainterval is Ai = [(Ri–Ri − 1)/(ti–ti − 1)]. As with R0, the initial tagging acceleration is assumed to be A0 = 0 tag/min². Changes in tagging frequency (Ai) can be positive or negative. When Ai > 0, it indicates an increased frequency, when Ai < 0, it indicates a decreased frequency, when Ai = 0, it indicates a frequency equal to the one in the previous time interval (Ri = Ri − 1). Thus, if a particular 5′ end is a TSS, the Ai values correspond to changes in transcription initiation frequency (changes in TIF). Ai values for a PSS correspond to changes in the RNA processing frequency (changes in PF). Supplemental Table S5 summarizes mean values for Ri and Ai for each 5′ RNA end mapped. Since only three biological replicates were used and that RNA treatments feature large technical variability, measured tag-counts and calculated Ri and Ai values were typified by high dispersion, as expected and shown by standard deviations (Supplemental Tables S4, S5). Consequently, the robustness of mean values was evaluated by statistical methods applied to small numbers of samples (Supplemental Tables S4, S5). We calculated confidence values pAi, expressed as meta-analysis P-value using Fisher's method. pAi values <1.45 × 10−5 are considered significant, with Ai values at least 95% reliable. Out of 3447 total 5′ termini mapped, this level of confidence applies for ∼40% of RNA ends.

TABLE 1.

5′ RNA ends showing changes in frequency

Assessment of a threshold for biological significance

Deducing biological relevance from changes in frequency (tagging acceleration) calculations is a challenge without a similar guiding precedent. Here, we leveraged the well-studied σE system and the fact that an increased TIF was anticipated for σE-dependent TSSs (σE-TSSs) in order to establish an empirical “biological significance” threshold for changes in Ai. Out of 88 reported (confirmed experimentally or predicted) σE-TSSs (Gama-Castro et al. 2016; Keseler et al. 2017), 40 were mapped, including 22 assigned as TSSs (Table 2). A5 values (5 min following σE induction) for this set of promoters varied greatly from 0.1 tag/min2 (P) to 21.1 tag/min2 (P) and prevented direct determination of a threshold. This result was unsurprising given the large dynamic range of bacterial promoters, including σE-dependent ones (Mutalik et al. 2009; Rhodius and Mutalik 2010). However, since it may be that certain σE-TSS did not respond to the σE induction, we imposed an arbitrary threshold based on the lowest A5 value: |Ai–Ai − 1| > 0.7 tag/min2, or 7 times the value for P. While arbitrary, such a threshold was effectively stringent as it excluded σE-dependent P from promoters called as transcriptionally activated during the first 5 min (t5–0). Empirical significance of this threshold is strengthened as over 80% of RNA ends with a significant biological change at t5 were retained as statistically significant (Supplemental Table S5). This cutoff was subsequently applied to 5′RNA ends over the course of the experiment, including RNA cleavage sites (PSSs). Out of all mapped 5′ RNA ends, approximately one third (38%) showed a significant change over 20 min of measurements. The bulk of significant TSSs changes were found at the earliest interval, t5–0. In comparison, the bulk of responding PSSs peaked later at t10 (Table 1). The observation that most gene reprogramming occurred within the first 10 min following induction fits with the conclusion inferred from RNA levels (Supplemental Section S1; Fig. 1B). Significant changes in UND 5′ends mirrored the pattern of PSSs, most likely due to our assignment method that favors the presence of true PSSs in this group (Supplemental Section S2).

By applying the threshold value of >|0.7| tag/min2, we retrieved tendencies observed for changes in RNA levels, indicating that such an empirical threshold is suited for a global analysis. To pinpoint general features and highlight specific outputs from changes in TIF and PF, we then focused on patterns described by changes in A values. These are presented as selected examples. Unless otherwise mentioned, all satisfy the biological significance and the statistical confidence.

Pattern of changes in TIF at σE-dependent promoters

Significant changes in Ai values were recorded for 21/22 assigned σE-TSSs, the sole exception being P (Table 2). Nine TSSs (P, P, P, P, P, P, P, P, and P) exhibited a biphasic response, with increased TIF observed at intervals t5–0 and t15–10, and decreased values at t10–5 and t20–15. The second increase at t15-10 was of similar or higher magnitude. We also observed that decreases in TIF were similar or higher compared with increases, indicating that transcription initiation tended to return to the t0 state (Table 2; Fig. 1C,D). The biphasic pattern reflects a “transcriptional oscillation” within the bacterial population, initially synchronized by addition of the inducer at t0. This pattern could reflect transcriptional bursting occurring at the single cell level for highly expressed genes where initiation occurs in surges and promoters exhibit intermittent periods of inactivity possibly due to topological constraints (Golding et al. 2005; So et al. 2011; Chong et al. 2014).

Four additional promoters (P, P, P, and P) showed a delayed response with increased TIF at t10 instead of t5, followed by decreased- and increased TIFs (Fig. 1C,E). Three others (P, P, and P) displayed an increased/decreased TIF pattern during t15–0, but then continued in increase in t20-15 (Fig. 1C,F). The remaining five, from the original 21 (P, P, P P, and P), showed an increased TIF that peaked at t10 or t15, followed by a drastic decrease (Fig. 1C,G). Variations from the biphasic response most likely reflect additional controls modulating the use of σE-TSSs by σE-RNAP. Changes in TIF for P and P, for instance, were rather weak compared to other σE-TSSs, indicating low transcriptional activation by σE induction (Supplemental Table S5). In contrast, RNA levels for dsbC and ydhIJK were highly increased at t5 and remained 20-fold above baseline from t10 through t20 (Supplemental Table S2). These observations strongly suggest that in addition to transcription, stabilization of these transcripts occurs and can be responsible for the variations in the biphasic pattern.

Among other known σE-TSSs, 13 were assigned as UND in our study (Table 2). Six showed no significant variations in TIF, suggesting little to no activation. The remaining seven generally followed the biphasic pattern, indicating that these promoters were activated by σE induction (P, P, P, P, P, P, and P) (Table 2; Supplemental Fig. S3).

Five previously mapped TSSs (Keseler et al. 2017) were assigned as “PSSs” (Table 2), and only three of them (P P, and P) showed significant variations in A, with similar patterns to those observed for mapped σE-TSSs (Fig. 1E,G; Supplemental Fig. S3). The assignment of these 5′ RNA ends as PSSs, may be due to the transcriptional organization revealed by tagging density at these loci (Supplemental Table S5). lpxD is embedded in a complex operon where at least three promoters are located upstream of the mapped promoter (coordinate 200960; Table 2). Longer transcripts containing lpxD may favor the PSS assignment of the corresponding nucleotide due to RNA processing or degradation generating 5′ monophosphate ends. In line with this possibility, a TSS was mapped two nucleotides downstream, which may be the TSS of lpxD (coordinate 200962; Supplemental Table S5). Similarly, reported TSSs for degP and fkpA (Table 2) were embedded within strongly transcribed regions that featured an abundance of tags ligated to upstream and downstream nucleotides: 28 and 12 consecutive nucleotides were tagged for degP and fkpA, respectively. Based on our data, we predict TSSs at 180840/42 for degP, and 3477528 for fkpA.

In general, the tag-counts dynamics revealed specific patterns of changes in TIF at each σE-TSS under σE induction. Although each σE-TSS did not respond with the similar intensity and can exhibit variations, transcription activation features a general biphasic pattern originating with an increased TIF.

Biphasic patterns at σE-independent promoters

Half of mapped TSSs responded to σE induction, with ∼75% of these responding within the first 10 min (Table 1). We observed two major TIF variations profiles for promoters responding to σE induction but reportedly not transcribed by σE-RNAP (Keseler et al. 2017). A biphasic pattern with increased TIFs at intervals t5–0 and t15–10 and decreases at t10–5 and t20–15 was observed (dps, grxA, rpmHp2, rpsMp2, rpsT, secG, tufB, yabI, and ybhQ2) (Fig. 2A; Supplemental Fig. S4A). As with σE-TSSs, this indicates transcription activation, a conclusion corroborated by increased RNA levels at t5 for most of the corresponding CDSs (Supplemental Table S2). The second pattern was typified by decreases at intervals t5-0 and t15-10 (cspD, gatY, glpT, hupA, lpp, polA, sodA, ssrSp1, and treB) (Fig. 2B; Supplemental Fig. S4B). This latter profile, with decreased TIF at the interval t5–0 is likely to indicate transcriptional repression. The half-life of most E. coli mRNAs is <5 min and so RNA degradation has a prominent impact on RNA levels (Chen et al. 2015; Dar and Sorek 2018). Reduced synthesis during t5–0 would, therefore, be sufficient to decrease the amount of RNA. Decreased levels of transcript corroborate this conclusion, except for relatively stable transcripts (lpp and ssrS) (Supplemental Table S1; Supplemental Section S1).

FIGURE 2.

Changes in transcription initiation frequency at σE-independent TSSs. (A,B) Heatmap correlating patterns of changes in TIF (A values) at selected TSSs mapped and also reported in Gama-Castro et al. (2016) and Keseler et al. (2017). (A) Transcription activation; (B) Transcription repression. (C) Graphical representation of changes in TIF indicating transcription activation and repression for TSSs dps-848948 and gatY-2177231, respectively. (D) RNA levels for dps and gatY. Y-axes (RNA): log2 of the ratio between RNA amounts at an experimental time point (ti) and t0. Legends are otherwise identical to Figure 1.

From the data, we hypothesize that the opposing patterns of changes in TIF could result from competition between sigma factors binding to core-RNAP, with contribution from promoter strengths, transcription regulators, and transcriptional burst that temporarily prevent promoters from being reused by the RNAP (Jishage et al. 1996; Golding et al. 2005; Chong et al. 2014; Mauri and Klumpp 2014). In the simplest model, promoters with an increased activity at t5 would outcompete those displaying a decreased TIF. During the first minutes after σE induction, the activated promoters would be occupied and/or refractory to reinitiation and the completion of ongoing transcription would release core-RNAP. The situation would then reverse, with the binding of free sigma factors to core-RNAP and the utilization of the “unoccupied” promoters during the interval t5–0. This would generate opposite changes in TIF between the two groups of promoters and transcription oscillation to the bacterial population, a hypothesis consistent with current modeling of transcriptionally active promoters indicating ON/OFF states at the single cell level (Jones and Elf 2018).

About 15% of total mapped TSSs showed significant changes in TIF after t10 (Table 1). These promoters may themselves be a response to cellular variations provoked by the increase of σE. However, many of these TSSs are flanking nucleotides to TSSs detected in previous time intervals. For instance, the TSS for yccA, reported at 1031443 (Keseler et al. 2017) and encoding a modulator of FtsH protease, was also tagged at 1031444/5/6/7/8: Two positions showed no significant changes in TIF, two others decreased at t5–0 and one decreased at t15–10 (Supplemental Table S5). Other examples were also observed (gatY-2177230/1/2/3, yobF-1907615/6/7/8, yceDp1-1146648/49/50) (Supplemental Table S5), and may suggest that TIF changes can modify stringency of the holo-RNAP to utilize specific nucleotides within the DNA promoter.

Further analysis combining changes in TIF and RNA levels revealed additional transcriptional and post-transcriptional effects. For instance, the comparison of TSSs and RNA levels of dps and gatY (Fig. 2C,D) showed that RNA levels track changes in TIF for both genes in the interval t5–0, but then diverge. dps RNA synthesis is accompanied by an increased degradation as marked by decreasing RNA levels. The decreased TIF for gatY during the first 5 min would be sufficient to decrease RNA levels without a change in RNA stability. After t5, changes in TIF marked a global increase that may be responsible for the increasing gatY levels observed onward from t10 (Fig. 2C,D). Both dps and gatY reach similar RNA levels at t20 (two- to threefold lower) but their respective paths likely involved differing controls. This highlights the interest of combining kinetics and 5′ tagging approaches with RNA-seq to shed light on the mechanistic aspects of gene expression control.

Selectivity of RNA processing sites

We examined the activity of RNA cleavage at PSSs and identified three distinct classes corresponding to PF changes. The first class, representing ∼60% of total PSSs (Table 1), included PF changes below the threshold of biological significance (|A–A| ≤ 0.7 tag/min²), indicating cleavage independent of any variation in RNA amounts. For activated genes, most of the newly synthesized RNA will be not cleaved at these sites, which suggests stabilization. The second class, ∼30% of total PSSs, included PSSs with late changes in PF relative to variations of measurable RNA amounts (Table 1, three right columns). At these PSSs, cleavage activity responded to RNA accumulation and PF peaked when RNA reached maximal amounts, primarily at t10, t15, and t20. The remaining third PSSs, ∼10% of total PSSs, included early changes in PF during the first 5 min (column “t5” in Table 1). These changes paralleled patterns of increasing RNA levels or variations in TIF when transcription activation was counterbalanced by RNA degradation. For transcriptionally activated genes, RNA cleavage most likely occurs concomitantly to RNA synthesis. Each of the three classes, with examples is discussed in the following sections.

RNA dynamics at the rpsU-dnaG-rpoD operon

In bacteria, polycistronic transcripts encode proteins which may be required in different amounts. A common solution is separation via cleavage coupled to differential RNA decay for each CDS component (Rochat et al. 2013). The E. coli σE-independent rpsU-dnaG-rpoD operon is one such example that encodes the ribosomal protein S21, the primase DnaG, and the vegetative sigma factor σD (Fig. 3A). The dnaG mRNA is rendered unstable in comparison to rpoD by RNase E cleavage, maintaining the amount of DnaG at about 100 molecules, versus a few thousands for σD (Burton et al. 1983; Yajnik and Godson 1993). Interestingly, amounts of σD increase under σE induction (Rhodius et al. 2006). We hypothesized that increased dnaG-rpoD RNA levels should be accompanied by a commensurate increased PF at dnaG PSSs such as to preserve the DnaG and σD protein ratios.

FIGURE 3.

RNA dynamics at the rpsU-dnaG-rpoD locus. (A) Locus organization. Black thick arrows indicate CDSs and transcription orientation. Thin horizontal arrows show TSSs mapped and coordinates on the E. coli K12 chromosome MG1655. The blue color for TSSs indicates significant changes in TIF during the σE-mediated adaptation; gray color marks the absence of significant changes in TIF. Sigma factors known to be involved in the activity of promoters are in brackets. Vertical arrows indicate PSSs mapped, in bold and in color those with significant changes in PF, in gray, those with no significant changes. (B) RNA levels over the course of the experiment. (C) Pattern of changes in TIF at TSSs mapped; colors correspond to those used in panel A. (D) Pattern of changes in PF for PSSs mapped. Colors correspond to those used in panel A. Legends are otherwise identical to Figures 1 and 2.

Three TSSs were mapped at the rpsU-dnaG-rpoD locus: two σD-TSSs upstream of rpsU at coordinates 3210646/7 and 3210716/7/8, and one σH-TSS within dnaG, assigned as UND, at position 3212688 (Fig. 3A). Only TSS-3210716 displayed a significantly increased TIF during the interval t10–5, which was insufficient to explain the increased abundance of RNA levels observed at t5 (Fig. 3B,C).

Eleven PSSs were mapped: Nine nested in transcripts harboring dnaG and rpoD sequences and two within the 5′ untranslated region (5′UTR) of rpsU (Fig. 3A). Four (-3212541, -3212942, -3213927, and -3214086) had no significant changes in PF. PSSs-3213710 and -3213842 within rpoD and PSS-3212941 in between dnaG and rpoD CDSs showed modest changes in PFs at t15 and t20 suggesting that RNA cleavage varies due to RNA accumulation (Fig. 3D). Changes in PF for PSS-3210759, within the 5′UTR of rpsU, also had a delayed activity relative to RNA synthesis, peaking at t15 once the global rpsU RNA amounts were declining (Fig. 3B,D). The final two (PSSs-3212852/53) within the dnaG stop codon, showed progressive and significant increased PF that paralleled the increasing amounts of dnaG-rpoD RNA. It should be noted that the RNase E site acting at the end of dnaG was previously reported at position 3212852 (Burton et al. 1983). This indicates that cleavages at these PSSs affects the newly synthesized dnaG-rpoD transcripts (Fig. 3B,D). Upon closer investigation, PF for PSS-3212853 also manifested a response that peaked once dnaG-rpoD RNAs had reached their maximal amounts at t15, (Fig. 3B,D), strongly suggesting a cleavage activity linked to RNA accumulation. However, PF for PSS-3212852 mirrored the increasing amounts of dnaG-rpoD RNAs and peaked with RNA levels at t10, which parallels the increased TIF observed for TSS-3210716 at t10–5 (Fig. 3C,D). The concomitant increase of RNA amounts, PF and TIF at t10–5 strongly suggests that PSS-3212852 responded to RNA synthesis and, possibly, to transcription activation.

Our results confirm the importance of the RNase E site (3212852) in the control of rpsU-dnaG-rpoD operon (Burton et al. 1983; Yajnik and Godson 1993). Beyond this, we are able to detect that PSS-3212852 responds to the increasing RNA synthesis, in contrast to other PSSs which depend on accumulated RNA or act independently of RNA levels. To our knowledge, this is the first observation of differential activity of PSSs within a bacterial operon using genome scale RNA-seq. Most importantly, detection was performed without altering transcription (by rifampicin) or by shifting environmental conditions to inactivate an essential RNase (RNase E).

Changes in PFs for PSSs in dnaG-rpoD transcripts revealed that, within a transcript, certain PSSs respond specifically and differently to changes in RNA levels and RNA synthesis. The evidence for such differentiation is supported by the independence of RNA level and responding PFs across the transcriptome. For instance, mreC, phoBR, ptrA, sbmA, or yabI RNAs reached their highest amounts at t5 and harbor PSSs whose changes in PF were not significant and/or occurred at later time intervals (Supplemental Tables S1, S5). A further example of different classes of PSSs in the σE-dependent bicistronic operon bepA-yfgD is presented in Supplemental Section S4.

RNA dynamics at the ahpCF operon

We wished to test reports that transcription activation of certain genes is overbalanced by increased RNA degradation resulting in net lower levels of RNA (Redon et al. 2005; Esquerré et al. 2014; Nouaille et al. 2017). If such a control occurred during σE induction, we expected: (i) increasing RNA levels at early time points, followed by a degradation-mediated decrease, and (ii) changes in TIF featuring activation coupled with increased PFs of newly synthesized RNA molecules. This implies synchronous changes in TIF and PF for the 5′ ends of the two distinct RNA molecules, which originate from a common primary transcript.

Several transcripts satisfy these expectations, and the ahpCF operon is presented here as we mapped TSSs at locations previously reported (Fig. 4A; Keseler et al. 2017). ahpCF encodes a hydroperoxidase involved in the detoxification of hydrogen peroxide. Under σE induction, the abundance of each CDS (ahpC and ahpF) increased greater than fourfold in the interval t5–0 and then fell to approximately twofold below that of t0 (Fig. 4B). In response to σE induction, a biphasic pattern was observed for only one of the σD-TSSs (638921) indicating that it was responsible for the majority of RNA increase measured at t5 (Fig. 4C). Six PSSs were assigned within the transcription unit ahpCF. Four PSSs showed no significant (638762, 638907) or mild and late (638763, 639731) changes in PF indicating that RNA cleavage occurred independently of RNA amounts at these sites. In contrast, PSSs upstream of ahpC (638922) and between ahpC and ahpF (639732) displayed changes in PFs that paralleled changes in TIFs measured for TSS-638921 (Fig. 4C,D). However, we are cautious in our conclusions here as TSS-638921 and PSS-638922 are separated by one nucleotide, and so we cannot fully rule out that the later might be a TSS explaining why we observed a biphasic pattern. In contrast, for PSS-639732, the similarity between changes in PF and in TIF patterns strongly suggests that RNA cleavages responded to increased TIF, and most likely occurred on nascent RNA molecules. Given the requirement for RNA degradation, we remarked on an absence of PSSs mapped inside CDSs ahpC and ahpF. This indicates that our method “captured” only certain RNA cleavage sites and it is probable that short-lived or smaller RNA fragments eluded the tagging method or the RNA-seq protocol.

FIGURE 4.

RNA dynamics at the ahpCF locus. Legends are identical to Figure 3. Note that the proximity of PSS-638922 and TSS-638921 does not allow us to rule out that the PSS is not a TSS. This PSS is presented in bold and gray in panel A, and with dashed lines in panel D.

Both rpsU-dnaG-rpoD and ahpCF operons attest to the specific response of certain PSSs to RNA synthesis during the adaptation process. The difference between patterns presented in Figure 3D and Figure 4D can be explained by the accumulation of cleaved RNA molecules between experimental time points for dnaG-rpoD, and the elimination of cleaved RNA molecules between time points for ahpCF. Most importantly, changes in PF at specific PSSs that parallel changes in TIF indicates that cleavage activity is somehow coupled to transcription. A few previous studies touched upon the relationship between RNA synthesis and stability in E. coli (Chow and Dennis 1994; Chen et al. 2015; Nouaille et al. 2017). Here we were able to confirm this observation while providing additional insight and examples.

A major role for RNA processing and degradation has been long recognized in single gene expression studies (Burton et al. 1983; Gorski et al. 1985; Newbury et al. 1987; Haugel-Nielsen et al. 1996; Ludwig et al. 2001; Repoila and Gottesman 2001; Winkler et al. 2004; Urban and Vogel 2008; Prevost et al. 2011). The contribution of RNA degradation in adjusting gene expression during bacterial adaptation has been highlighted and, more recently, differential mRNA decay was visualized in a genome wide manner as a key actor reshaping the expression of operonic organizations (Dar and Sorek 2018). In E. coli and Lactococcus lactis, Cocaign-Bousquet and colleagues established the major impact of RNA stability to counteract RNA synthesis and adjust growth to nutrient availability (Redon et al. 2005; Esquerré et al. 2014; Nouaille et al. 2017; Dressaire et al. 2018). In S. aureus responding to diverse cues, Dunman and colleagues showed that more than 80% of stress-modulated transcripts have modified stability during adaptation and that factors of RNA processing and degradation can be targets of antimicrobial molecules (Anderson et al. 2006; Olson et al. 2011). In Salmonella, Vogel and colleagues assigned ∼22,000 RNase-E-mediated PSSs resulting from a 28°C to 44°C temperature shift, reinforcing the role of RNA processing and turnover in gene expression control (Chao et al. 2017). Our study also highlights the strong impact of RNA processing and stability on the bacterial adaptation response. Moreover, we improve on existing techniques through the dual mapping of TSSs and PSSs in a single RNA sequencing run. Combined measurement of changes in frequencies (TIF and PF) allow us to discriminate transcriptional and post-transcriptional events and correlate these fundamental events to variations of RNA levels. In the case of PSSs, changes in PF unveils cleavage sites that figure prominently in the evolving RNA landscape, and hence, in the bacterial response to environmental cues.

Conclusion

Here we describe, for the first time in the literature, the direct and simultaneous observation of changes in the efficacies of transcription initiation and RNA cleavage at the genomic scale. This work provides the first comprehensive and dynamic view of fundamental processes modulating gene expression.

The σE regulon in E. coli served as an ideal model by providing a list of 5′ RNA ends where frequencies were expected to vary upon σE induction: the σE-TSSs. From σE-TSSs, general patterns were inferred for activation and repression of transcription initiation at σE-independent TSSs. We revealed that a significant number of assigned TSSs showing changes in TIF displayed a biphasic pattern, conferring transcriptional oscillation to the bacterial population sensing and responding to an environmental signal. Within this same evolving transcriptome, we discovered three classes of PSSs including one whose activity responds to transcription. The factors providing selectivity to PSSs or able to couple transcription initiation to RNA cleavage remain to be established.

In the light of the data and the results, it is evident that additional experiments will be required to strengthen and refine conclusions on transcriptional and post-transcriptional events shaping the evolving RNA landscape. Yet, absent the insights provided by tagRNA-seq in this study, these avenues of investigation would have remained obscured. Directed mutagenesis on certain PSSs, or swapping promoters and cleaved RNA sequences, may clarify the coupling observed between transcription initiation and certain cleavage sites. Such a work will also further clarify the activity of PSSs as a function of RNA amounts or transcriptional activity.

Together, our data show that tagRNA-seq is suited for the direct and simultaneous analysis of transcriptional and post-transcriptional events taking place at sub-generation timeframes in response to an environmental cue. We have laid the foundation to probe, at the genomic scale, major processes controlling gene expression during physiological adaptations in bacteria. This generalized approach will provide dynamic and mechanistic insights to studies that exploit changes in RNA levels to decipher and understand gene expression reprogramming or for engineering purposes.

MATERIALS AND METHODS

Bacterial growth and RNA preparation

Three individual clones of the E. coli strain MG1655/pZE21-rpoE (Bury-Mone et al. 2009) were grown 18 h in LB medium at 37°C with agitation (200 rpm). Each culture was diluted 1/1000 in prewarmed LB medium and grown to an OD600 ∼ 0.3. Then, the t0 sample was collected and aTc was immediately added to a final concentration of 10 nM. Samples of each culture were collected at 5, 10, 15, and 20 min after the addition of aTc. The sampling period was deliberately shorter than the doubling time of the bacterial strain and total RNA was prepared as previously described (Bury-Mone et al. 2009).

RNA tagging and sequencing

RNA tagging was performed as previously described (Fouquier d'Herouel et al. 2011; Innocenti et al. 2015). Briefly for each sample, 20 µg of total RNA was ligated with the RNA adaptor PSS-tag (5′-GCAUAGGGGUAAA-3′) using the T4 RNA ligase I (New England Biolabs). Samples were then treated with the tobacco alkaline phosphatase (TAP; EPICENTRE Biotechnologies) and ligated to the TSS-tag RNA adaptor (5′-GCGAGACUGAGAA-3′). Since 5′ RNA ends can be tagged by both RNA adaptors and may provide ambiguities in assignments (5′ termini coined “UNDs”) (Innocenti et al. 2015), we performed at t0 for each of the three experimental series, a transcriptome where the TAP treatment was omitted. The comparison between samples treated with TAP and untreated enabled us to assign TSSs otherwise classified as “UND” (Supplemental Section S2). The 18 RNA samples (five time points for each of the three kinetic series, and the three untreated TAP at t0) were sequenced on a SOLiD 5500W technology instrument (Thermo Fisher Scientific).

Alignment, coverage, gene expression level, 5′ tag detection, and normalization

We applied an updated version of our previously described procedure (Innocenti et al. 2015). 5′ and 3′ sequencing adapters (Supplemental Table S6) were first stripped in colorspace using cutadapt v1.16 (Martin 2011), retaining sequences of ≥30 nt (options -c -m 30 -g file:Solid_5prime.fasta -a file:Solid_3prime.fasta). Remaining sequences were detagged in colorspace of TSS- and PSS-tags using “cutadapt,” retaining sequences ≥20 nt (option -c -m 20), resulting in reads sorted as TSS, PSS, or unknown. These reads were then aligned to the genome of the E. coli K12 strain MG1655 (GenBank: U00096.3) using Bowtie v1.2.2 (Langmead et al. 2009). Alignments were used for calculating expression levels with Cuffdiff v2.2.1 (Trapnell et al. 2010). Ribosomal RNA regions were masked and appropriate sequence corrections were applied (options -M U00096_rRNA.gtf –-multi-read-correct –library-type fr-secondstrand –frag-bias-correct). The remaining analysis was done as in Innocenti et al. (2015). Normalized data are available in numerical format in Supplemental Table S1 (RNA levels) and Supplemental Table S5 (tag counts).

Assignments of transcription starts and RNA processing sites

A 5′ RNA end was considered mapped when at least 42 tags in total were found in the 18 transcriptomes. As RNAP can initiate transcription at a few consecutive nucleotides within a unique promoter (Robb et al. 2013; Vvedenskaya et al. 2015), each 5′ RNA end mapped was treated individually. Determination of the nature of 5′ RNA ends, assigned as “TSS,” “PSS,” or “UND,” was based on tag-counts and TAP treatments. Further details are provided in Supplemental Section S2.

Changes in transcription initiation and RNA processing frequencies

Changes in the frequency of processes generating 5′ RNA ends was postulated as null before the induction of σE at t0. Tag-counts reflect the relative abundance of a given 5′ RNA end (Innocenti et al. 2015). The tagging rate, Ri, is the number of tags appearing or disappearing per time unit and reflects the frequency of the process generating a given 5′ RNA end. Changes in the frequency was inferred by the average variation of Ri between two consecutive time points, a magnitude that we termed “tagging acceleration” (Ai).

Statistical analysis

Reported TSS and PSS numbers gave rise to corresponding empirical cumulative distribution functions, yielding empirical P-values for the number of sites detected at each genomic position. Rate P-values were obtained from Fisher's combined probability test by establishing the test statistic which then was evaluated against the χ2 distribution with 4 degrees of freedom. nɛ{1,2,3,4} designates each time point following the initial one, that is, 5, 10, 15, and 20 min, and p0 = 0. Acceleration P-values were similarly calculated by aggregating rate P-values. Results with P-values <1.45 × 10−5 were considered significant (Bonferroni correction).

Consistency between kinetic series

Log-ratios of total tag-counts between consecutive time points were calculated for each reported position within each kinetic series. Observations at a given position and time point were called “point-wise consistent” if ratios from all three kinetic series showed the same trend (i.e., were all positive or zero and negative or zero, respectively, with at least one non-zero observation). We allowed for a pseudocount of one in the calculation of total tag log-ratios to avoid numerical divergence without introducing bias. Observations at positions with all four time-points reproducible were called “globally consistent.” Thus, for all 3447 mapped ends, we analyzed the consistency of tag-counts between kinetic series as provided in Supplemental Table S4.

DATA DEPOSITION

The data are available in the NCBI SRA repository (PRJNA561076).

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.