Overlapping stimulons arising in response to divergent stresses in Escherichia coli

Cellular responses to stress can cause a similar change in some facets of fitness even if the stresses are different. Lactose as a sole carbon source for Escherichia coli is an established example: too little causes starvation while excessive lactose import causes toxicity as a side-effect. In an E. coli strain that is robust to osmotic and ionic differences in growth media, B REL606, the rate of antibiotic-tolerant persister formation is elevated in both starvation-inducing and toxicity-inducing concentrations of lactose in comparison to less stressful intermediate concentrations. Such similarities between starvation and toxification raise the question of how much the global stress response stimulon differs between them. We hypothesized that a common stress response is conserved between the two conditions, but that a previously shown threshold driving growth rate heterogeneity in a lactose-toxifying medium would reveal that the growing fraction of cells in that medium to be missing key stress responses that curb growth. To test this, we performed RNA-seq in three representative conditions for differential expression analysis. In comparison to nominally unstressed cultures, both stress conditions showed global shifts in gene expression, with informative similarities and differences. Functional analysis of pathways, gene ontology terms, and clusters of orthogonal groups revealed signatures of overflow metabolism, membrane component shifts, and altered cytosolic and periplasmic contents in toxified cultures. Starving cultures showed an increased tendency toward stringent response-like regulatory signatures. Along with other emerging evidence, our results show multiple possible pathways to stress responses, persistence, and possibly other phenotypes. These results suggest a set of overlapping responses that drives emergence of stress-tolerant phenotypes in diverse conditions.


Introduction
Fitness and survival of a single-celled species across diverse environments incur classic trade-offs between metabolic costs and improvement of survival probability: an effective response in one environment may incur penalties in another. Mesophilic bacteria have adapted to life in intermediate conditions where several dimensions of the cellular environment can readily cause stress, such as excessively high or low temperatures, acidity, nutrients, osmotic pressure, chemical concentrations, and so on. In the lifecycle of an enteric bacterium, drastic changes arise within and between the host and in natura conditions. One such stress is antibiotics. The widespread use of antibiotics to treat livestock for enhanced meat production as well as the rise of antibiotic-resistant bacterial strains in medical contexts lends the question practical importance. Genetic resistance may also have an important relationship with transient antibiotic There is an envelope of survivable conditions in mesophilic bacteria. The volume labeled "Growth Rate-Optimized" denotes conditions to which the microbes are well-adapted for fast growth. In the peripheral volume, "Stress Survival-Optimized" (e.g., the stringent response and bet-hedging), the combination of environment and response improves colony survival but not growth rate. b. The mechanism to perturb fitness in this study is titration of lactose as a sole carbon source in minimal media. Lactose processing involves intrinsic physiological costs and may cause downstream toxicity in certain conditions via the Leloir pathway (galactose processing; not shown). tolerance: persistence may allow resistance to evolve more quickly (3,4), though this hypothesis has been contested. Multiple mechanisms induce persister formation (5)(6)(7)(8)(9), including starvation or loss of metabolic activity (10)(11)(12)(13)(14).
We previously demonstrated that excessively high or low lactose concentrations (as a sole carbon source) can predispose populations of bacteria derived from a strain of Escherichia coli, B REL606, to lowered death rates in antibiotics (15). Varying the concentration of lactose in a minimal medium drives a non-monotonic relationship with the exponential growth rate in culture, with a fast-growth plateau at intermediate concentrations (15). Lactose has established toxic effects on E. coli cultures, often attributed to membrane depolarization via excessive proton symport with lactose through the permease LacY, a member of the major facilitator superfamily of permeases (16)(17)(18). In B REL606, toxic lactose levels create a heterogeneous population dynamic with a chance of fast-growing cells to enter growth arrest, yet both starving and toxified cultures exhibit increased stress tolerance (15, 19; Fig. 2). Growth-arrested cells represent a persister-prone subpopulation in both conditions such that the toxified culture has the effect of a bet-hedging dynamic with average population growth rate sacrificed in favor of stress tolerance in a minority of the population. Does the global transcriptional program of stress responses have overall similarities between these conditions, or are key alternate pathways activated? The mechanisms for these conflicting stresses to attain a similar phenotype are unknown. We interrogated this system with the bulk (population-level) RNA-seq under the hypothesis that a common stress response is conserved between starvation and toxicity. While making use of bioinformatic, biochemical, and physiological data about specific mechanisms, we avoided presuming that any particular mechanism was at work beyond the proven phenomena arising in these culture conditions. We    ed with ampicillin demonstrates a lowered rate of death in starving and toxified cultures (Fig. 2). To confirm this, we fitted the data to a mixed linear-exponential statistical model in logarithmic coordinates on the y axis: . All three parameters were fit with low error ‫ݎ(‬ ଶ 0 . 9 9 ). Taking the time derivative of the statistical model revealed the estimated rate of killing for each culture condition ( Fig. 2b): highest in non-stressful conditions. Thus, both the starving and toxified cultures are prone to produce higher levels of antibiotic tolerance than less stressful intermediate conditions.

