Gene Networks and Pathways Involved in Escherichia coli Response to Multiple Stressors

Stress response helps microorganisms survive extreme environmental conditions and host immunity, making them more virulent or drug resistant. Although both reductionist approaches investigating specific genes and systems approaches analyzing individual stress conditions are being used, less is known about gene networks involved in multiple stress responses. Here, using a systems biology approach, we mined hundreds of transcriptomic data sets for key genes and pathways involved in the tolerance of the model microorganism Escherichia coli to multiple stressors. Specifically, we investigated the E. coli K-12 MG1655 transcriptome under five stresses: heat, cold, oxidative stress, nitrosative stress, and antibiotic treatment. Overlaps of transcriptional changes between studies of each stress factor and between different stressors were determined: energy-requiring metabolic pathways, transport, and motility are typically downregulated to conserve energy, while genes related to survival, bona fide stress response, biofilm formation, and DNA repair are mainly upregulated. The transcription of 15 genes with uncharacterized functions is higher in response to multiple stressors, which suggests they may play pivotal roles in stress response. In conclusion, using rank normalization of transcriptomic data, we identified a set of E. coli stress response genes and pathways, which could be potential targets to overcome antibiotic tolerance or multidrug resistance.


Introduction
As a common human pathogen, Escherichia coli faces and survives numerous stresses, either in the environment or inside the host during infection. Bacterial stress response systems sense, respond, and enable bacterial adaptation to physical and chemical challenges, such as shifts in pH, temperature, antibiotics, nutrition limitation, and oxidative stress. In certain cases, an adaptive response toward one stress provides an acquired resistance to a second one, a phenomenon known as cross-protection. In E. coli, stress response systems are mainly centered around two major mechanisms: (i) protein degradation through different sigma factors and (ii) phosphorylation cascades, functioning through two-component regulatory systems. Regardless of the original mechanisms, activation of stress response systems results in a transcriptional modulation [1].
Gene expression reprogramming under different stress conditions has been explored in many studies [2][3][4]. Stress-induced response generates changes in efflux systems, virulence determinants [5], membrane integrity [6], motility, translation, and transcription [3]. It also induces biofilm formation and leads to the generation of more antibiotic-resistant bacteria [7]. Although there are shared features in all stress responses [1], several studies identified specific stress responses with their cognate protein machinery and tried to identify central stress response genes which might be linked to cross-protection [3].
The bacterial adaptive stress response is rapid and reproducible because it involves the induction of a common set of genes or pathways, which lead to a phenotypic plasticity [3]. Although genes involved in specific stress responses have been widely studied, underlying mechanisms in E. coli remain to be explored. Most studies describe the alteration of E. coli transcriptional profiles in response to specific stresses; however, many gaps still exist in our knowledge. For example, only a handful of studies explore stress response as a complex network, with central proteins controlling underlying mechanisms [8]. Network studies elucidate novel conclusions by summarizing big information into one picture [9]. In particular, a few comparative studies explored stress tolerance across multiple stress conditions for the commonly used E. coli strains [10].
In this work, we conducted a thorough investigation of transcriptomic data of the responses of E. coli K-12 MG1655 to five different stressors. Through an integrated systems biology approach, we identified common gene sets connecting and regulating important biochemical pathways and biological functions. Therefore, these gene sets can be considered part of the general stress response and candidates to improve tolerance towards multiple stressors (cross-protection) or factors promoting multidrug resistance. In addition to providing a comprehensive overview of the stress responses to understand the control of this complex system, this study highlights the genes and pathways that are common among different stressors and treatment with antibiotics. Therefore, pathways involved in bacterial stress response could be potential targets to overcome antibiotic tolerance or multidrug resistance.

