Elucidation of Regulatory Modes for Five Two-Component Systems in Escherichia coli Reveals Novel Relationships

BaeR and CpxR are activated by a wide range of stressors (10–14) and play a regulatory role in the envelope stress response (ESR) system (10, 15–17). We stimulated the wild-type, ΔbaeR, and ΔcpxR strains with 5% ethanol stress, which has been identified as a particularly effective inducer (18) to elucidate their regulons.

(i) BaeR. BaeR is involved in multidrug resistance through regulation of the MdtABC efflux pump (19–21). However, BaeR has also been shown to impact genes related to many other functions, such as signal transduction, chemotactic responses, flagellar biosynthesis, and maltose transport (22). Despite this range of functionality, there are only eight genes in the previously published BaeR regulon. The BaeR iModulon exhibited 37.5% (3 out of 8) overlap with these published regulon genes, consisting of spy, baeR, and mdtA genes (Table 2). Five of the published regulon genes (mdtB, mdtC, mdtD/iceT, baeS, and acrD) were not identified in the BaeR iModulon. This discrepancy was most likely due to differences in BaeR induction conditions across the published results (10–13) or the possibility that the three iModulon genes are much more significantly regulated by BaeR than the other five genes in the published regulon. To test this claim, ChIP-exo results were analyzed (Data Set S2) and then compared to the iModulon. ChIP-exo results indicated 21 BaeR binding peaks upstream of 17 operons consisting of 32 genes. We identified binding sites upstream of spy and mdtA genes (mdtABCD-baeSR operon). Since a peak was identified upstream of mdtA gene, which is part of the mdtABCD-baeSR operon (see Materials and Methods), we considered all the genes in that operon to be under the direct regulatory network of BaeR. In total, we designated six genes as functional direct targets of BaeR. No binding peaks were identified upstream of acrD, indicating that it is indirectly regulated by BaeR.

Differentially expressed genes were analyzed to further expand the global regulatory network of BaeR. A total of 328 genes were found to be differentially expressed due to the baeR knockout under 5% ethanol stress condition (Data Set S3). ChIP-exo peaks were detected upstream of nine of these genes, expanding the number of direct targets from seven to nine. Of these nine targets, two genes, intF and tnaA, have not been previously identified as a part of the BaeR regulon. However, tnaA has previously been shown to contribute to BaeR activity under indole stress (23), as it is involved in metabolizing tryptophan to indole. The remaining 319 differentially expressed genes could either be indirect targets of BaeR (Fig. 1 and Data Set S4) or other transcription factors responsive to ethanol such as CpxR. Functional analysis of indirect targets through COG (cluster of orthologous group) categories confirmed BaeR’s manifold activity (Fig. 2).

Overall, 2.74% of DEGs were under the direct regulatory network of BaeR. The remaining 23 genes found through ChIP-exo results could be classified as nonfunctional targets or cases of spurious binding (Fig. 1 and Data Set S4). Binding of these genes may have not led to expression if the induction period or strength was insufficient; alternatively, these genes might have shown ChIP-exo peaks as a result of off-target binding.

