Combined Transcriptomic and Proteomic Profiling of E. coli under Microaerobic versus Aerobic Conditions: The Multifaceted Roles of Noncoding Small RNAs and Oxygen-Dependent Sensing in Global Gene Expression Control

Adaptive mechanisms that facilitate intestinal colonization by the human microbiota, including Escherichia coli, may be better understood by analyzing the physiology and gene expression of bacteria in low-oxygen environments. We used high-throughput transcriptomics and proteomics to compare the expression profiles of E. coli grown under aerobic versus microaerobic conditions. Clustering of high-abundance transcripts under microaerobiosis highlighted genes controlling acid-stress adaptation (gadAXW, gadAB, hdeAB-yhiD and hdeD operons), cell adhesion/biofilm formation (pgaABCD and csgDEFG operons), electron transport (cydAB), oligopeptide transport (oppABCDF), and anaerobic respiration/fermentation (hyaABCDEF and hycABCDEFGHI operons). In contrast, downregulated genes were involved in iron transport (fhuABCD, feoABC and fepA-entD operons), iron-sulfur cluster assembly (iscRSUA and sufABCDSE operons), aerobic respiration (sdhDAB and sucABCDSE operons), and de novo nucleotide synthesis (nrdHIEF). Additionally, quantitative proteomics showed that the products (proteins) of these high- or low-abundance transcripts were expressed consistently. Our findings highlight interrelationships among energy production, carbon metabolism, and iron homeostasis. Moreover, we have identified and validated a subset of differentially expressed noncoding small RNAs (i.e., CsrC, RyhB, RprA and GcvB), and we discuss their regulatory functions during microaerobic growth. Collectively, we reveal key changes in gene expression at the transcriptional and post-transcriptional levels that sustain E. coli growth when oxygen levels are low.