Identification of Differentially Expressed Genes under Selected Stressors
In a given data set, we defined differentially expressed genes (DEGs) as those genes whose transcript levels are statistically significantly higher or lower at a given condition. We used the GEO-embedded tool, GEO2R, which relies on the Benjamini-Hochberg procedure with a false discovery rate > 0.05 and generates the top 250 DEGs in each sample. To further verify the validity of the tool results, we conducted a two-tailed unpaired t-test to compare the normalized log 2 expression value of treated (stress-induced) to untreated (control) samples. Genes with ratios with P values below 0.05 were considered as significantly differentially expressed.
To generate the heatmaps for each stress condition, we used the online heat mapper (http://www.heatmapper.ca/expression, accessed and used on 12 September 2021), which provided hierarchical clustering of top 250 DEGs based on the study/sample included through average linkage as a default setting. We subsequently combined the top 250 up-or downregulated genes in each sample for each stress condition to form a one list of putative significantly differentially expressed genes under each stress condition.

Identification of Genes That Are Differentially Expressed under Multiple Stressors
The Venn diagram online tool (http://bioinformatics.psb.ugent.be/webtools/Venn, accessed on 6 June 2021) was used for combining gene lists for each stressor to determine genes affected by more than one stress.

Generation of Protein-Protein Interaction Networks
The Search Tool for the Retrieval of Interacting Genes/Protein (STRING) database [12] was used for generating a protein-protein interaction network (PPINs) by integrating DEGs common to three or four stress conditions. The tool was used with medium (cutoff score: 0.4) confidence. Each of these networks was exported to Cytoscape 3.8.2. Highly connected nodes in the network representing major sub-clusters were identified using MCODE (Molecular Complex Detection) clustering algorithm [13].

Determination of Functional Categories from Gene Lists
To determine the functional categories of the shared DEG, we retrieved the Uniprot accession numbers of shared DEGs between at least four (as a restrictive threshold) or at least three stressors (as a permissive threshold), based on Venn diagram-based visualization. The corresponding list was used as the input list for the Database for Annotation, Visualization, and Integrated Discovery (DAVID) [14,15].

Details of Data Sets Used in Our Analysis
During host infection, E. coli encounter different sources of stress, such as low gastric pH, sudden neutralization and relatively high intestinal pH, detergent effect of bile salts, nutritional restrictions, and oxidative stress when it passes to the environment. At some points of their life span, they face one or more of these stresses simultaneously and are adapted to change their transcriptional makeup based on the stressor(s) they face.
Accordingly, we used a computational workflow ( Figure 1) to investigate the E. coli K-12 MG1655 transcriptome under five different stress conditions, with integrated transcriptional data from multiple data sets (Table 1), as detailed in the Methods section.

Differentially Expressed Genes (DEGs) under Selected Stressors
A typical challenge in such high-throughput studies is how to shortlist transcriptionally altered genes without missing biologically important ones. The reverse challenge is to be overwhelmed with data or to include noise that may mask the signal. As a first step to prioritizing and shortlisting genes, we only included statistically significant DEGs. We then established an arbitrary, rather inclusive, cutoff of the top/bottom 250 DEGs based on their fold change. Through the above-mentioned criteria, we aimed to identify key genes affected by each stressor and common between many stressors. The list of 250 top upregulated and downregulated genes per stressor in each sample are listed in Tables S1-S5.
First, we compared the top 250 DEGs under each stress condition reported by individual studies; some studies involved different samples which we included as well. We observed a noticeable variation among the profiles of DEGs between studies describing the response to each stressor, probably as a result of experimental and technical differences (time of exposure, inoculum size, use of different stressor subtype, e.g., the oxidative response to tellurite and peroxynitrite vs. H 2 O 2 exposure) (Figure 2A-D). When we combined the list of DEGs retrieved from all included studies (Tables S6 and S7)

Identification of Genes Affected by Multiple Stressors
Next, we examined the common DEGs when E. coli K-12 MG1655 were subjected to different stressors. The overlap between transcriptional responses to multiple stressors could provide better understanding of general stress response and suggest targets for decreased tolerance towards groups of stressors.
We found a significant overlap of up-and downregulated genes within the five stressors. Within our list of upregulated genes, we identified 683 genes that are common between antibiotic treatment and oxidative stress, 431 between antibiotic treatment and heat stress, 436 between antibiotic treatment and cold stress, and 120 between nitrosative and antibiotic treatment stress. Regarding downregulated genes, we identified 620 common genes between antibiotic treatment and oxidative stress, 466 between antibiotic treatment and heat stress, 430 between antibiotic treatment and cold stress, and 124 common downregulated genes between the nitrosative and antibiotic treatment stress ( Figure 3A,B). To this end, we sought to analyze the functional potential of the DEGs shared by different stressors to unveil the functions and pathways playing pivotal roles in stress response. Identification of key functional categories is particularly important for understanding E. coli physiology, mechanism of adaptation to different ecological niches, and pathogenicity.

Functional Categories Involved in E. coli Response to Multiple Stressors
Functional categorization, using DAVID, of the DEGs that are shared between at least four stressors (Tables 2 and 3) revealed that 26.5% of the identified upregulated genes are involved in survival, general stress response pathways, and response to DNA damage and various stimuli, 15.5% in transport, 12% in various metabolic processes, 12% in biofilm formation, 6 % DNA repair, 5% in transcription regulation, and 1% in motility ( Figure 4A). Intriguingly, the analysis delineated 15 genes with uncharacterized functions that are upregulated by at least four out of the five stressors (representing 18% of the total number of upregulated genes). In parallel, the downregulation of various metabolic pathways is a shared feature between the response to all included stressors and constituted 57.5% of the 127 gene expression changes, 22.8% involved in transport, 5.5% in motility, 4% in response to various stimuli, 3% in peptidoglycan biosynthetic processes, and 2.4% in transcription regulation of cellular processes ( Figure 4B). Three uncharacterized genes (out of 127), ymfI, ydiJ, and yedE, were found to be downregulated under all stressors except nitrosative stress. Of note, nitrosative stress had generally fewer represented functional categories than other stressors among transcriptionally altered genes (Tables 2 and 3).     The functional categories of upregulated genes common between at least three stressors are presented in Table 4.  Branched-chain amino acid and phenylalanine ABC transport system (LivJHMGF) in E. coli K-12 are downregulated. livM, livG, and livF are downregulated in four out of five stressors. livJ and livH are downregulated in three out of five stressors. artPIQMJ genes of E. coli encoding a periplasmic arginine transport system are also downregulated. Five genes involved in sugar transport are downregulated in all the five stress conditions vs. one upregulated gene. This is consistent with the observation of the decreased transcription of energy-requiring processes related to growth.
Ten genes involved in biofilm formation and regulation are upregulated, suggesting that biofilm formation can be a consequence of bacterial stress response.
The transcription of the flagellar protein-encoding gene, flgL, was increased under heat, cold, oxidative, and antibiotic stress; in contrast, seven flagellar genes were downregulated by the same stressors to conserve resources for cell survival in sustained stress.

Discussion
In this study, we used a systems biology approach by analyzing hundreds of transcriptomic data sets to identify key genes and pathways involved in the stress tolerance of E. coli K-12 MG1655, as a model for bacterial response to multiple stressors. Although genes involved in specific stress responses have been frequently studied in E. coli, the underlying mechanisms behind global stress response remain to be explored, as only a handful of studies explore multiple stress responses as a complex network, with central proteins controlling underlying mechanisms.
Our approach offers a thorough investigation of transcriptomic data in response to five different stressors: heat, cold, oxidative, and nitrosative stress, along with antibiotic treatment-providing a comprehensive overview of E. coli stress response. Our work is a step toward understanding the control of stress response, and toward identifying common gene sets connecting and regulating important biochemical pathways. Consequently, these pathways could be potential drug targets for overcoming antibiotic tolerance or multidrug resistance.
While looking for DEGs under five different stress types, we observed a substantial set of common DEGs. In most cases, all or most genes of an operon clustered together or showed the same expression profile; however, it was also common that some genes within the same operon were not detected in every stress, while their chromosomal neighbors were. This is expected with the use of inclusion/exclusion thresholds for DEGs (e.g., statistical significance thresholds, or a certain percentile rank to include). Among the other reasons are batch effects, inter-experiment variabilities, and the different sensitivities of different probes or primers used for transcript quantification. The variability between media types (e.g., rich LB vs. defined M9 or MOPS media) had different effects on sample clustering depending on the applied stressor. For example, in the heat stress category, samples from the GSE42675, GSE15534, and GSE20305 studies were clustered in the right panel of the heat map. The first two studies shared the same type of medium (M9) and time of exposure (10 min) with a small difference in the applied temperature. Samples from GSE15534 and GSE2030 studies shared the same time of exposure (10 min) and temperature (45 • C), but used different media types (M9 and MOPS, respectively).
In the oxidative stress data sets, samples from studies GSE20305 and GSE61736 were clustered together in the left panel of the heatmap, as they shared the source of oxidative stress (H 2 O 2 ) and the same culture medium (MOPS), while GSE56133 sample was clustered with them as a different culture medium was used (LB) with the same type of stress (H 2 O 2 ). Among the sets of upregulated genes, six genes were upregulated under all the five studied conditions: • osmB encodes an outer membrane lipoprotein, known to be upregulated by osmotic shock and stationary phase [16][17][18]. Although the function of osmB is not fully characterized [19], its transcriptional regulation was thoroughly studied and found to depend on the response regulator, RcsB, and its transcription was found to increase following the activation of the RcsCDB phosphorelay [20]. Our analysis indicated that osmY, osmF, and osmC were downregulated under three, two, and one stressor, respectively. • nrdH and nrdI: NrdH, a gluatredoxin-like protein, serves as the electron donor for NrdEF and is reduced by thioredoxin reductase [21]. nrdH is part of the nrdHIEF operon, while nrdEF encodes class Ib ribonucleotide reductase (RNR), catalyzing the reduction of ribonucleotides (NTP) to deoxyribonucleotides (dNTP), a key reaction for DNA synthesis, replication, and repair. nrdE encodes the alpha 2 subunit, harboring the site for nucleotide reduction, and nrfF encodes the beta 2 subunit, which contains the metallocofactor required for catalysis at alpha 1 subunit [22]. Our analysis indicated nrdE, nrdF to be upregulated under at least four stressors. NrdI is a flavodoxin that mediates the generation of the tyrosyl radical cofactor of NrdF. Oxidative stress was reported to induce nrdHIEF expression (up to 23.4-fold), suggesting that E. coli overexpresses various reductases and electron donors to increase the cell's free radical scavenging capacity to cope with oxidative stress [23]. Our analysis indicated that this operon was upregulated under oxidative stress, antibiotic treatment, cold, heat, and nitrosative stress. • yqjI (nfeR) is part of the yqjH-yqjI operon. yqjI encodes a DNA-binding transcription repressor that regulates the expression of the NADPH-dependent ferric siderophore reductase yqjH [24]. YqjH is required for iron homeostasis in E. coli and it is also part of the Fur regulon [25]. Our analysis indicated that yqjH was upregulated under antibiotic and heat stress. In aerobic conditions, iron is present as insoluble iron hydroxides; consequently, bacteria produce siderophores, which are extracellular ferric chelators, to mobilize iron [26]. Once ferri-siderophores complexes are transferred through membranes to the cytoplasm, they are either degraded by esterase or reduced by ferric siderophore reductases to release ferrous ions [27]. • yhcN encodes a putative periplasmic protein involved in response to hydrogen peroxide and acid stress. A knockout strain lacking yhcN was found to be more sensitive to hydrogen peroxide and was more able to form biofilm than its parental strain [28].
YhcN was first reported to be upregulated in response to cytoplasmic acid stress by Kannan, et al. [29]. • ycfJ is similar to the Proteus mirabilis umoD gene, involved in the flagellar synthesis and cell elongation [30]. Its expression was reported to be induced in biofilm formation, while its null mutant is deficient in biofilm formation [31].
Both yhcN and ycfJ are upregulated under all the five stressors included in our analysis, including antibiotics. As their roles remain elusive, this observation may help in exploring their molecular functions, as it suggests leads to gene function and heavy involvement in global stress response pathways.
Based on the above findings, we conclude that all five stressors upregulate the transcription of various stress-responsive proteins, such as OsmB, and genes related to various reductases and electron donors, such as the nrdHIEF operon. More pathways and functions have been identified by functional analysis of upregulated genes common to at least three stressors ( Figure 5A).
• purH encodes the bifunctional AICAR transformylase/IMP cyclohydrolase, which catalyzes the last two steps of the de novo purine biosynthetic pathway converting AICAR to IMP [50].
Two extra genes involved in the purine biosynthetic pathway, purK and purE, were upregulated under four stressors (heat, cold, nitrosative stress, and antibiotic treatment) but not under oxidative stress.
Additionally, two genes are indirectly related to purine metabolism as they encode the enzymes that catalyze the formation of 5,10-methylenetetrahydrofolate (5,10-mTHF), which acts as a methyl group donor in the biosynthesis of various cellular components, including purines and methionine: • glyA encodes a serine hydroxymethyltransferase, which catalyzes the conversion of serine to glycine through forming 5,10-mTHF. It is activated by MetR, and repressed by MetR and PurR [51]. • gcvT is part of gcvTHP operon. It encodes the T-protein in the glycine cleavage system, which catalyzes glycine degradation and the formation of 5,10-methylenetetrahydrofolate 5,10-mTHF. The other two genes encode GcvH, or H-protein, and GcvP, or P protein [52][53][54]. Expression of the glycine cleavage enzyme system is induced by glycine, activated by GcvA, and repressed by PurR and GcvA [55,56]. Both gcvH and gcvP in our analysis were downregulated under three and four stressors, respectively, indicating the many-sided roles of GcvTHP beyond amino acid metabolism. GcvA is also the activator of gcvB. Both were reported to be related to stress response. The hdeAB operon, which encodes chaperone-like functions, was reported to be activated by GcvB and repressed by GcvA. Both HdeA and HdeB protect periplasmic proteins from aggregation by acid stress [57]. gcvB encodes a small regulatory RNA. gcvB knockout mutant was found to be more sensitive to oxidative stress and accumulate more endogenous reactive oxygen species than wild type. The role of gcvB in oxidative stress was also found to be conferred by increasing OxyR expression [58]. These findings suggest that gcvB and gcvA products allow E. coli to survive in presence of both oxidative stress and low pH. Our findings expand their regulatory functions to cold, heat, and antibiotic treatment.
Among downregulated genes, four are involved in the pyrimidine biosynthetic and salvage pathway: • pyrC encodes a dihydroorotase, which catalyzes the third reaction in the de novo pyrimidine biosynthesis pathway. By looking for other genes involved in pyrimidine metabolism in our data, we found seven more genes to be downregulated: pyrB, pyrD, and pyrI were downregulated under four stressors while pyrE, pyrF, pyrG, and pyrL were downregulated under two stressors. • carA is part of the carAB operon. It encodes the amidotransferase component, CarA, while the CarB subunit is the synthetase component of the carbamoylphosphate synthetase, involved in L-arginine synthetic pathway, along with the de novo uridine monophosphate (UMP) biosynthetic pathway (part of pyrimidine biosynthesis). We found carB to be downregulated under four stressors. • codB and codA form the codBA operon. CodB is a cytosine permease, which brings cytosine into the cell [59], while CodA is a cytosine deaminase (CDA), which catalyzes cytosine deamination into uracil. It is one of the enzymes in the pyrimidine salvage pathway, permitting the cell to utilize cytosine for pyrimidine nucleotide synthesis [60,61].
The downregulation of purine and pyrimidine metabolism, which we repeatedly observed here, has previously been reported as a response to acid exposure (pH 3.6) in E. coli O26:H11, which induces the accumulation of certain metabolites and enzymes involved in purine metabolism [62]. Pyrimidine levels were reported to decrease to 59%, after 20 min of HOCl stress, relative to the reference levels of untreated cells [63]. The biosynthesis of nucleotides, such as cytidine and uridine, was decreased after the combined ultrasound and mild acidic electrolyzed water treatment or electrolyzed water treatment alone, as well as the pool concentration of pyrimidines in the planktonic E. coli cells. Simultaneously, ribose-5-phosphate, the precursor for nucleotide synthesis, was downregulated in planktonic E. coli cells, after electrolyzed water treatment, indicating a depression in nucleotide biosynthesis [64]. In addition, levels of adenosine monophosphate (AMP) and uridine monophosphate (UMP), which are products of nucleotide degradation, were elevated in pulsed light treated E. coli, indicating DNA damage as well as blockage of RNA synthesis [65].
In addition to nucleotide metabolism-related genes, two genes are part of electron transport chain complexes: • cyoB encodes subunit I of the cytochrome bo 3 complex. It is a part of the cyoABCDE operon. Our analysis indicated that cyoA, cyoC, cyoD, and cyoE were downregulated under two out of five stressors. • nuoC is part of the nuoABCEFGHIJKLMN operon, representing the 13 subunits of NADH ubiquinone oxidoreductase. Our analysis detected the downregulation of nuoE, nuoF, nuoH, nuoI, and nuoJ under four stressors, and of nuoB, nuoG, nuoK, nuoL, nuoM, and nuoN under three stressors, while nuoA was detectably downregulated under two stressors.
E. coli was previously reported to adapt to stressors affecting cell envelope integrity through the Cpx-mediated repression of the nuo and cyo operons. The nuo and cyo transcriptional downregulation suggested that those proteins are toxic in the presence of envelope stress [66]. The membrane-bound systems for proton and electron transport Nuo (the NADH dehydrogenases I and II) and Cyo (cytochrome o oxidase) were downregulated at high pH [67]. Our analysis detected the downregulation of these systems, or part thereof, under all five investigated stressors.
Other genes that were downregulated under all studied stressors are: • fadL encodes the component of a channel involved in the import of long-chain (C12-C18) fatty acids (LCFA) across the bacterial outer membrane [68]. It is part of the fatty acid-degrading (fad), regulon which includes eight more genes involved in fatty acids catabolism [69]. Our analysis indicated that fadD, fadE, fadI, and fadR were downregulated under at least two stressors, suggesting that, during stress, E. coli downregulates the uptake of exogenous fatty acids along with its metabolism and degradation which generates various reduced cofactors. As previously reported, LCFA degradation generates oxidative stress with high levels of reactive oxygen species [70,71]. • cvpA encodes colicin V production protein [72] and is located directly upstream of purF and repressed by PurR, the main repressor of purine synthetic pathway [73]. Its deletion mutant in EHEC was highly sensitive to bile and was deemed important in cell envelope homeostasis in response to stressors, such as deoxycholate bile salt [74]. • ybhC encodes an outer membrane lipoprotein [75]. Genome-wide screening of an ASKA library found that overexpression of yhbC leads to an increase in the minimum inhibitory concentration (MIC) of hydrogen peroxide by more than three folds [76]. • gtrB (yfdH) and gtrS (yfdI) are part of the yfdGHI operon. Our analysis indicated that yfdG was downregulated under two stressors. yfdG, yfdH, and yfdI are homologous to the type IV O antigen modification genes (gtrAIV, gtrBIV, and gtrIV) in the genome of Shigella flexneri NCTC 8296 [77]. • rnb encodes the ribonuclease II enzyme (RNase II), involved in the specific degradation of mRNA in the 3 to 5 direction [78].
Network analysis identified upregulated hub genes, which included seven genes involved in biofilm formation (dosC, yhcN, yodD, dgcZ, mqsR, elfA, and glgS), suggesting that stress response can lead to biofilm formation. Hub genes also include genes involved in ferri-enterobactin transport (fepB, fepC, fepD, and fepG) and utilization (fes). entF was upregulated under three stressors (antibiotics, heat, and oxidative stress). It is part of the entABCDEF operon. Our analysis indicated that entA, entB, entC, entD, and entE were upregulated in antibiotic and oxidative stress. This operon is involved in enterobactin biosynthesis, which is a well-known siderophore. This highlights the importance of iron availability for E. coli to cope with different stressors. Downregulated hub genes include seven genes of flagellar biogenesis and assembly (flgH, flgG, fliK, fliM, flip, flhA, and flhB). A strong relation between motility suppression and high pH was previously reported [67]. Here, we believe that most of those genes are downregulated to preserve energy when E. coli is exposed to cold, heat, oxidative stress, or antibiotic treatment.
In general, the purpose of stress response to concentrate cellular resources on survival, as opposed to growth/multiplication. So, based on our analysis and several prior studies, we conclude that major stress responses involve maintaining or upregulating genes related to survival, biofilm formation, and stress response, while downregulating non-essential energy-requiring processes related to growth and replication. Overall, the globally upregulated pathways we identified for stress responses fall into three major categories: (i) cellular stress response, (ii) DNA repair, and (iii) cell adhesion and biofilm. On the other hand, the main downregulated pathways we identified for stress responses are: (i) de novo purine and pyrimidine biosynthesis pathways, salvage, and uptake, (ii) oxidative phosphorylation/aerobic respiration, and (iii) motility.
Of note, while this manuscript was in final preparation, a publication using a similar approach reported a set of genes and networks involved in a collection of stressors [8]. The studied stressors were nutrition limitations (rich M9 and poor M9), acidic stress (pH 5 ), antibiotics (trimethoprim (TMP) and chloramphenicol (CAM)), and oxidative stress conferred by growth in low oxygen (LOX) environments, yet the main emphasis of the paper was on transcription factors and regulatory networks. Although different computational approaches were used, the work agrees with ours in the following: nutrition limitations in E. coli indicated by the poor M9 medium were enriched in downregulation of the flagellar genes and upregulation of CsgD, a biofilm master regulator. TMP induced SOS response and cellular response to various stimuli. The SOS response inhibits cell division upon exposure to TMP by upregulating sulA. sulA was also upregulated in all five stressors in our analysis. TMP also induced upregulation of genes csgD, csgE, and csgF of the csgDEFG operon as well as csgB of the csgBA operon. In our study csgF and csgB were upregulated in three stressors, csgG was upregulated in two stressors, and heat stress upregulated csgA, csgD, and csgE. Those two operons involved in biofilm formation confer protection to bacteria against antibiotics and stressful conditions. Genes in the main five sub-networks generated in that study belonged to flagellar assembly, energy metabolism, SOS response and DNA repair, RNA binding proteins, and biosynthesis of amino acids and secondary metabolites. This is consistent with our finding as biofilm formation, DNA repair, and the bacterial response to stress and to various stimuli were enriched in upregulated genes. Functions such as energy-requiring metabolism, purine and pyrimidine biosynthesis, flagellar assembly, biogenesis, and export were enriched in downregulated genes. Individual genes in those main sub-networks were represented in our hub up or downregulated genes: nineteen genes of energy metabolism (elaB, dps, fbaB, fic, msyB, osmY, sra, wrbA, ybaY, ybgS, ycaC, yccJ, ycgB, yeaG, yeaH, yegP, yiaG, yodD, and yqjC; Figure 6A), four genes of DNA repair (lexA, recF, recN, and recQ; Figure 6D) and four genes of amino acids and secondary metabolites biosynthesis (entF, fepC, fepG, and fes; Figure 6B) were in hub upregulated genes; six flagellar genes (flgH, flgG, fliK, fliM, flip, and flhB; Figure 7B) and one gene of amino acids and secondary metabolites biosynthesis (cyoB; Figure 7C) were in hub downregulated genes.
Intriguingly, PurR was identified as a major regulator that is related to several of the genes and pathways detected in our analysis. PurR is considered one of the general transcription factors (TFs) as it controls the expression of genes related to connecting reactions as well as genes with diverse functional categories [79]. PurR function and size are considered limited when compared to other general TFs in E. coli, such as Crp [80], Fnr [81], and Lrp [82]. In Bacillus subtilis, the nucleotide messengers ppGpp and pppGpp (collectively (p)ppGpp) were reported to bind to PurR, therefore increasing its DNA binding and enhancing the downregulation of PurR-regulated genes involved in purine biosynthesis during nutrient starvation. This suggests that purine biosynthesis is repressed by PurR to cope with unfavorable conditions during starvation [83]. PurR has not been previously reported to be related to stress in E. coli, yet in our analysis, most of the downregulated genes are known to be repressed by PurR, and the purine and pyrimidine biosynthesis, salvage, and uptake pathways were highly represented among downregulated data. Taken together, these findings suggest that PurR has a critical role in balancing the metabolism E. coli under different types of stress.

Conclusions
In summary, this study aimed at using a systems approach, through integrating several data sets, to determine genes that are involved in the response to multiple stresses. Because multiple array platforms were used, we adopted a percentile rank-based approach to pick significantly altered genes, and we gave emphasis to gene sets or pathways rather than individual genes, as gene-level analysis is sensitive to noise, batch effects, and experimental variations.
The strengths of our approach are that it is unbiased and statistically robust and that it relies on network enrichment to overcome the false positive inclusion or false negative exclusion of individual genes. Most of the statistically significant results were also biologically relevant, or at least biologically explainable. For example, the downregulation of non-essential functions under stress and the upregulation of genes necessary for survival are consistent with an "emergency response" in which the cell manages its resources carefully. Meanwhile, a few genes remain unassigned with any functions and will be subject to further experimental validation. Like all computational approaches, the strategy followed here remains hypothesis generating rather than confirmatory and certainly requires future validation via experiments relying, for example, on gene deletion and insertion (loss and gain of function, respectively). In addition, like many high-throughput analyses, the study cannot establish causation, and some of the observed variations might be consequential while some may be collateral. In addition, the approach may overlook genes that may be biologically relevant but fall out of the statistical thresholds for one reason or another. This is exemplified by the absence of some operon members even when most of the operon is selected among stress response genes.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/microorganisms10091793/s1, Figure S1: Network of (A) 83 upregulated genes in at least four stressors with medium confidence (cut-off score: 0.4); Figure S2: The main sub-clusters of the upregulated genes in at least four stressors network; Tables S1-S5: Top 250 up-or downregulated genes in each study/sample in heat stress, oxidative stress, cold stress, nitrosative stress, and antibiotic stress, respectively; Table S6: Combined list of upregulated in all samples of each stress; Table S7