(ii) CpxR. CpxR is among the most extensively studied of the response regulators, which is unsurprising when considering its host of functions. CpxR has been shown to have a regulatory role in the envelope stress response system (10, 15–17), protein folding and degradation (24–27), pilus assembly and expression (28–30), secretion (15), motility and chemotaxis (15, 31), biofilm development (14, 32), adherence (33, 34), multidrug resistance and efflux (12, 35), porins (36), and copper response (37, 38), among others. Two iModulons were linked to the CpxR regulon: CpxR and CpxR-KO (KO stands for knockout). The CpxR iModulon consists of genes that were regulated specifically under wild-type conditions, and CpxR_KO consists of genes differentially expressed due to CpxR deletion. Collectively, both CpxR iModulons consist of 16 genes that cover diverse functions such as motility, inorganic ion transport and metabolism, and carbohydrate metabolism, among others. Nine of the 16 iModulon genes were previously published regulon genes, which is only a fraction (14%) of CpxR’s previously identified 66 genes (9) (Table 2). This provides extra incentive to eventually validate the CpxR iModulon with different induction conditions. Interestingly, the CpxR iModulon consisted of seven new genes that were not previously identified to be part of the CpxR regulon (Table 2). These seven new genes may indicate expansion of the CpxR regulon specific to ethanol stress. To validate this hypothesis, ChIP-exo peaks were examined. We were able to identify 44 binding peaks upstream of 16 operons (consisting of 30 genes). Of these genes, five were in the CpxR iModulon (cpxP, cpxR, raiA, yagU, and yccA). Two genes (raiA and yagU) in the CpxR iModulon with associated binding peaks had not been previously identified as part of the CpxR regulon. The five overlapping ChIP-exo/iModulon/regulon genes were classified as direct targets of CpxR under ethanol. For the remaining five iModulon genes that are not previously published regulon genes and do not have binding peaks identified, 60% are part of the uncharacterized “y-ome” (39). These genes may be the result of cross-regulation between different pathways responding to ethanol and put under the category of cross-regulatory network.

We inspected DEGs to further understand the nature of the five new iModulon genes and to expand the global regulatory network of CpxR. Differential gene expression analysis expanded the indirect targets of CpxR to include 368 genes. Comparison of ChIP-exo binding peaks and differentially expressed genes of the cpxR knockout versus the wild type under 5% ethanol stress condition revealed seven additional genes (cpxA, fimA, fimI, lldD, lldR, ppiA, and yceI) that fall under the direct regulatory network of CpxR, expanding this network to 12 genes. The remaining 357 DEGs formed the indirect regulon of CpxR, including eight genes from the previously published regulon (Fig. 2 and Data Set S4). A diverse group of COG categories corroborated CpxR’s regulatory role.

Interestingly, among the 368 differentially expressed genes, 194 were differentially expressed in both knockout strains (baeR and cpxR), suggesting nonspecific targets of BaeR and CpxR and high cross-regulation between these two pathways. The remaining 17 ChIP-exo binding peaks overlapped with one previously published regulon gene (ompC) which was not part of either the DEGs or iModulons. This suggested that ompC was neither directly nor indirectly affected by cpxR knockout under ethanol stress, and therefore is considered to be part of the “hypothetically functional binding” network (Data Set S4). Five iModulon genes (alx, mdtJ, yjfN, yliI, and yobB) did not form part of direct and indirect targets but instead formed part of the network that is involved in cross-regulation. The remaining 52 previously published regulon genes that did not have binding sites upstream and were not differentially expressed also formed part of a cross-regulatory network (Data Set S4).

Decoding the role of metal and nutrient sensors in the transcriptional regulatory network may provide a deeper understanding of how bacteria utilize these systems to sense nutrients and ionic strength of the environment. These mechanisms play a crucial role in the organism’s ability to adapt to various environmental niches. The KdpE, PhoB, and ZraR two-component systems were selected for further analysis in this category.