Introduction
Escherichia coli is a Gram-negative commensal bacterium that commonly inhabits the intestines of humans and other animals under microaerobic or anaerobic conditions. Previous studies have shown that E. coli growth at different concentrations of oxygen involves substantial reprogramming of the gene expression controlled by several transcription factors, E. coli MG1655 strain was grown on minimal medium containing glucose as a carbon source under continuous aerobic or microaerobic conditions. These experiments were carried out in a benchtop fermentor (Winpact Parallel Fermentation System FS-05-220) to obtain multiple biological replicates (namely, biological replicates of 5 aerobic and 10 microaerobic cultures) for further analysis. Cell doubling times of aerobic and microaerobic cultures were 77.9 ± 8.6 and 245.5 ± 24.7 min, respectively. To prepare total RNA [17] or protein samples, cells were grown to mid-logarithmic phase, corresponding to OD 460 0.5-0.6. Samples of purified RNA were used for the RNA deep-sequencing analyses, and the same batches of RNA were used for Northern blot validation. Details of the bioreactor culture conditions, RNA deep-sequencing procedures, and reagents used are provided in the Materials and Methods.
More than 9 million raw sequencing reads were generated for each RNA sample. After trimming the raw sequencing reads, the high-quality unique sequences were mapped to the E. coli K12 substrain MG1655 reference genome (NC_000913) [18]. On average, mapping coverages of 94.6% and 88.8% were obtained for total unique sequence reads under aerobic and microaerobic growth conditions, respectively (Table S1). With a threshold of ≥1 transcript per million mapped reads (TPM), we detected expression of 4388 and 4385 genes under aerobic and microaerobic growth conditions, respectively (NCBI GEO accession # GSE189154).
To assess the robustness of our datasets, we calculated correlations across biological replicates. In Figure S1A, we show values for the correlation coefficients for all pairwise scatterplots obtained for aerobic or microaerobic samples. We identified strong correlations between biological replicates representing the same growth condition. More specifically, correlation coefficients for aerobic and microaerobic transcriptomes were greater than 0.97 and 0.91, respectively, implying that our results are highly reproducible. In addition, we performed principal component analysis (PCA) on a combined dataset that successively maximized variance across all datasets (i.e., datasets O-1~O-5 and N-1~N-10). Our PCA revealed two distinct groups (i.e., O-1~O-5 and N-1~N-10) corresponding to aerobic and microaerobic cultures, respectively ( Figure S1B). Together, our correlation and PCA analyses justify the use of all our RNA-seq datasets for further gene expression analyses.

Defining Major Gene Clusters Involved in Adaptation to Microaerobiosis
To identify differentially expressed genes (DEGs), we ranked the processed data on detected gene transcripts under aerobic (4388 genes) and microaerobic (4385 genes) conditions according to false discovery rate (FDR)-adjusted p-values and log 2 fold-change (FC).

Identification of DEGs under Microaerobic Versus Aerobic Conditions
To compare the gene expression patterns of microaerobic and aerobic cultures, we calculated the normalized transcript expression values (expressed as TPM) for each gene with a p-value ≤ 0.05 and selected 280 DEG transcripts displaying log 2 fold-change ≥ 2, representing 176 upregulated and 104 downregulated genes, respectively ( Figure 1A). Next, we filtered these DEGs using a FDR threshold of −log 10 (p-value) ≥ 1 ( Figure 1B), which resulted in 176 DEGs (105 upregulated and 71 downregulated genes) that showed significant differential expression. aerobic and microaerobic cultures, respectively ( Figure S1B). Together, our correlatio and PCA analyses justify the use of all our RNA-seq datasets for further gene expressio analyses.

Defining Major Gene Clusters Involved in Adaptation to Microaerobiosis
To identify differentially expressed genes (DEGs), we ranked the processed data o detected gene transcripts under aerobic (4388 genes) and microaerobic (4385 genes) con ditions according to false discovery rate (FDR)-adjusted p-values and log2 fold-chang (FC).

Identification of DEGs under Microaerobic versus Aerobic Conditions
To compare the gene expression patterns of microaerobic and aerobic cultures, w calculated the normalized transcript expression values (expressed as TPM) for each gen with a p-value ≤ 0.05 and selected 280 DEG transcripts displaying log2 fold-change ≥ representing 176 upregulated and 104 downregulated genes, respectively ( Figure 1A Next, we filtered these DEGs using a FDR threshold of -log10 (p-value) ≥ 1 ( Figure 1B which resulted in 176 DEGs (105 upregulated and 71 downregulated genes) that showe significant differential expression.

Functional Clusters of DEGs
We present a comprehensive overview of the identified DEGs in Table 1, reveal key gene clusters, genes/operons, biological functions, small regulatory RNAs, and lated transcription factors involved in adaptation to changing oxygen conditions. Upr ulated ( Figure 2A) and downregulated DEGs ( Figure 2B) were functionally classified i three general gene ontology (GO) categories, i.e., biological processes, cellular comp nents, and molecular functions. Based on significant fold enrichment (≥10 compared w reference genes in the same subcategory), the upregulated DEGs we identified w mainly found in the following GO subcategories: peptidoglycan-associated pept transport, oxidative phosphorylation, Ni-Fe hydrogenase complex, peptidoglycan p tide transmembrane transporter activity, peptidoglycan transmembrane transporter tivity, hydrogenase (acceptor) activity, and oxidoreductase activity acting on hydrogen donor (Figure 2A). Through the same process, we identified 76 GO subcategories downregulated DEGs, including TCA cycle, several transport processes (such as fer hydroxamate import into cell, iron import into cell, copper ion export), ion homeosta chemical homeostasis, enterobactin biosynthetic process, nonribosomal peptide biosy thetic process, secondary metabolite biosynthetic process, lactone metabolic process, tibiotic metabolic process, stress response to metal ion, detoxification of inorganic co . DEGs displaying statistical significance (i.e., meeting this FC criterion) are shown as red (176 upregulated genes) or green (104 downregulated genes) dots. (B) Volcano plot displaying FC plotted against the false discovery rate (FDR) p-value. The y axis represents the −log 10 FDR p-value and the x axis represents the log 2 FC value. The horizontal black line indicates the significance threshold (−log 10 p-value ≥ 1), and the vertical black lines indicate the FC threshold (absolute value of log 2 FC ≥ 2). DEGs displaying statistical significance (i.e., those meeting both criteria) are shown as 105 upregulated (red dots) and 71 downregulated (green dots) genes in the right-upper and left-upper areas of the panel delineated by black lines, respectively.

Functional Clusters of DEGs
We present a comprehensive overview of the identified DEGs in Table 1, revealing key gene clusters, genes/operons, biological functions, small regulatory RNAs, and related transcription factors involved in adaptation to changing oxygen conditions. Upregulated ( Figure 2A) and downregulated DEGs ( Figure 2B) were functionally classified into three general gene ontology (GO) categories, i.e., biological processes, cellular components, and molecular functions. Based on significant fold enrichment (≥10 compared with reference genes in the same subcategory), the upregulated DEGs we identified were mainly found in the following GO subcategories: peptidoglycan-associated peptide transport, oxidative phosphorylation, Ni-Fe hydrogenase complex, peptidoglycan peptide transmembrane transporter activity, peptidoglycan transmembrane transporter activity, hydrogenase (acceptor) activity, and oxidoreductase activity acting on hydrogen as donor (Figure 2A). Through the same process, we identified 76 GO subcategories for downregulated DEGs, including TCA cycle, several transport processes (such as ferric hydroxamate import into cell, iron import into cell, copper ion export), ion homeostasis, chemical homeostasis, enterobactin biosynthetic process, nonribosomal peptide biosynthetic process, secondary metabolite biosynthetic process, lactone metabolic process, antibiotic metabolic process, stress response to metal ion, detoxification of inorganic compound, energy transducer activity, signaling receptor activity, and Fe 2 S 2 cluster binding, among others ( Figure 2B). To better define the molecular interactions, reactions, and relationship network of the biogenesis pathways that are differentially affected under microaerobic versus aerobic conditions, we conducted a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (p-value ≤ 0.05) on the DEGs. This analysis revealed four and seven KEGG pathways that were upregulated or downregulated, respectively ( Figure 2C,D, left panels, respectively). The upregulated KEGG pathways were β-alanine metabolism, nitrotoluene degradation, quorum sensing, and β-lactam resistance, whereas the downregulated KEGG pathways were propanoate metabolism, carbon metabolism, biosynthesis of secondary metabolites, ABC transporters, biosynthesis of antibiotics, TCA cycle, and biosynthesis of siderophore group nonribosomal peptides.
We also employed the UniProt Knowledgebase (UniProtKB) keyword database (EMBL-EBI, Cambridge, UK; SIB Swiss Institute of Bioinformatics, Geneva, Switzerland; PIR, Washington, DC, USA) to further characterize the retrieved pathways pertaining to DEGs. By assigning DEGs to functional, structural, or other UniProtKB keyword categories, we essentially generated a highly similar classification to those determined from our GO and KEGG pathway analyses. More specifically, UniProtKB keywords associated with upregulated DEGs were iron, electron transport, peptide transport, and membrane, whereas copper, phosphopantetheine, copper transport, transmembrane beta strand, ligase, Fe 2 S 2 , receptor, bacteriocin transport, TonB box, TCA cycle, transport, enterobactin biosynthesis, iron, ion transport, and iron transport were all associated with downregulated DEGs ( Figure 2C,D, right panels, respectively).

Prophage-and Phage-Related Genes
Prophage-related genes constitute up to 13.5% of the E. coli genome [34,35], contributing to bacterial survival in hosts by increasing cell fitness and virulence. Recent studies have revealed that the expression of such prophage-related genes can (i) increase E. coli resistance to adverse conditions such as exposure to antibiotic, acid, oxidative, or osmotic stress; and (ii) influence metabolic remodeling, biofilm formation, cell movement, and growth [36][37][38][39][40]. Although our recent works revealed that one such prophage gene, dicF [41], plays a critical role in regulating cell division under anaerobic conditions [42], how other prophage-related genes are expressed and function under oxygen-limited conditions remains unclear.
To identify other prophage-related genes that potentially contribute to cell fitness and survival in microaerobic environments, we compared the expression of prophage-and phage-related genes under microaerobic versus aerobic conditions. Using our RNA-seq transcriptomic dataset and based on 245 known prophage-related genes within 134 operons [2,33], we detected the expression of 200 prophage-related genes within 118 operons, reflecting 132 upregulated and 68 downregulated genes. As summarized in Table 2, highly upregulated (log 2 FC ≥ 1.5) prophage-and phage-related genes could be assigned to 15 operons, whereas the downregulated ones were solely localized in the fhuACDB operon that codes for proteins involved in ferrichrome transport. Notably, FhuA protein can also serve as a phage receptor [43].

Prophage-and Phage-Related Genes
Prophage-related genes constitute up to 13.5% of the E. coli genome [34,35], contributing to bacterial survival in hosts by increasing cell fitness and virulence. Recent studies have revealed that the expression of such prophage-related genes can (i) increase E. coli resistance to adverse conditions such as exposure to antibiotic, acid, oxidative, or osmotic stress; and (ii) influence metabolic remodeling, biofilm formation, cell movement, and growth [36][37][38][39][40]. Although our recent works revealed that one such prophage gene, dicF showing the fold change log 2 value. Protein (increased/decreased) showing the abundance ratio value. Undetected transcripts or proteins are in black or indicated by "X", respectively. Non-prophage related gene marking within  Interestingly, some of these transcripts encode TFs that control metal ion homeostasis; for instance, CusR for copper, and LscR, Dps, and FecI for iron. Theoretically, when the level of a TF increases, its target genes should be controlled accordingly, depending on whether the TF is a positive or negative regulator. We selected the IscR regulon for further analysis. Levels of the iscR mRNA decreased~4.4-fold (log 2 = 2.13; Table 3) under microaerobic conditions relative to those under aerobic growth, and the level of the respective protein, lscR, consequently also diminished (~1.3-fold; Table S2). This reduction in IscR abundance, along with an anticipated reduction in [Fe 2 S 2 ] iron-sulfur cluster availability (Table 1) required for IscR activity, could be accountable, at least in part, for the nearly 3-fold upregulation of torT (log 2 = 1.67; Table 4). The changes in IscR abundance had the most pronounced effect on the nrdHIEF operon, the genes of which exhibited a 10.8-59.3 foldchange (log 2 = 3.43 to log 2 = 5.89 in Table 4) in downregulated DEG expression.
52 Log 2 fold change values for up-and down-regulated mRNA and abundance ratio value of increased and decreased protein (in red and green, respectively) are indicated. Undetected molecules are marked by "X". *: No detected value under aerobic condition leading to a huge protein abundance ratio that is marked as "100".
Log 2 fold change values for up-and down-regulated mRNA and abundance ratio value of increased and decreased protein (in red and green, respectively) are indicated.
The RyhB sRNA acts as a global regulator of iron homeostasis [10]. We observed extremely low abundance of this sRNA under microaerobic conditions (~8-fold less compared to levels under aerobic conditions; log 2 = 3, Table 5) and it was not detectable by Northern blotting (see Result Section 2.3), an outcome consistent with Fur-dependent repression of ryhB transcription.

Identification of Differentially Expressed sRNAs and Northern-Blot-Based Validation
sRNAs are common in bacteria, where they play critical roles in regulating a wide range of cellular functions [45]. Our RNA-seq dataset also revealed differential expression of a number of sRNAs under microaerobic and anaerobic conditions. Of the 64 known sRNAs in E. coli [2,33], 18 exhibited >1.5-fold difference in abundance under microaerobic versus aerobic conditions ( Figure 3A). We employed Northern blot analysis to validate these results. Consistently, we observed that the abundance of RyhB was dramatically reduced under oxygen-limited conditions, rendering it almost undetectable under microaerobiosis ( Figure 4A, second panel from right). In contrast, the levels of several other sRNAs (e.g., CsrC, GcvB and RprA) were considerably higher under microaerobiosis ( Figure 4A). We detected two or more species of some sRNAs (e.g., GadY and RprA ( Figure 4A), GlmY, and RyeA ( Figure 4B)). To test whether the increase/decrease in abundance of certain sRNAs could be attributable to their higher metabolic stability, we also determined the half-lives of three sRNAs-namely CsrB, CsrC, and RyhB-by inhibiting their transcription by means of rifampicin treatment, and then determined their time-dependent decreased abundance using Northern blot analysis. As shown in Figure 4A-C,E-G, the half-lives of both CsrC and CsrB increased moderately under microaerobiosis (from 4.3 to 6.9 min and from 3.8 to 5.9 min, respectively), though only the steady-state level of CsrC increased dramatically ( Figure 5H). This outcome indicates that it is more efficient transcription, rather than increased RNA stability, that is responsible for the higher CsrC abundance under microaerobic conditions. Given that RyhB was almost undetectable in cells cultured under oxygen-limited conditions, we could not directly compare the half-life of this sRNA under different conditions ( Figure 5I-L). In addition to validating sRNA levels by Northern blotting, we also assessed levels of Hfq [46]-an RNA chaperone that plays an important role in facilitating sRNA/mRNA interactions-by means of Western blotting, which revealed only minor differences in abundance under aerobic and microaerobic conditions ( Figure 5M,N).

Identification of Differentially Expressed sRNAs and Northern-Blot-Based Validation
sRNAs are common in bacteria, where they play critical roles in regulating a wide range of cellular functions [45]. Our RNA-seq dataset also revealed differential expression of a number of sRNAs under microaerobic and anaerobic conditions. Of the 64 known sRNAs in E. coli [2,33], 18 exhibited >1.5-fold difference in abundance under microaerobic versus aerobic conditions ( Figure 3A). We employed Northern blot analysis to validate these results. Consistently, we observed that the abundance of RyhB was dramatically reduced under oxygen-limited conditions, rendering it almost undetectable under microaerobiosis ( Figure 4A, second panel from right). In contrast, the levels of several other sRNAs (e.g., CsrC, GcvB and RprA) were considerably higher under microaerobiosis (Figure 4A). We detected two or more species of some sRNAs (e.g., GadY and RprA ( Figure  4A), GlmY, and RyeA ( Figure 4B)). To test whether the increase/decrease in abundance of certain sRNAs could be attributable to their higher metabolic stability, we also determined the half-lives of three sRNAs-namely CsrB, CsrC, and RyhB-by inhibiting their transcription by means of rifampicin treatment, and then determined their time-dependent decreased abundance using Northern blot analysis. As shown in Figure 4A-C and 4E-G, the half-lives of both CsrC and CsrB increased moderately under microaerobiosis (from 4.3 to 6.9 min and from 3.8 to 5.9 min, respectively), though only the steady-state level of CsrC increased dramatically ( Figure 5H). This outcome indicates that it is more efficient transcription, rather than increased RNA stability, that is responsible for the higher CsrC abundance under microaerobic conditions. Given that RyhB was almost undetectable in cells cultured under oxygen-limited conditions, we could not directly compare the halflife of this sRNA under different conditions ( Figure 5I-L). In addition to validating sRNA levels by Northern blotting, we also assessed levels of Hfq [46]-an RNA chaperone that plays an important role in facilitating sRNA/mRNA interactions-by means of Western blotting, which revealed only minor differences in abundance under aerobic and microaerobic conditions ( Figure 5M,N).    (left and middle panels, respectively) and csrA (right panel) mRNAs within the E. coli genome. The y axis represents the number of RNA-seq reads for the sRNAs and csrA mRNA on the largest scale of 400,000 and 4000, respectively. The coding region of each gene is shown in blue at the top of each panel, and expression is shown in blue and red for aerobic and microaerobic growth conditions, respectively.  conditions. Mean values for CsrB (C), CsrC (G), and RyhB (K) half-lives under aerobic and microaerobic conditions are shown (encompassing three biological repeats, bars represent standard error). The dotted gray line indicates 50% of total RNA remaining. Black circles and blue squares represent the signal intensities corresponding to RNA samples from aerobic and microaerobic cultures, respectively. CsrB, CsrC, and RyhB half-lives under aerobic conditions were calculated as 3.8 ± 0.2, 4.3 ± 0.4, and 7.4 ± 0.1 min (C,G,K), respectively, whereas under microaerobic conditions they were 5.5 ± 0.3 min, 6.1 ± 0.4 min, and no detectable signal (see panels (B,F,J)), respectively. Bar graph shows the relative steady-state levels of small RNAs (time 0) normalized to their levels under aerobic conditions, which were arbitrarily set as 1. Experiments were performed with three biological replicates and representative images are shown. The steady-state level of CsrB under microaerobiosis relative to aerobiosis was 1.03 ± 0.05-fold (p-value = 0.49) (D), whereas for CsrC it was 5.32 ± 0.51-fold (p-value < 0.0001, indicated as ****) (H). Expression of RyhB was not detectable (nd) under microaerobic conditions (L). (M) Hfq protein abundance analyzed via Western blotting. Equal amounts of total protein were fractionated in 20% SDS polyacrylamide gels and transferred to a membrane, and the lower part of the membrane was probed with anti-Hfq antibody. The upper part of the membrane was used to detect GAPDH as a loading control. Experiments were performed with three biological replicates and representative images are shown. (N) Quantification of Hfq level. The signal obtained with anti-Hfq antibody was normalized using GAPDH and further processed to calculate the relative protein expression level, plotted as vertical bars. Hfq level under microaerobiosis was normalized to its level under aerobiosis, which was arbitrarily set as 1. The difference in Hfq level under these conditions was not statistically significant (p-value = 0.26).

Proteome Analysis Corroborates Differential Protein Abundance under Changing Oxygen Availability
To further validate our DEG analyses, we adopted a quantitative proteomic approach to analyze differential protein abundance under aerobic and microaerobic growth conditions. We conducted this analysis on aliquots of the same batches of cultured cells used for our above-described RNA-seq assays, encompassing two biological repeats for both growth conditions (samples O-1 and O-2 for aerobiosis; samples N-2 and N-3 for microaerobiosis).

Identification of Differentially Abundant Proteins under Microaerobic Versus Aerobic Conditions
We deployed commercially available isobaric iTRAQ * mass tags [47,48] to simultaneously analyze multiple biological samples. The identical masses and chemical properties of these isobaric tags enabled co-elution of various isotopologues. The isobaric tags of peptides were cleaved by collision-induced dissociation (CID) during tandem mass spectrometry (MS/MS), before assessment of peptide fragment ions to define their sequence and quantitation of the isobaric tags. Peptide identification and relative quantitation were determined concurrently. In total, we identified 1498 and 1488 proteins from cells grown under aerobic or microaerobic conditions, respectively (Table S2). Pairwise scatterplots revealed strong correlations between two biological repeats for the same growth condition (r = 0.98 and 0.97 for aerobic and microaerobic growth, respectively) ( Figure S2).
We used protein abundance ratios to identify proteins that were differentially abundant under microaerobic growth conditions relative to aerobic growth. Using 95% confidence intervals with two standard deviations (SD), we assumed that ratios >1.39 or <0.698 (i.e., ratios with more than 1.39-fold increase or more than 1.43-fold decrease in protein abundance) indicated significant changes. Accordingly, we set the protein abundance ratio threshold to 1.5 and found 113 and 92 proteins that had increased and decreased in abundance, respectively, in the E. coli MG1655 cells grown under microaerobic versus aerobic conditions (Table S3). We explore the consistency among our proteomic and RNAseq datasets in the Discussion.

Biological Function Operon
Translational Regulator *

Activator Inhibitor Attenuator
Gene name rpsMKD-rpoA-rplQ   Operon gene name showed red letter: increased; bold red letter with protein ratio ≥ 1.5; black letter: not detected. Operon gene name showed green letter: decreased; bold green letter with protein ratio ≤ −1.5; black letter: not detected. *: According to the RegulonDB & EcoCyc databases.  The protein-protein interaction networks for increased (A) and decreased (B) differentially abundant proteins were generated using the STRING platform (https://string-db.org/). Abundance-increased proteins were involved in processes such as glycolysis, ATP metabolism, and coenzyme/small-molecule metabolism, for which proteins are represented in red, blue, and green, respectively. Abundancedecreased proteins were involved in ribosome biogenesis, post-transcriptional regulation of gene expression and peptide metabolism, and are indicated in red, blue, and green, respectively.

Discussion
The ability of enterobacteria to colonize the digestive systems of mammals is dependent on their capacity to adapt and thrive in low-oxygen environments. Previous studies [1-16] of E. coli revealed that the transition from aerobic to microaerobic/anaerobic conditions requires a substantial reprogramming of gene expression, which greatly affects the bacterial lifestyle and major cellular functions such as metabolism, transport, and energy production. Nevertheless, many details of the underlying regulatory networks, including post-transcriptional mechanisms that coordinate E. coli adaptation and survival under oxygen-limited conditions, remain poorly defined and merit further analysis.
In this study, we employed a combination of transcriptomic and proteomic approaches to analyze differences in RNA and protein abundance for E. coli grown in minimal medium under microaerobic versus aerobic conditions. Significantly, the E. coli cells were grown in a Bench-Top Fermentor, allowing us to control key parameters (i.e., temperature, composition of the medium, oxygen level, and pH) and thereby ensuring that growth conditions were equivalent for each experimental culture.
Our transcriptomic analysis of aerobic and microaerobic cultures uncovered numerous upregulated and downregulated genes. Annotation and functional clustering of DEGs using the web-based tools available on the STRING, GO, KEGG, UniProtKB, RegulonDB,  . (A,B) The proteinprotein interaction networks for increased (A) and decreased (B) differentially abundant proteins were generated using the STRING platform (https://string-db.org/). Abundance-increased proteins were involved in processes such as glycolysis, ATP metabolism, and coenzyme/small-molecule metabolism, for which proteins are represented in red, blue, and green, respectively. Abundancedecreased proteins were involved in ribosome biogenesis, post-transcriptional regulation of gene expression and peptide metabolism, and are indicated in red, blue, and green, respectively.

Discussion
The ability of enterobacteria to colonize the digestive systems of mammals is dependent on their capacity to adapt and thrive in low-oxygen environments. Previous studies [1-16] of E. coli revealed that the transition from aerobic to microaerobic/anaerobic conditions requires a substantial reprogramming of gene expression, which greatly affects the bacterial lifestyle and major cellular functions such as metabolism, transport, and energy production. Nevertheless, many details of the underlying regulatory networks, including post-transcriptional mechanisms that coordinate E. coli adaptation and survival under oxygen-limited conditions, remain poorly defined and merit further analysis.
In this study, we employed a combination of transcriptomic and proteomic approaches to analyze differences in RNA and protein abundance for E. coli grown in minimal medium under microaerobic versus aerobic conditions. Significantly, the E. coli cells were grown in a Bench-Top Fermentor, allowing us to control key parameters (i.e., temperature, composition of the medium, oxygen level, and pH) and thereby ensuring that growth conditions were equivalent for each experimental culture.
Our transcriptomic analysis of aerobic and microaerobic cultures uncovered numerous upregulated and downregulated genes. Annotation and functional clustering of DEGs us-ing the web-based tools available on the STRING, GO, KEGG, UniProtKB, RegulonDB, and EcoCyc internet platforms revealed that oxygen level directly influences carbon metabolism, energy production, metal ion homeostasis, and cell envelope functions.
Previous studies have shown that the generation of proton motive force by cytochrome bo oxidase facilitates ATP production by E. coli ATP synthase under aerobic conditions [49]. However, this process can become less efficient at low oxygen concentrations, potentially requiring the action of a second cytochrome oxidase (i.e., cytochrome bd-1 oxidase) with a much higher affinity for molecular oxygen. Indeed, our transcriptomic and proteomic data clearly show that the cydABX operon, which encodes this secondary oxidase, was strongly upregulated in E. coli grown under microaerobic conditions (see Table 1). Moreover, the same operon hosts another upregulated gene, namely ndh, which encodes NADH ubiquinone oxidoreductase (a co-factor of cytochrome bd-1 oxidase), responsible for 2a ubiquinol regeneration.
ATP production through oxidative phosphorylation is mainly efficient during aerobic growth, with the TCA cycle greatly contributing to this process by producing the NADH and FADH that feed into the respiratory cycle. This latter occurs when succinate:quinone oxidoreductase (encoded by sdhCDA) converts succinate to fumarate and, concurrently, reduces ubiquinone to ubiquinol. Given the diminished role of oxidative phosphorylation under microaerobic conditions, the expression of genes coding for TCA cycle enzymes is likely reduced due to their repression by the FNR and ArcA TFs. Consistently, we observed downregulation of several operons coding for enzymes involved in steps 4 (sucAB; 2-oxoglutarate decarboxylase), 5 (sucCD; succinyl-CoA synthetase), 6 (sdhCDAB; succinate dehydrogenase), and 7 (fumAC; fumarase) of the TCA cycle. Furthermore, reduced production of these enzymes that host numerous Fe-S clusters implies a reduced cellular need for iron. Indeed, we found that many genes involved in iron homeostasis (fhuACDB, fepA-entD, fes-ybdZ-entF-fepE, tonB, feoABC) and the production of protein complexes responsible for iron incorporation into Fe-S clusters (sufABCDSE and iscRSUA) were repressed during microaerobic growth. Their repression is likely mediated by the TF Fur, which is activated under microaerobic conditions in the presence of free Fe 2+ ions [10]. Simultaneously, Fur-dependent upregulation of ftnA (ca. 19-fold) and bfr (ca. 2.5-fold) elevates levels of the iron-storage proteins ferritin (2.4-fold) and bacterioferritin (1.6-fold), which efficiently sequester free iron atoms. Higher abundances of both these proteins are likely attributable to the lack of RyhB-mediated translational repression, since concentrations of this sRNA are extremely low during microaerobic growth (see below for further details).
Our transcriptomic analysis also highlighted an enhanced role for mixed-acid fermentation (MAF) during microaerobic growth, arising from ArcA-and FNR-mediated upregulation. E. coli cells employ MAF to convert glucose into various end-products such as formate, succinate, acetate, lactate, and ethanol. We detected upregulation of pflB under microaerobiosis, suggesting an increased production of formic acid and its subsequent conversion to hydrogen (H 2 ) and carbon dioxide (CO 2 ) by the formate hydrogenlyase complex encoded by the hycABCDEFGHI operon, with the latter also being upregulated under oxygen-limited conditions. Moreover, H 2 production appears to be coupled to reduction of menaquinone and the periplasmic protons responsible for the protein motive force that drives ATP production. This reaction is carried out by hydrogenase 1, which is encoded by another operon (i.e., hyaABCDEF) that is strongly upregulated in oxygen-poor environments. Similarly, we observed increased expression of several genes involved in the production of other known products of the MAF pathway, namely lactate (ldhA), ethanol and acetate (adhE), and succinate (i.e., fumB and frdB) ( Table 1).
Although we anticipated observing some changes in central carbon metabolism and energy production under microaerobiosis, the differential expression of some of the other major gene clusters is somewhat puzzling. For instance, the reasons for upregulation of multiple genes involved in acid (low pH) responses (i.e., gadAXW, gadBC, hdeAB-yhiD, hdeD) and oligopeptide transport (i.e., oppABCDF) are unclear. Since our cell cultures were continuously grown in a fermentor, the pH of the medium and its content was consistent throughout both aerobic and microaerobic growth, indicating that enhanced expression of these operons was not attributable to any other environmental factor except oxygen limitation. Thus, the low oxygen concentration in the environment may serve as a signal for E. coli to adapt to acidic environments and oligopeptide availability. Both these scenarios are encountered by enterobacteria upon entering the mammalian digestive system, which is characterized (at least in some regions) by low pH and the presence of oligopeptides produced from food digestion (i.e., polypeptide digestion by proteases). This observation raises an intriguing hypothesis that low oxygen concentrations might serve as a universal signal to alert bacterial cells that they have entered a host digestive system. Apart from microaerobiosis promoting expression of respiratory and acid stress response genes, we also detected clear upregulation of operons involved in biofilm formation (i.e., csgDEFG and pgaABCD) under this condition. Biofilm production could be considered an adaptive strategy allowing E. coli to survive in low-oxygen environments.
The expression patterns revealed by our transcriptomic analysis were largely confirmed by our proteomic data (Table 1). An unexpected exception was the regulation of the cusCFBA operon. Despite a decrease in the abundance of this polycistronic mRNA, the levels of each of the proteins encoded by this operon were increased (Table 1). The CusCFBA copper/silver efflux system contributes to maintaining copper homeostasis in low-oxygen environments [50]. The individual components of the tripartite CusCBA complex exist in a disassembled form to maintain the plasticity of the periplasm and its dynamic functions [51]. Although an increase in CusCFBA protein levels under microaerobic conditions is consistent with the documented role of this complex in copper tolerance at low oxygen concentrations [50], the exact transcriptional and post-transcriptional mechanisms responsible for the observed changes in cusCFBA expression at the RNA and protein levels are currently unknown.
In addition, our proteomic data revealed a considerable reduction in the abundance of many ribosomal (r-) proteins under microaerobic conditions (Table 6 and Table S3). E. coli r-proteins are encoded by polycistronic operons, and they are normally autoregulated at the translational level [52]. This regulatory mechanism involves the respective free r-proteins binding to their own polycistronic mRNAs and inhibiting their translation (e.g., autoregulation of the rpsJ-rplCDWB-rpsS-rplV-rpsC-rplP-rpmC-rpsQ operon by L4). As the structures of the L4 binding sites in the polycistronic mRNA and ribosomal RNA closely mimic each other, L4 can act as an efficient inhibitor of its own mRNA only when it is present in excess relative to ribosomes and, therefore, is available for interaction with its cognate mRNA. In other words, a decrease in the concentration of ribosomes under microaerobic conditions should release r-proteins to inhibit translation of their cognate mRNAs, thereby reducing their abundance in vivo. Moreover, an additional extraribosomal function of L4 is to change the abundance of numerous mRNAs, mainly by inhibiting RNase E-dependent mRNA decay during bacterial adaptation to adverse environments [53]. Another extraribosomal function of L4 is to post-transcriptionally regulate Tna expression in the stationary phase of growth through its direct binding to the tna intergenic region [54]. Imbalances in ribosomal synthesis can release ribosomal proteins to perform other extraribosomal functions [55]. Thus, under microaerobic growth conditions, the abundance of many ribosomal proteins is reduced, possibly leading some free ribosomal proteins to perform other extraribosomal functions to maintain cell fitness.
Interestingly, although the FNR TF potentially activates the transcription of multiple genes in oxygen-limited environments, our data suggest that many FNR-dependent genes remain silent, apparently due to repression by other factors acting at the transcriptional and post-transcriptional levels. For example, FNR-mediated gene activation of nitrate reductase does not occur in the absence of nitrate and may additionally be inhibited by sRNAs such as RprA [56]. Indeed, we found that levels of RprA were considerably higher during microaerobic growth, supporting its role in controlling nitrate respiration. In fact, apart from TF-mediated gene expression (such as through Fur, FNR, and ArcA, among others), sRNAs are also widely employed by E. coli to exert post-transcriptional control. sRNAs in E. coli range from~50 nucleotides (nt) (e.g., DicF; 53 nt) to >300 nt (e.g., CsrB; 369 nt) in length. However, detection of very short sRNAs by RNA-seq can be achieved only by including additional steps (i.e., specific size selection) in standard protocols, which are performed after fragmentation of the purified RNA and prior to cDNA library construction. Indeed, the shortest sRNA we detected was RdlD (66 nt), and we did not detect DicF (53 nt).
Our assessment of sRNA abundance uncovered several that were differentially expressed under microaerobic conditions, and their expression patterns were confirmed by Northern blotting. Particularly notable was the substantial decrease in RyhB concentration under those conditions. It is conceivable that the reduced abundance of this sRNA is inversely correlated with levels of its targets. Indeed, we observed higher abundances of the sodB, ftnA, and bfr mRNAs and their translational products (i.e., superoxide dismutase B and the two iron storage proteins FtnA and Bfr, respectively). However, low RyhB abundance did not similarly increase expression of other known RyhB targets located in the iscRSUA, sucCDAB, and sdhCDAB operons (Table 1), which are known for their roles in assimilating iron and homeostasis of that ion in many essential metabolic enzymes. In fact, we detected diminished abundance of their respective transcripts (see downregulated clusters in Table 1), which may be attributable to transcriptional repression by other global regulators such as the TFs FNR [57] and ArcA [58]. The latter regulators are known to downregulate the sucCDAB and sdhCDAB operons under microaerobic conditions. Moreover, the decreased expression of sdhCDAB we report was likely due to its repression by Fur, another TF that greatly impacts gene expression during anaerobic growth [10].
Unlike RyhB, we identified a number of sRNAs as being more abundant under microaerobic conditions (i.e., CsrB, CsrC, GcvB, and RprA) (Figure 3). CsrB and CsrC exert their regulatory functions by binding to the translational inhibitor CsrA, thereby preventing interaction of the latter with the translation initiation regions of numerous transcripts, including pgaABCD (see upregulated clusters in Table 1), under microaerobic conditions. Translational activation of this operon via competitive binding of CsrB and CsrC to CsrA enhances polysaccharide biosynthesis, thereby promoting biofilm formation. Similarly, CsrB-and CsrC-competitive binding mechanisms are likely involved in the translational activation of other genes (e.g., iraD [59] and glgS [60]) that are likewise upregulated under microaerobic conditions. The iraD gene encodes an anti-adapter protein that inhibits RssB-mediated degradation of the sigma stress factor RpoS, whereas GlgS is known as an inhibitor of cell motility [61]. CsrA often acts as a translational repressor, but it can also activate gene expression [60,61]. Although previous integrated transcriptomic data [62,63] have indicated that CsrA globally controls the levels of a large number of transcripts, the specific role of this translational repressor under microaerobic growth, i.e., when sRNAs CsrB/C are much more abundant than csrA mRNA (>60-fold; Figure 4C), remains to be determined. Interestingly, in comparison to the relatively short half-lives of CsrB (1.4 min) and CsrC (2.2. min) in Luria-Bertani (LB) medium [64], our substitution of rich medium (LB) with minimal medium and depletion of oxygen synergistically increased the stability of these small RNAs. Since both CsrB and CsrC are regulatory RNAs and substrates of RNase E, their increased stability could be attributable, at least in part, to lower RNase E levels under microaerobic conditions [42].
Notably, levels of another sRNA, GcvB, also increased (~1.8-fold) under microaerobic conditions. Despite the well-documented role of GcvB in downregulating the oppABCDF operon under aerobic conditions [65], our transcriptomic and proteomic data indicate that this sRNA does not inhibit oppABCDF expression under microaerobic conditions. Our finding that levels of aspartate 1-decarboxylase (PanD), a GcvB target previously reported to be involved in pantothenate biosynthesis [66], were decreased under microaerobiosis support the idea that GcvB may downregulate this biosynthetic pathway in oxygen-limited environments. That same study [66] reported other elements of the GcvB targetome, including the csgDEFG operon. It is conceivable that inhibited translation by GcvB-as well as by another sRNA, RprA, also upregulated under microaerobic conditions-likely results in reduced csgDEFG transcript levels. The notion that RprA exerts an active repressional role is supported by its contribution to inhibition of dgcM translation [2,33], which reduces levels of DgcM, a member of the signaling cascade that controls Curli biosynthesis.
It is well documented that SsrA RNA (also known as tmRNA) is involved in ribosomeassociated quality control. It releases stalling ribosomes from truncated mRNAs lacking stop codons through a tmRNA-mediated mechanism, termed trans-translation, which includes peptide-tagging of incompletely synthesized polypeptides for degradation [67]. SsrA processing that leads to functional SsrA-tmRNA translation activity requires cleavage of the SsrA precursor by RNase E [68]. We found that SsrA was the most abundant sRNA under both aerobic and microaerobic growth conditions ( Figure 3B,C). Moreover, its abundance was 1.8 times higher under microaerobiosis ( Figure 3A), suggesting that protein quality control plays an important role in conferring cellular fitness under low-oxygen conditions.

Bacterial Strain and Growth Conditions
To prepare subcultures from fresh overnight cultures, E. coli K-12 strain MG1655 was grown overnight at 37 • C for 16 ) 2 ). The 16 h fresh overnight culture was diluted into 750 mL fresh M9 medium (to OD 460 = 0.04 to 0.05) in a 1 L fermentation vessel chamber (Winpact Parallel Fermentation System FS-05-220, Saratoga, CA, USA). For aerobic culture conditions, air was continuously pumped into the chamber at 0.4 LPM (liters per min). For microaerobic culture conditions, oxygen levels in the chamber containing fresh M9 medium were initially decreased by supplying N 2 at 0.4 LPM until dissolved oxygen (DO) reached 0. N 2 was then pumped for a further 30 min before diluting the overnight culture in the chamber to OD 460 = 0.04 to 0.05, and finally the N 2 supply was turned off. The chamber was completely sealed and the culture was allowed to grow under microaerobic conditions without any additional gas supply. Aerobic and microaerobic cultures were both grown at 200 rpm, 37 • C, and maintained at pH 7.0 by automatic titration with sterile 1 M KOH. Cultures were harvested at OD 460 = 0.5 to 0.6 for transcriptomic or proteomic analyses.
For RNA half-life analysis, rifampicin at a final concentration of 50 mg/mL was used to inhibit new RNA synthesis, and an aliquot of the culture was collected for halflife determinations. In brief, 42 mL of culture from multiple biological replicates for each time-point was collected into 50 mL tubes with 7 mL (1/6 volume) of ice-cold stop solution (5% phenol and 95% ethanol (v/v)) for RNA isolation (see details below). Bacterial pellets were harvested following centrifugation at 4000× g, 4 • C for 15 min, and stored at −80 • C before use. We prepared 5 and 10 biological repeats for aerobic and microaerobic conditions, respectively.

RNA Isolation
Total RNA was extracted as described previously [17]. In brief, bacterial pellets were resuspended in 4 mL KJ medium (50 mM glucose, 25 mM Tris-HCl pH 8.0, 10 mM EDTA pH 8.0, 100 mM NaCl), lysed by placing into boiling 4 mL buffer (0.2 M NaCl, 20 mM Tris-HCl pH 7.5, 40 mM EDTA, 0.5% SDS), and boiled in a boiling water bath for 45 sec before adding 4 mL of acidic phenol (pH = 4.5) and mixing gently by slowly inverting the tube~20 times. Total RNA was extracted in aqueous phase by centrifugation at 4000 g, 4 • C for 1 h. The RNA was precipitated in 1 volume of isopropanol and 1/10 volume of 3 M sodium acetate (pH 7.8) at −20 • C. All RNA samples were maintained in isopropanol at −20 • C before use. When RNA isolation was performed on aliquots, the same volume of culture from the same batch of biological replicates was used for protein isolation, Western blot analysis, and proteomic analysis.

RNA-seq Analysis
Total RNA in isopropanol was pelleted by centrifugation at 15,000 rpm for 15 min, washed with 70% ethanol, and then resuspended in DEPC H 2 O. DNase I (Thermo Fisher antibody (produced in our laboratory, rabbit anti-His-tag Hfq purified protein) was diluted at 1:10,000, whereas α-GAPDH monoclonal antibody (SignalChem, Richmond, VA, USA) was diluted at 1:1000. Membranes were washed once with 1 × TBST buffer for 5 min, before they were incubated with secondary antibody. Secondary mouse-or rabbit-HRP antibody (GE Healthcare, Chicago, IL, USA) was diluted at 1:10,000 and incubated for 1 h at room temperature. Finally, membranes were washed three times with 1 × TBST buffer for 5 min at room temperature. Signals were detected by means of an ECL Western Blotting Detection kit (GE Healthcare) and captured by a BioSpectrum 815 system (UVP).

Protein Extraction, iTRAQ Labeling, and LC-MS/MS
Sample proteins (O-1 and O-2 for aerobiosis; N2 and N-3 for microaerobiosis) were purified using a commercialized B-PER ® Bacterial Protein Extraction Reagent (Thermo Scientific) according to the manufacturer's instructions. Extracted protein digestion, iTRAQ labeling of peptides and subsequent LC-MS/MS analysis, iTARQ signal normalization, and protein quantitation were performed by the core services of the Mass Spectrometry Facility (Academia Sinica, Taipei, Taiwan), as previously described [69]. Briefly, the extracted proteins were subjected to trypsin digestion at 37 • C overnight, lyophilized, and then reconstituted in iTRAQ reaction buffer. Equal amounts of peptides from each sample were individually labeled by adding iTRAQ Reagent 113, iTRAQ Reagent 114, iTRAQ Reagent 117, or iTRAQ Reagent 118, and vortexing the resulting mixtures at room temperature for 1 h. The iTRAQ-labeled peptides were then desalted using a ZipTip concentrator (Merck, Kenilworth, NJ, USA) and mixed. The multiplexed samples were further analyzed by LC-ESI-Q-TOF mass spectrometry. The resulting MS/MS spectra were exported using Mascot Distiller with default parameters. Mascot search results that satisfied the standard criteria [69] revealed the qualified peptides. Their normalized iTRAQ signals were used to quantify the relative abundances of each peptide as well as their fold-changes.