Differential gene expression analysis reveals growth-mediated shifts in lac operon expression, toxin-antitoxin responses, and global regulator responses
To analyze gene expression profiles, we purified total RNA from early-mid-exponential cell culture in different lactose concentrations. We mapped sequencing reads to the reference genome E. coli B REL606 NC_012967.1 (20) with kallisto (21) and subsequently analyzed the count data using DESeq2 in R (22) with subsequent processing in Python. Setting the moderate lactose concentration (2.5 mg/ml) condition as the reference, we defined differentially expressed genes (DEGs) in starvation (lactose conc. 0.1 mg/ml) and toxified (lactose conc. 50 mg/ml) conditions as genes with an adjusted p-value (false detection rate; FDR) of below 0.05 and a log 2 fold change (LFC) greater than 1 ( Figure 3).
In Figure 3a, the base mean is the normalized mean expression level for each gene in all replicates in the culture condition. The standard errors from the shrunken log 2 fold change to corresponding maximumlikelihood estimates are well controlled (Figure 3c). They aligned well for a proper size factor calculation.
Overlaying the differential expression log 2 fold change shows that the starving cultures had wider log 2 fold changes, suggesting a more extensive regulon, and perhaps more severe stress, compared to toxified cultures. This result is consistent with our previous results (15) and the model shown in Figure 1, where cells have a higher death rate under low lactose compared to toxic lactose conditions.  activated promoters, and effects on mRNA concentrations by differential loss from either growth or a shift in mRNA degradation. The particular strain we used also expresses a recombinant GFP constitutively, which we have previously exploited as an optical proxy for growth rate because slower-growing cells accumulate GFP and become brighter (15,19).
Following the steps of lactose processing is informative to see the difference between molecular gene regulation and the consequences of fitness, growth, and selection on gene expression in culture. LacY is the lactose permease, transferring extracellular lactose across the inner membrane into the cytoplasm (Fig.   1b). The starved cells contain substantially increased expression of LacY mRNA, suggesting that surviving cells in the culture have increased lactose uptake. In toxified cultures, lacY is not differentially expressed, suggesting relaxed selection for high lactose permease levels. The gene lacZ encodes β galactosidase, the catabolic enzyme that directly assimilates lactose into the metabolic network. lacZ is downregulated in toxified cultures, which may reflect the fact that most cells in the toxified culture randomly avoided toxicity by lowered lactose assimilation. Galactose degradation (the Leloir pathway) is directly downstream of lactose degradation by β -galactosidase, and UDP-galactose and galactose-phosphate intermediates in the Leloir pathway can cause stress and reduce cell growth (24,25). Differential expression of lacA has unclear functional consequences and may simply reflect co-expression with lacZ and lacY from the operon. We next explored differential expression of global regulators of growth rate: toxin-antitoxin systems ( Fig. 5b) and transcriptional regulators (Fig. 5c) (26). In starving cultures, 6 global regulators out of 9 are downregulated, including stress response regulator genes dnaKJ, DNA protectionrelated regulator genes hupAB, a folate-dependent enzyme inhibitor gene ygfA, and persister cell formation regulator gene yigB. Though the recombination-related regulator genes ihfAB are upregulated, suggesting a possible higher mutation rate, dksA is also upregulated; the latter's product can promote DNA recombination repair (27). The product of dksA is also an rRNA repressor that mediates stress responses and inhibits DnaK. Thus, we observe a pattern of gene regulation with counterbalancing effects that may be consistent with phenotypic decisions being made post-translationally. In toxified cultures, ihfB and the chemotaxis signaling genes cheZ and cheY were upregulated.  Our results show a significant decrease in σ 70 (rpoD) in toxified, but not starving, cultures (Fig. 6c,d).
This σ factor initiates transcription in a set of approximately 1723 genes involving cell proliferation-  As starvation and toxicity are both intertwined with metabolism, we performed metabolic pathway gene set enrichment analysis in peripheral (Section 2.4) and central (Section 2.5) metabolism. We used the normalized enrichment score (NES), a statistic that reflects the degree that a pathway is overrepresented at the top or bottom of a ranked list of genes (Methods Section 4.5). An overview of our pathway enrichment analysis is presented in Figure S2.