(i) KdpE. The KdpE response regulator’s sole target is the kdpFABC operon, which has been shown to activate the potassium uptake system during low potassium conditions (40, 41) or salt stress (42, 43). To induce KdpE, a small amount of potassium (0.1 mM KCl) was added to the cells in Tris-maleic acid minimal medium (TMA), which aligns with previous experiments where successful induction was achieved (44, 45). The KdpE iModulon included the expected kdpA, kdpB, and kdpC potassium uptake genes; however, the iModulon does not include the fourth gene in the kdpFABC operon, kdpF (Table 2). This result provides validation of previous studies that have instead described the target operon as kdpABC (43, 46). kdpF may function only as a stabilizer of the transporter complex (47), and the gene may also just be a less vital component of the regulon. Apart from the kdpABC operon, the only other KdpE iModulon genes include those that code for the TCS itself, kdpD and kdpE. The small KdpE iModulon confirms the existing knowledge of the response regulator as targeting a single locus in the bacterial genome. Despite the small iModulon size, there were still 31 ChIP-exo binding peaks upstream of 37 genes, including kdpABC genes from KdpE iModulon observed. Surprisingly, no binding peaks were observed upstream of two iModulon genes, kdpDE; however, kdpDE is downstream of the kdpFABC operon, which may indicate that kdpDE is in the same transcription unit as kdpABC. Consistent with this, when comparing 197 DEGs in response to kdpE knockout under potassium limited condition (0.1 mM KCl), we did find kdpDE to be differentially expressed. Comparing ChIP-exo and DEG results further added three extra genes (adeP, uxaA, and uxaC) to the direct regulatory network of KdpE, amounting to six genes under direct regulation of KdpE. The remaining 191 DEGs can be interpreted as potassium-dependent genes under indirect KdpE regulation. As expected, the COG category “inorganic ion transport” was overrepresented. In addition, the categories “energy production and conversion” and “metabolism and nucleotide transport and metabolism” were also overrepresented. A wide array of COG categories suggests that KdpDE may indirectly regulate genes involved in metal efflux pump and electron transport systems to maintain homeostasis (Fig. 3).

Notably, there was evidence of kdpF being down-expressed in the KdpE knockout (log fold change of −3.22) but due to an insignificant P value, this was not regarded as differentially expressed. Therefore, kdpF appears to be the sole gene from the KdpE published regulon that was not identified by DEG analysis in this study (placing it in the category of “cross-regulated” genes). Apart from kdpF and genes categorized as direct or indirect targets for KdpDE, there were also 31 genes identified with ChIP-exo that did not overlap with other data sets and were thus categorized as potential spurious binding genes.

(ii) PhoB. The PhoB response regulator has been demonstrated to regulate phosphate uptake and metabolism under phosphate-limiting conditions (48–51) and is linked to virulence in pathogenic E. coli (4). To induce PhoB via phosphate deficiency, cells were first grown in M9 minimal medium and then washed with and further grown in M9 minimal medium without phosphate (M9-P). The PhoB iModulon consisted of the expected pstSCAB operon for regulation of the phosphate uptake system (52–54) and the genes coding for PhoRB itself. In contrast to the modest six-gene iModulon, the previously published PhoB regulon is considerably larger and is thought to consist of 60 genes (55, 56), including all six of the iModulon genes. For the binding site analysis, we found 12 ChIP-exo binding peaks upstream of 17 genes. Six of these genes were included in the iModulon (the pstSCAB phoBR operons), indicating direct response to phoB knockout. When comparing the 113 DEGs under phosphate-limited condition to the ChIP-exo binding peaks, an additional gene, phoU (encoding a phosphate-specific transporter) (which was also part of the previously published regulon) was added to the direct regulatory network of PhoB. The remaining 106 DEGs can be identified as phosphate-dependent genes, including six previously determined regulon genes (cusB, phnC, phnD, phnN, phoA, and ugpB) or indirect targets of PhoB. Interestingly, the COG category “cell motility” was more pronounced than “inorganic ion transport and metabolism,” indicative of the involvement of PhoRB in cell migration in addition to adapting to a phosphate-limited environment. Aside from the direct and indirect targets, the remaining 57 genes found to be associated with PhoB included 10 genes that were identified as a potential result of spurious binding, two genes that fell under the hypothetically functional binding category, and 45 genes that were possibly related to the KdpE network via cross-regulation.

(iii) ZraR. Finally, we chose to examine ZraR (formerly called HydG), which had been previously characterized as a σ⁵⁴-dependent response regulator that activates chaperone ZraP to provide tolerance to high zinc concentrations (57–59). More recently, it has been proposed that ZraR is not directly involved in zinc resistance (60) but rather is simply activated by zinc to achieve its broader role in the envelope stress response, similar to CpxAR (61). Since the commonality of these differing theories centers on zinc, we chose to study ZraR regulatory activity on LB medium supplemented with 1 mM ZnCl₂ (57). The ZraR iModulon included only two of the three previously published regulon genes, zraR and zraP. One of the regulon genes, zraS, was not identified by the iModulon. However, ChIP-exo binding peaks were observed upstream of all three regulon genes (zraSR and zraP), confirming that these genes are under the direct control of ZraR. The ChIP-exo results indicated only one additional direct target gene, mgtA, and the remaining 12 binding peaks might be attributed to spurious binding. DEGs were examined to identify indirect targets and expand direct targets (if any) of ZraR. In the presence of 1 mM ZnCl₂, 97 genes were found to be differentially expressed between the zraR knockout and wild-type strain. Four of these genes (zraSR, zraP, and mgtA) had corresponding upstream binding peaks found by ChIP-exo or were included in the ZraR iModulon. The remaining 93 DEGs can be attributed to indirect ZraR regulation. Therefore, ZraSR seems to be involved in cell migration and carbohydrate metabolism (Fig. 3) similar to CpxAR, further strengthening the cross-regulation of these two TCSs.

Overall, integrating iModulons with ChIP-exo and DEGs enabled us to expand the global regulatory network contribution of each selected TCS and categorize genes in different modes (direct or indirect) of transcription regulation.

To identify potential novel inducers and interconnection between TCS pathways, we compared the activity levels of each TCS iModulon across all PRECISE experiments. To this end, the iModulon activities of each TCS were analyzed across their expected induction condition(s) in our experiments, as well as the previous conditions of PRECISE experiments. By calculating each TCS iModulon activity for the experiments under control conditions (e.g., LB medium for BaeSR, CpxAR, and ZraSR), a baseline was established for expected TCS behavior in the absence of an induction condition. From there, iModulon activities were calculated for all other experimental conditions across PRECISE and the other TCS experiments of this study, and then these activities were evaluated against the control condition iModulon activity. While many experimental conditions yielded iModulon activities indiscernible from those of the control, several key environmental inducers were shown to have significant differential activity for each TCS (Fig. 4). These results are important for understanding which types of environments are most likely to trigger TCS action. The ability to succinctly examine the activity levels of an iModulon over an ever-growing compendium of conditions enables the broad characterization of the iModulon’s response (or lack of response) to each of the conditions. Figure S2 shows the activity profiles for each of the five TCS response regulators (RRs) across all conditions, including PRECISE conditions and the new conditions tested in this study.

When observing the activity levels of ethanol sensors (BaeR and CpxR), we noted that a high average activity level occurred when the wild-type strain was cultured on LB plus 5% ethanol. Interestingly, both the BaeR and CpxR iModulons exhibited increased activity for a previous set of experiments that studied the effect of osmotic stress on E. coli K-12 MG1655 in relation to the OmpR regulon (62). Although Bury-Moné et al. (18) had suggested that 0.6 M NaCl did not activate a CpxR and BaeR target gene reporter, our results suggest that osmotic stress (via the addition of 0.3 M NaCl) may induce BaeR and CpxR to some extent. A prior study attempted to identify cross-regulation between CpxAR and OmpR-EnvZ systems (63). In addition to osmotic stress induction, previous PRECISE iron starvation conditions (64) also yielded a significant increase in BaeR iModulon activity, which indicates that BaeR may also be involved in response to iron deficiency (Fig. 4). For CpxR, iModulon activity peaks also occurred across experiments testing for phosphate starvation. The array of conditions inducing the BaeR and CpxR iModulons is reasonable considering the response regulators’ extensive functionalities and roles in responding to cell envelope stress (Fig. 4).