Pathway regulatory similarities between starving and toxified cultures
Gene set enrichment analysis calculates trends for defined gene sets. The similarity of pathway regulation was calculated using rank correlations (Spearman's ρ ) with differential gene expression profiles (see Methods Section 4.5.1). Figure S3 shows the correlation between pathway regulation in starvation and toxicity. Each node in Figure S3 represents a pathway, with a total of 362 pathways annotated in E. coli B REL606. Nodes are connected by metabolite edges as annotated in EcoCyc. We found that most metabolic pathways are either similarly regulated or have no correlation between treatments. We observed conserved regulatory pathways for persister-prone phenotypes, including central metabolic pathways such as the pentose phosphate pathway, chorismate biosynthesis I, glycolysis I, and the superpathway of glyoxylate bypass and TCA. Glycolysis is the hub node for the pathway network with similarly regulated pathways linking to glycolysis with degree 1 node distances. Statistically significant oppositely regulated pathways are rare.
The only oppositely regulated pathway with statistical significance is tetrahydrofolate (FH 4 ) biosynthesis, which is upregulated in starvation but downregulated in toxicity. FH 4 biosynthesis produces vitamin B9 (folic acid), a cofactor leading to the biosynthesis of methionine, purines, thymidylate and pantothenate.

Common enriched pathways in persister-prone phenotypes
Among 362 enriched pathways, 44 were upregulated in starvation and toxification, 69 were downregulated in starvation and toxification, 14 were uniquely upregulated in starvation, 31 pathways were uniquely downregulated in starvation, 12 pathways were uniquely upregulated in toxification, and 25 pathways were uniquely downregulated in toxification. Detailed pathway enrichment analysis results can be found in Supplementary Tables S1-S6. Figure 7 shows the pathway hierarchy structure based on EcoCyc annotations. The directed edges point to the parent pathway which is often consisted of multiple pathways to form superpathways. In Figure 7 we mapped the normalized enrichment score calculated with FGSEA (34) onto the pathway hierarchy graph. The result is many clustered enriched pathways. The top 3 clusters in Figure 7 are (1) (36). Biosynthesis of ubiquinone is reported to accumulate pathway intermediates (37,38) with the effect of improving osmotic stress-tolerance.
We found that 17 out of 20 amino acid biosynthesis pathways are downregulated in persister-prone phenotypes accompanying upregulated amino acid degradation pathways (Table S2.). Alanine and arginine biosynthesis are uniquely downregulated in toxification, while asparagine is uniquely downregulated in starvation, implying lower NAD synthesis (39) and likely resulting in lower energy levels. Amino acid starvation can lead to accumulation of uncharged tRNAs that enter the ribosomal A site, halting transla-tion (40). We observe downregulated tRNA charging in both toxicity and starvation. The pyrimidine, purine and pyridine nucleotide synthesis-related pathways are downregulated in both starvation and toxicity except for phosphoribosylpyrophosphate (PRPP) biosynthesis II, which provides PRPP as a pivotal metabolic precursor to pyrimidine, purine, and pyridine nucleotide synthesis. Other downregulated pathways include sulfate assimilation downstream of cysteine biosynthesis, lipid A, lipopolysaccharide (LPS) and peptidoglycan biosynthesis, spermidine biosynthesis, UDP-glucose biosynthesis, and cytochrome b o oxidase electron transfer pathways from proline succinate.
Upregulated pathways have a sparser network structure (Table S1). Among 44 commonly upregulated pathways, 30 are degradation pathways targeting amino acids, other carbon sources, and electron transferrelated metabolites, such as L-ascorbate (vitamin C) and putrescine. Broadly upregulated nutrient assimilation pathways reflect supervening carbon source uptake and balancing energetic electron transfer. Other upregulated pathways enhance cellular fitness through energy metabolism, overflow metabolism, and membrane component regulation. In both persister-prone conditions of this study, the pyruvate oxidation pathway is upregulated; this pathway generates the transmembrane potential for constructing respiratory chains consisting of pyruvate oxidase, ubiquinone-8, and the cytochrome bd complex (41). Thiamine diphosphate biosynthesis I produces cofactor vitamin B1, which plays a fundamental role in energy metabolism. Trehalose pathway upregulation affects the osmotic stress response, notably present in both stress conditions. In response to acetate and osmotic pressure, the arsenate efflux pump, GDAR system, and aerobic utilization of acetate are upregulated. GDP-L-fucose biosynthesis I produces LPS components in the membrane, and membrane structure may be further stabilized by enrichment in putrescine biosynthesis III, which is the precursor for spermidine biosynthesis.

Uniquely enriched pathways in starvation
Pathways uniquely upregulated in starvation (Table S3) include nutrient assimilation from metabolites mannitol, D-arabinose, acetoin, glycerol, glycolate, glyoxylate, trehalose, and fatty acids. The arginine-dependent acid resistance pathway was also observed.
Though downregulation (Table S4)  Finally, asparagine biosynthesis pathways are uniquely downregulated in starvation, with unclear functional significance. (Table S5) Two pathways forming the E. coli anaerobic respiratory chain are upregulated: the formate-todimethyl sulfoxide electron transfer pathway and nitrate reduction III, implying increased proton-motive force across the cytoplasmic membrane (42)(43)(44)(45). We observed upregulation of taurine degradation IV, which provides alternative sulfur under sulfate starvation induced by cysteine starvation (we note that our The superpathway of fatty acid biosynthesis initiation is uniquely downregulated in toxified conditions (Table S6). Related pathways such as fatty acid elongation, palmitoleate biosynthesis, and cisvaccenate biosynthesis are also uniquely downregulated. Fatty acids are key building blocks for phospholipid components of cell membranes and are determinants of intracellular communication where palmitoleate is a common unsaturated fatty acid. Farnesyl pyrophosphate (FPP) biosynthesis is downregulated, possibly implying less membrane attachment in posttranslational modifications.

Gene set enrichment analysis of central metabolism
To complete the analysis of regulated metabolic pathways in starvation and toxicity, we mapped differential expression to central metabolism (Fig. 9). As the culture conditions in this study demand that the initial entry to central metabolism occurs via lactose degradation, we define central metabolism herein to be composed of the lactose degradation pathway, pentose phosphate pathway, glycolysis, Entner-Doudoroff shunt, TCA cycle, and the glyoxylate bypass. As galactose is a direct product of lactose degradation, we included the Leloir pathway (galactose degradation I) as well.
In starvation, the glycolysis pathway is mostly downregulated except for one gene encoding a phosphate transfer reaction. Components of glycolysis are more consistently upregulated in toxicity. The Entner-Doudoroff shunt linking glycolysis to the TCA cycle is regulated oppositely between starving and toxified cultures. The shunt is upregulated in starvation, possibly providing paths to alternative metabolites. Though the TCA cycle is largely downregulated in both conditions, the glyoxylate cycle is upregulated by various degrees in both starvation and toxicity. The glyoxylate cycle is a bypass for the TCA cycle to skip steps that remove CO 2 , thus conserving carbon intermediates during growth (47)(48)(49).
The pathway is repressed during growth on glucose, and induced by growth on acetate, the byproduct of overflow metabolism.

Discussion
Variation in fitness-relevant environmental properties affects cellular gene expression patterns in quantitative and qualitative ways. Our previous demonstration that an excess of a required nutrient drives the formation of antibiotic tolerance (15) provided an opportunity to re-evaluate the nature of integrated microbial stress responses. This phenomenon appears to arise from the robust nature of the cell wall in B strains of E. coli, which permits survival in adverse osmotic conditions (50). We were specifically able to enrich antibiotic-tolerant cells in conditions both below and in excess of optimal concentrations for lactose as the sole carbon source (Fig. 1). We created such conditions in batch culture to test hypotheses about the nature of antibiotic tolerance in E. coli: what kind of signal may predispose cells to stress tolerance, and do conditions that could be described as "opposite" of each other from osmotic, nutrient concentration, and media toxicity perspectives induce similar or largely different responses?
We conceptualized analysis of the stress response by analogy to a statistical model, with linearly independent (normalized) stress response components is a high-dimensional vector describing the phenotype. Independence between components may arise in some contexts and not others -we did not presume to comprehensively describe this space in our study, but rather to sample from it for exploration of the specific nature of stress responses. Nonzero interaction terms may arise from several sources, including pleiotropic effects (e.g., temperature or growth rate), co-expression (operon membership, transcriptional co-regulation, etc), and catalysis or synthesis of metabolites with broad effects (amino acids, nucleotides, tRNAs, etc).
Single cell dynamics play an important role in physiology and have a strong effect on lactose-toxified cultures. In the E. coli strain used in this study, a threshold drives toxicity to arrest growth in a subset of cells while subthreshold cells maintain growth. Thus, the results presented here represent an average across the subthreshold growing cells and the growth-arrested toxified subpopulation. The subthreshold subpopulation comprises approximately 10 percent of the entire population (15). This quantity arises from exponential growth of sub-threshold cells causing them to quickly overtake the population.
We exploited extensive mechanistic knowledge about the direct effects of pathways to interpret transcriptional signatures. Targeted experiments to further test other aspects of responses, such as posttranslational effects, are beyond the scope of this study. Our results recapitulated previous interpretations of pathway analysis in stress but added crucial new results that broaden some, and narrow other, key interpretations. This study focused on the nature of the average gene expression in cultures known to be enriched for persister cells, but with no specific selection for them. In a following study, we will present an analysis of persister cells in the starving and toxified conditions that survived antibiotic treatment.

Leloir pathway intermediates
One hypothesis for non-starving persistence is that metabolic toxicity is induced by critical proteomic concentrations. Galactose degradation I (the Leloir pathway) is directly downstream of lactose degradation by β -galactosidase and may accumulate UDP-galactose and galactose-phosphate intermediates causing stress and loss of growth. With high metabolic rates, E. coli (and virtually every other species) undergoes a metabolic shift from primarily aerobic metabolism to incomplete oxidation of metabolites, including ATP synthesis (51). The cause seems to be linked to proteomic optimization, as anaerobic ATP synthesis requires a smaller fraction of the proteome to synthesize an equivalent amount of ATP at the cost of more sugar (52). The smaller proteome allows for higher growth rates due to the reduced size of the necessary metabolome, allowing more transcriptional/translational machinery to be devoted to ribosome synthesis (52,53).
We found that toxic culture conditions that produce persister cells have apparent utilization of the Entner-Doudoroff shunt, which connects pyruvate to phosphoenolpyruvate, thus linking to the glyoxylate cycle and downregulation of the enzymes responsible for oxaloacetate and acetyl-CoA entering the citric acid cycle, consistent with cells undergoing overflow metabolism. As the only carbon source initially pre-sent in the medium is lactose, GalE fluctuation may lead to UDP-Galactose toxicity (24), and indeed, galE is downregulated in high-lactose antibiotic-tolerant cells, but not significantly differentially expressed in untreated cultures.

Different phenotypes may arise from sigma factor competition
Regulation of stress responses in bacteria is generally considered to be robust to a variety of conditions. Despite being exposed to opposite stresses, we showed that differential expression of genes in starvation and toxicity is positively correlated. This may be due to leaky expression of genes in the different loci on the genome or generalized stress responses. As expected, regulated differentially expressed genes have wider fold change distributions compared to constitutively expressed genes and more genes are regulated in starvation than toxicity. We found that the global proliferation regulator

Non-transcriptional and undetected factors at play in stress-responsive gene expression
Survivorship bias, as opposed to targeted transcriptional regulation, may play a role in the observed transcriptome profiles obtained from stressed cells. Noisy gene expression dynamics creates a distribution of cell states, including some in sub-optimal transcriptional states. With exponential population growth, fitter cells quickly overtake the population and dominate the measured transcriptional signal. Furthermore, it was impractical to examine every possible transcriptional function. Some that we did not emphasize, such as prophage induction, may play an uncharcterized role in distinguishing responses to di-vergent stresses. Thus, studies of subpopulations and deeper targeted analysis may further clarify the overall similarities between divergent stresses.

Summary
The relationship between bacterial responses to divergent stresses contains four primary components: (i-ii) specific responses to each of the two stresses, (iii) a common stress-responsive regulon, and (iv) noise or functionally irrelevant responses. Overlapping responses include transcripts relating to nutrient assimilation, tRNA charging, and utilization of the glyoxylate bypass while condition-specific responses make sense in terms of unique properties of the stress in question (Fig. 10).

Persister Enriched RNA-Seq Experiments
E. coli B REL606 lacItransformed with Tn7::PlacO1GFP(KanR) was inoculated in LB medium from a -80 C bacterial stock and grown for 16 hours in a 37 C shaking incubator. The LB culture was then resuspended (1:1000) into 5mL of Davis Minimal medium (DM; Difco) supplemented with thiamine and one of three lactose concentrations (0.1 mg/mL, 2.5 mg/mL, and 50 mg/mL). The cultures were allowed to acclimatize for 24 hours before being resuspended (1:1000) into 5mL of the same Davis Minimal medium and lactose concentration. Cultures were grown long enough to provide enough biomass for RNA-seq after antibiotic treatment (8 hours in 2.5 mg/mL lactose, 10 hours in 50 mg/mL lactose, and 12 hours in 0.1 mg/mL lactose). After the initial growth phase, 1.5 mL cell culture aliquots were subjected to RNA isolation according to the following procedure.
Cell cultures were pelleted in a microcentrifuge (10,000 G for 2 minutes), washed in PBS buffer twice, resuspended in 500 μ L of RNA-Later (ThermoFisher) and stored in at -20 C for up to one week.
RNA isolation was performed using Direct-Zol (Zymo) and TRIzol reagent (ThermoFisher) and stored in

RNA-seq sequence alignment and genome annotation with Ecocyc and RegulonDB
RNA transcript quantification was performed using kallisto (21)  RegulonDB (55). Lacking extensive annotation of gene regulation in REL606, we used the K-12 MG1655 strain annotation based on gene names and gene product similarities.

Differential Expression Analysis
We used R package DESeq2 (56) for gene differential expression analysis. The RNA transcription quantification data were first clustered to isolate outlier replicates using principal component analysis (PCA) ( Figure S2). Two samples with a high number of missing transcripts were dropped from subsequent analysis. To confirm reproducible outcomes of sample treatment, hierarchical clustering using the Wald significance test was performed on all samples; sample treatment was retrieved perfectly ( Figure   S1).
The DESeq2 pipeline includes size factor estimation, dispersion estimation, and DEG tests. Lowcount RNA quantifications are noisy and may decrease the sensitivity of DEGs detection (57). Size factors were calculated with a subset of control genes: non-regulated genes according to RegulonDB (55) with expression higher than a threshold (10) across all replicates. Setting transcriptome quantification from the moderate lactose condition, we applied the adaptive-T prior shrinkage estimator "apeglm" and used Wald significance tests for detecting DEGs and obtaining the log 2 fold changes (LFC2).

Enrichment analysis by FGSEA
At the time of analysis, there were 368 pathways in EcoCyc for B REL606. 6 pathways that are not apparently regulated by gene products were discarded. Pathway enrichment was analyzed using FGSEA (58). Differentially expressed genes were pre-ranked by their log 2 fold change. Pathway gene sets were defined with reference to the EcoCyc database. The pre-ranked gene data and pathway gene sets were then processed by the fgseaMultilevel function to obtain final results. Minimum gene set size was set to 3 with 200 bootstrap replicates. The threshold for pathway significance was 0.05. Figure 7 are linked by key metabolites and their flow direction. The pathway ontology in Figure 7 is annotated in the Ecocyc database. The pathway link and pathway ontology was extracted from EcoCyc and visualized with Cytoscape (59).

Pathway regulatory mechanism similarity
To calculate the similarity between pathway regulation in different conditions, we used Spearman's ρ rank correlation between the differential expression profile for each pathway. The similarity score for cases where all gene enrichment in a pathway are of the same direction is set to 1.
Metabolic pathway visualization used the Python package Escher (60).