KdpE was shown to be directly phosphorylated by PhoR under simultaneous potassium and phosphate limitation (45). Consistent with this, we observed high peak activity of the KdpE iModulon for our TCS experiments conducted on phosphate-depleted M9 medium, in addition to the expected induction of the wild-type strain on TMA plus 0.1 mM KCl (Fig. 4). Another notable iModulon activity increase was observed for a subset of PRECISE experiments that were conducted under osmotic stress (62), which is consistent with KdpE’s documented induction under salt stress (42, 43). Compared to KdpE, the PhoB iModulon was active in several of our other TCS experimental conditions, with an especially notable trend of activity increases for all our samples on LB plus 5% ethanol stress. Although PhoB is not generally thought to be a component of the envelope stress response, these iModulon activity increases could potentially be linked to previous research that has found ethanol to induce calcium phosphate crystallization (65, 66) that could limit the amount of freely available phosphate and mirror starvation conditions. For the PRECISE experiments, the PhoB iModulon activity was significantly lower than the PhoB iModulon activity under the phosphate starvation condition in this study. For example, the experiment with the highest PhoB iModulon activity in PRECISE, a yddM knockout experiment (67), had a PhoB activity of 2.8, which is only 23% of the phosphate starvation condition PhoB iModulon activity of 11.8. The lower PhoB iModulon activities in the PRECISE experiments is consistent with the absence of phosphate starvation conditions in any of the previous experiments (Fig. 4).

Notably, the ZraR iModulon includes only the two expected regulon genes, zraR and zraP, yet this iModulon exhibits significantly heightened activity over a surprisingly wide range of experimental conditions apart from LB plus 1 mM ZnCl₂ (Fig. S2). For this study’s experiments, the ZraR iModulon activity was high for all experiments conducted in phosphate-depleted M9 medium and most experiments conducted in LB plus 5% ethanol. The high activity of the iModulon across diverse conditions supports the view of ZraR as a broader regulator past simple zinc resistance; specifically, ethanol’s activation of the iModulon could be indicative of ZraR’s relation to ESR, as previously suggested (61) (Fig. 4).

Overall, combining our new data set with PRECISE and applying ICA enabled us to identify potential cross-regulation between TCS pathways, which allowed the identification of novel or noncognate TCS-inducer pairs. Both the known induction conditions validated by this study and the novel induction conditions illuminated by this study are summarized in Table 3 below.

Before using the various gene data sets to scrutinize the regulatory network of each TCS, the quality of the newly generated data was first evaluated. Namely, the peaks found with ChIP-exo were shown to contain the published consensus motif of each TCS, and the RNA-seq data were analyzed to verify TCS regulon induction under their respective stimulation conditions.

To validate the binding peaks found by the ChIP-exo analysis, the peaks associated with direct targets from this study were compared to known binding peak motifs using AME (analysis of motif enrichment) from the MEME Suite (68). AME computes position-weight matrix (PWM) score by comparing sequence to the motif. Strong binding sites may have higher PWM scores, while weak binding sites have low PWM scores. The experimentally identified peaks that were found to match the consensus motif as represented by RegulonDB (69) are summarized below in Table 4.

All five of the examined TCSs were shown to have at least one direct target ChIP-exo peak that matches the consensus motif using AME. In accordance with expectations, the matching peaks pinpointed by AME included canonical primary targets of each TCS response regulator. A visualization of the peak genomic location and sequence match to the consensus motif is shown below in Fig. 5 for the pstSCAB-phoU peak (identifier [id], peak_7) that was identified as a match to the PhoB consensus motif.

The results of this study rely upon the adequate stimulation of the various TCSs under their various induction conditions. As a verification of these induction conditions, a closer look is given to differential gene expression in both the unstimulated and stimulated conditions. For the histidine kinase and response regulator genes of each TCS, the response regulator knockout versus wild-type RNA-seq gene expression is compared under both the unstimulated and stimulated conditions to confirm successful induction (Table 5). A broader comparison of the response regulator knockout versus wild-type expression for the entire direct regulon of each TCS is available in Data Set S5.

Both the histidine kinase and response regulator genes were shown to be differentially expressed at their respective stimulation condition with a log fold change greater than 1.5 and a P value less than 0.05 (Table 5). This observation confirms that each of the TCS genes was more expressed in the wild-type strains versus the response regulator knockout (KO) strains under induction. All histidine kinase and response regulator genes other than cpxR exhibited a greater magnitude of log fold change than the unstimulated condition. Overall, this serves as an independent validation of successful induction of the histidine kinase elements to the respective stimuli.

Strains used in this study are E. coli K-12 MG1655 and its derivatives (deletion strains and myc-tagged strains). As previously described (75), strains retaining 8-myc were generated by a λ red-mediated site-specific recombination system targeting the C-terminal region of each selected response regulator. Knockout mutants (ΔbaeR, ΔcpxR, ΔkdpE, ΔphoB, and ΔzraR) were generated according to the procedure of Datsenko and Wanner (76), using pKD13 and pKD46 as suggested by the authors (please refer to Table S1 in the supplemental material for a full list of oligonucleotides). All knockouts are complete deletions, except for KdpE; the KdpE locus could potentially produce an 18-bp peptide, since the sequence after the stop codon contains too many A/T-rich repeats for proper oligonucleotide annealing so the binding site had to be shifted into the coding region (please see Table S1 for the oligonucleotide list). If not stated otherwise, cells were grown in precultures prepared from glycerol stocks overnight and then transferred to fresh medium the next morning. For ethanol sensors (BaeR and CpxR), cells were grown at 37°C in liquid LB medium to an optical density at 600 nm (OD₆₀₀) of 0.5, and then ethanol was added to a final concentration of 5% (wt/vol). To ensure proper binding of RRs to their target genes, cells were grown for 30 min in the presence of ethanol before being collected for ChIP-exo and RNA-seq. This growth condition was chosen according to Bury-Moné et al. (18) (see Table 1 in reference 18) who showed that 5% ethanol led to the strongest induction of β-galactosidase reporters for CpxR and BaeR among all tested conditions (i.e., 2 mM and 4 mM indole, 3% and 5% ethanol, 0.5 mM dibucaine, 5 mM EDTA, and 0.6 M NaCl). For RNA-seq controls, cells were grown in liquid LB medium and collected at an OD₆₀₀ of 0.5. The KdpDE TCS induces expression of the high-affinity K⁺ transporter KdpFABC under K-limiting conditions (Schramke et al. [45]). According to Schramke et al. (45), KdpD acts as a phosphatase on KdpE-P and prevents production of the high-affinity K⁺ transporter at high extracellular K⁺ concentration (>5 mM). When environmental levels of K⁺ fall below the threshold for autokinase activation, kdpFABC expression is initiated; however, as long as the intracellular K⁺ concentration remains high, the KdpD phosphatase activity remains stimulated (Schramke et al. [45]). Therefore, to induce KdpDE, cells were grown in K-sufficient conditions, washed, and then grown under K-limiting conditions to ensure drops in intracellular K⁺ levels below the threshold level and thus proper induction of the TCS. To do this, cells were grown at 37°C overnight in liquid Tris-maleic acid minimal medium (TMA) (44, 45) supplemented with 115 mM KCl and 0.4% (wt/vol) glucose and then washed twice with TMA containing 0.1 mM KCl and 0.4% (wt/vol) glucose. Cells were inoculated in TMA with 0.1 mM KCl and 0.4% (wt/vol) glucose and collected at an OD₆₀₀ of 0.5 for ChIP-exo and RNA-seq analysis. For RNA-seq controls, cells were grown in liquid TMA supplemented with 115 mM KCl and 0.4% (wt/vol) glucose and collected at an OD₆₀₀ of 0.5. For induction of the PhoRB TCS, we chose to use M9 minimal medium without any phosphate to ensure proper induction of the system. To our knowledge, this condition was not used before. To induce phosphate-limiting conditions, cells were grown in liquid M9 minimal medium (containing phosphate) until an OD₆₀₀ of 0.5 was reached. Growing cells until the mid-exponential growth phase was necessary, because the lack of phosphate led to growth arrest. The cells were then washed three times with M9 minimal medium without phosphate (M9-P) (without Na₂HPO₄ and KH₂PO₄) and incubated in M9-P for 60 min (about the doubling time of strain MG1655 in M9 minimal medium) at 37°C. Cells were then collected for ChIP-exo and RNA-seq analysis. For RNA-seq controls, ΔphoB cells were grown in liquid M9 minimal medium and collected at an OD₆₀₀ of 0.5 (WT samples were taken from reference 7; see Table 1). For ZraSR, cells were grown at 37°C in liquid LB medium containing 1 mM ZnCl₂ by the method of Leonhartsberger et al. (57) but without the addition of glucose to the LB medium. At an OD₆₀₀ of 0.5, cells were collected for ChIP-exo and RNA-seq analysis. For RNA-seq controls, cells were grown in liquid LB medium and collected at an OD₆₀₀ of 0.5. Samples for ChIP-exo and RNA-seq were taken independently.

For RNA-seq, 3 ml of culture was mixed with 6 ml of RNAprotect bacterial reagent (Qiagen) and processed according to the manufacturer’s instructions. Cell pellets were frozen and stored at –80°C until processed. RNA was extracted using the Zymo Research Quick RNA fungal/bacterial microprep kit, according to the manufacturer’s instructions. rRNA was removed from total RNA preparations using RNase H. First, traces of genomic DNA were removed with a DNase I treatment. Then, secondary structures in the rRNA were removed by heating to 90°C for 1 s. A set of 32-mer DNA oligonucleotide probes complementary to 5S, 16S, and 23S rRNA subunits and spaced nine bases apart were then annealed at 65°C followed by digestion with Hybridase (Lucigen), a thermostable RNase H. Hybridase was added at 65°C, the reaction mixture was incubated for 20 min at that temperature, then heated again to 90°C for 1 s to remove remaining secondary structures, and finally returned to 65°C for 10 min. The reaction was quickly quenched by the addition of guanidine thiocyanate while still at 65°C before purifying the mRNA with a Zymo Research RNA Clean and Concentrator kit using their 200-nucleotide (nt) cutoff protocol. Carryover oligonucleotides were removed with a DNase I digestion which was started at room temperature and gradually increased to 42°C over a half hour. This was followed up with another column purification as stated above.

Paired-end library preparation was done using the KAPA RNA HyperPrep kit following the manufacturer’s instructions with an average insert size of 300 bp. The libraries were then analyzed on an Agilent Bioanalyzer DNA 1000 chip (Agilent). RNA-seq libraries were sequenced on the NextSeq 500 (Illumina) at the Salk Institute for Biological Studies, La Jolla, CA, USA. RNA-seq experiments were performed in duplicates.

To identify binding sites for each response regulator under their respective induction conditions (see above and Table 1), ChIP-exo experiments were performed as described previously (64). Specifically, cells were cross-linked in 1% formaldehyde and DNA bound to each response regulator was isolated by chromatin immunoprecipitation (ChIP) with the specific antibodies that recognize the myc tag (9E10; Santa Cruz Biotechnology) and Dynabeads Pan Mouse IgG magnetic beads (Invitrogen). This was followed by stringent washing as described earlier (77). To perform on-bead enzymatic reactions of the ChIP-exo experiments, ChIP materials (chromatin-beads) were used (5). The method has been described in detail previously (64, 67). ChIP-exo experiments were performed in duplicates. ChIP-exo libraries were sequenced on the NextSeq 500 (Illumina) at the Sanford Burnham Prebys Medical Discovery Institute, La Jolla, CA, USA, or on the NextSeq 500 (Illumina) at the Salk Institute for Biological Studies.

The ChIP-exo sequence reads were mapped onto the reference genome (E. coli MG1655 GenBank accession no. NC_000913.3) using bowtie (78) with default options as also described previously (64, 67). This generated SAM (sequence alignment/map) output files. Peak calling from biological duplicates for each experimental condition was done using the MACE program (79). Two-level filtering of ChIP-exo data was done: first, we removed any peaks with a signal-to-noise (S/N) ratio of less than 1.5 to reduce false-positive peaks. The top 5% of the signals at genomic positions was set as the noise level as previously described (64, 67). MetaScope (https://sites.google.com/view/systemskimlab/software) was used to visualize the peaks. To further reduce the false-positive peaks, a cutoff analysis was done for each regulator on the S/N ratio (see Fig. S1 in the supplemental material). We removed any peaks that had a S/N ratio below the threshold. Thereafter, the peak was assigned to the nearest operon on both strands. Operons located within 500 bp of the peak were considered. Specifically, if a peak was upstream of an operon, we considered all the genes in that operon to be regulated by that response regulator.

After the final S/N thresholding and 500-bp cutoffs were applied to filter the ChIP-exo data, the resulting binding peaks were buffered with a margin of 20 bp on either end of their sequence and then compared to the published consensus motif to identify matching sequences (Table 4). The search was performed by analysis of motif enrichment (AME) from the MEME suite (68) using the options “--scoring avg --method fisher --hit-lo-fraction 0.25 --evalue-report-threshold 10.0 --control --shuffle-- --kmer 2.” The “--scoring avg” option specifies that the method of scoring a sequence compared to the motif’s PWN consists of averaging the motif odds score across all positions in the sequence. The “--method fisher” option selects the one-tailed Fisher’s exact test method to evaluate motif enrichment. The fraction of maximum log odds for scoring a potential match was specified as 0.25 using “--hit-lo-fraction 0.25,” and the E value for reporting was arbitrarily set at 10 (no TCS had more than 10 identified sequences) using “-evalue-report-threshold 10.0.” A zero-order Markov sequence model was created by AME to normalize for biased distributions; additionally, AME created a control sequence by shuffling the letters in each of the input sequences while preserving 2-mers using “--control --shuffle-- --kmer 2.”

RNA-seq sequence reads were mapped onto the reference genome (E. coli MG1655 GenBank accession no. NC_000913.3) using bowtie 1.1.2 (78) with the following options “-X 1000 -n 2 − 3 3,” where X is maximum insert size, n is number of mismatches, and −3 3 denotes trimming of three base pairs at the 3′ end. Read count was performed using summarizeOverlaps from the R GenomicAlignments package, using options “mode=“IntersectionStrict,” singleEnd=FALSE, ignore.strand=FALSE, preprocess.reads=invertStrand” (80). Differential gene expression was identified by DESeq2 (81). Genes with log₂ fold change of >1.5 and a false discovery rate (FDR) value of <0.05 were considered differentially expressed genes. Genes with P values assigned “NA” (for not available) based on extreme count outlier detection were not considered potential DEGs.

To identify independent sources of gene expression control, independent component analysis (ICA) was applied to the combined data set that consisted of all 278 RNA-seq gene expression data from PRECISE (7) and 32 data sets from this study which included gene expression data of knockout mutants and wild-type strains under stimulated and unstimulated conditions. Genes with less than 10 fragments mapped per million reads across the entire data set were removed from the RNA-seq data sets to reduce noise. FastICA was performed 100 times with random seeds and a convergence tolerance of 10⁻⁶. For each iteration, we constrained the number of components to the number that generates 99% variance in principal-component analysis (PCA). We further used DBSCAN to cluster the resulting components to identify a set of robust independent components using a minimum cluster seed size of 50. We repeated this process five times, and only the components that occurred consistently in each run were taken. iModulons were extracted as previously described (7).

We checked for the significant genes in each iModulon and mapped it to the set of genes regulated by a specific regulator or transcription factor. This comparison was done using one-sided Fisher’s exact test (FDR < 10⁻⁵). The regulon/regulator with the lowest P value was given the name of the iModulon. The regulon/regulators that were not characterized by this method were further analyzed through Gene Ontology (GO) enrichment. In this case, significant genes in each iModulon were compared to genes in each GO term, and a P value was assigned using the one-sided Fisher’s exact test (FDR < 10⁻⁵). The GO term with the lowest P value was assigned as the name of the iModulon.

The ChIP-exo and RNA-seq data sets have been deposited into GEO with the accession number GSE143856. All other data are available in the supplemental material.