Analysis of Growth Phases of Enterotoxigenic Escherichia coli Reveals a Distinct Transition Phase before Entry into Early Stationary Phase with Shifts in Tryptophan, Fucose, and Putrescine Metabolism and Degradation of Neurotransmitter Precursors.

Enterotoxigenic Escherichia coli (ETEC) is a major cause of diarrhea in children and adults in endemic areas. Gene regulation of ETEC during growth in vitro and in vivo needs to be further evaluated, and here we describe the full transcriptome and metabolome of ETEC during growth from mid-logarithmic growth to early stationary phase in rich medium (LB medium). We identified specific genes and pathways subjected to rapid transient alterations in gene expression and metabolite production during the transition from logarithmic to stationary growth. The transient phase was found to be different from the subsequent induction of early stationary phase-induced genes. The transient phase was characterized by the repression of genes and metabolites involved in organic substance transport. Genes involved in fucose and putrescine metabolism were upregulated, and genes involved in iron transport were repressed. Expression of toxins and colonization factors were not changed, suggesting retained virulence from mid-logarithmic to the start of the stationary phase. Metabolomic analyses showed that the transient phase was characterized by a drop of intracellular amino acids, e.g., l-tyrosine, l-tryptophan, l-phenylalanine, l-leucine, and l-glutamic acid, followed by increased levels at induction of stationary phase. A pathway enrichment analysis of the entire combined transcriptome and metabolome revealed that significant pathways during progression from logarithmic to early stationary phase are involved in the degradation of neurotransmitters aminobutyrate (GABA) and precursors of 5-hydroxytryptamine (serotonin). This work provides a comprehensive framework for further studies on transcriptional and metabolic regulation in pathogenic E. coli. IMPORTANCE We show that E. coli, exemplified by the pathogenic subspecies enterotoxigenic E. coli (ETEC), undergoes a stepwise transcriptional and metabolic transition into the stationary phase. At a specific entry point, E. coli induces activation and repression of specific pathways. This leads to a rapid decrease of intracellular levels of certain amino acids. The resulting metabolic activity leads to an intense but short peak of indole production, suggesting that this is the previously described "indole peak," rapid decrease of intermediate molecules of bacterial neurotransmitters, increased putrescine and fucose uptake, increased glutathione levels, and decreased iron uptake. This specific transient shift in gene expression and metabolome is short-lived and disappears when bacteria enter the early stationary phase. We suggest that these changes mainly prepare bacteria for ceased growth, but based on the pathways involved, we could suggest that this transient phase substantially influences survival and virulence.

Escherichia coli is a facultative anaerobic Gram-negative bacterium that normally inhabits the intestines of mammals and reptiles as a commensal bacterium. Pathogenic E. coli have acquired extrachromosomal genetic properties that enable them to colonize and adhere to the epithelium, thereby delivering toxins or virulence factors that harm the host (1). The virulence factors can either be located on plasmids or inserted in the chromosomes as pathogenicity islands, and pathogenic E. coli can be found in most E. coli phylogroups (2–4). Enterotoxigenic Escherichia coli (ETEC) is characterized by the production of the heat-labile toxin (LT) and/or the heat-stable toxin (ST). In most cases, each bacterium expresses one to three colonization factors that mediate adhesion to the epithelium in the small intestine (5, 6). The toxins and colonization factors of ETEC are mainly encoded on extrachromosomal plasmids that have been acquired by horizontal transfer to ancestral commensal E. coli. Successful combinations of host and plasmid may lead to global transmission of virulent clones with optimal colonization and survival abilities (2, 4).

The molecular events that govern the virulence of ETEC are less well characterized than several other intestinal pathogenic E. coli. Studies on ETEC have so far focused on a few global regulators such as Crp, H-NS, and members of the AraC family such as Rns and CsvR that were shown to regulate the toxins and/or CFs (7–11). However, virulence factors are not only the expression of toxins or colonization factors but can also be expanded to involve the ability to persist host-induced stress or to facilitate spreading or colonization ability. Studies using real-time PCR, microarrays, and RNA-Seq have recently begun to elucidate the global and/or specific transcriptional regulation in ETEC in response to host environmental factors such as bile and glucose and adhesion to epithelial cells (8, 11–15).

In this study, we explored the transcriptional and metabolomic profile of two clinical ETEC isolates belonging to a globally spread linage during growth from logarithmic cell division to early stationary phase to elucidate the effect of growth on transcriptomic profile and virulence gene regulation. We used a multi-omics approach including RNA-seq and mass spectrometric global metabolomics techniques to analyze the global regulation and provide a framework of E. coli genes, transcription factors, and analysis of metabolites involved in different growth phases in LB medium. The results can be applied not only on ETEC but also on other enterobacterial species and provide insights into gene and metabolite shifts at different growth phases.

RESULTS

Expression profiling of two ETEC revealed characteristic transcriptional patterns during transition from log phase to induction of stationary phase.

To profile the ETEC/E. coli transcriptome during bacterial growth transition from mid-exponential to early stationary phase, we performed RNA-seq analysis on RNA isolated from two clinical isolates of ETEC (E1777 and E2265) grown in LB media (single replicates per strain). We sequenced the transcripts expressed after 3 h, 4 h, and 5 h of growth corresponding to mid-log phase (3 h, OD600 = 1.1 to 1.3), late log phase (4 h, OD600 = 2.8 to 3.1), and early stationary phase (5 h, OD600 = 4.6 to 4.8) (Fig. 1).

FIG 1

Growth curve of ETEC strains. Total bacterial counts in CFU (CFU/mL) and optic density (OD600) of the ETEC strains E2265 and E1777 in LB media. Samples were measured every hour for 6 h and after 12 h of incubation.

The two strains express the toxins heat-labile toxin (LT) and heat-stable toxin (STh), and colonization factors coli surface antigens 5 and 6, (CS5+CS6), and belong to the globally distributed L5 ETEC lineage (4). The time points were chosen to reflect the transition from active growth to carbon and nutrient starvation, an environment that enteropathogens face in the human gut and growth in LB (16). First, Illumina sequence reads were assembled and assigned as unigene sequences, translated, and compared with protein databases for annotation. As a result, a total of 4,166 unigenes with expression in at least one sample were detected (Table S1).

Next, in order to identify the differentially expressed genes shared by both isolates E2265 and E1777 and display the dynamic of the transcriptome during the transition from the log phase (3 h) to the early stationary phase (5 h), we performed a differential gene expression analysis by DESeq with a fold change cutoff of 4 (Log2 > 2) and P-value < 0.05. Overall, 617 genes (Table S2, Table S3) showed a significant change (up/down) between two time points, and these results were displayed and clustered in a heatmap in Fig. 2. The following comparisons, 3 h versus 4 h and 3 h versus 5 h, resulted in a number of 486 (257 upregulated and 229 downregulated) and 392 (209 upregulated and 183 downregulated) differently expressed genes (DEG), respectively (Table S3). A total of 495 filtrated differentially expressed genes (DEG) in at least one condition in both bacterial transcriptomes were included to perform the K-means clustering analysis of the DEG heatmap (Fig. 2a; Table S4) resulting in four clusters of specific gene expression patterns (cluster I to IV). In parallel, to further determine the most significant differential gene expression patterns shared among genes with similar dynamics in transcriptional changes over time (Fig. S1), we employed the short time-series expression miner (STEM) clustering method (17).

FIG 2

Transcriptomic response of ETEC during bacterial growth transition from mid-exponential to early stationary phase. (a) Heatmap of the differential expressions of the two strains E1777 and E2265 after 3, 4, and 5 h of growth in LB medium. K-mers analysis indicated four main clusters (I-IV). (b) STEM analysis identifying the most common patterns of gene expression. (c) GO biological process enrichment analysis of significant genes per cluster and fold enrichment. (d) Examples of genes with transient up or downregulation.

STEM analysis identified six gene expression profiles, of which four (labeled in red, white, green, and blue) enclosed the majority of the DEGs and matched with the profiles obtained by k-means clustering analysis. STEM clusters 1, 6, and 9 (Fig. 2b) that corresponds to cluster I, III, and IV (Fig. 2a) had a statistically significantly higher number of genes assigned using a permutation test (P < 0.05). Both methods highlight four distinctive temporal dynamics of the transcriptome in response to the transition from log to early stationary phase (Fig. 2a and b, Table S5). In summary, clusters 1 and III were the largest, enclosing 166 and 110 genes, respectively. Cluster I showed a decreasing expression toward the entry to the stationary phase (3 h to 4 h to 5 h), while cluster III showed the opposite trend. Interestingly, clusters II and IV included 19 and 74 genes, respectively, with a significant transient down or upregulation at 4 h compared with 3 and 5 h. Therefore, the data indicate that a specific transient phase in gene regulation occurs after 4 h and OD600 around 3 h when ETEC enters the late log phase and starts transit into the stationary phase.

In order to gain biological insights from the temporal transcriptomic dynamics, we performed gene ontology (GO) enrichment analysis for biological processes and metabolic pathways of the significantly expressed genes from each cluster (Fig. 2c; Table S6). Significant gene enrichment (FDR < 0.05) of cluster I indicated downregulation of genes involved in the de novo purine biosynthesis pathways such as “de novo” IMP biosynthetic processes as well as nitrogen metabolic pathways, and protein targeting (intracellular protein transport) when growth slowed down upon entry to late log phase and early stationary phase. In contrast, cluster III, which represents a progressive transcriptional gene activation toward the stationary phase, was significantly enriched in glyoxylate metabolism, tricarboxylic acid cycle (TCA) cycle, carboxylic acid metabolic processes, aerobic respiration, and small molecule metabolism. The results indicate that the cells were growing and had an active TCA cycle with increasing expression of glyoxylate bypass genes aceABK at this time point. Organic substance transport, which involves the movement of organic substances that contains carbon in, out, or within a cell, was identified as the only significantly enriched biological process among genes of cluster IV. No significant enrichment was identified for cluster II. Even though clusters II and IV represent a minor proportion of the DEG data set, they exhibited an interesting transient transcriptional response prior to switching into the stationary phase.

Transient transcriptomic activation of putrescine and fucose utilization and reduction of iron transport prior to entry into stationary phase.

To increase our understanding of the biological role of the transiently altered genes (4 h of growth at the beginning of the stationary phase) identified in clusters II and IV reads mapped to the gene of each temporal gene expression were plotted in Fig. 2d. The trendline confirmed the down (cluster II) and up (cluster IV) regulation of these genes during the transient shift to stationary phase. Among the activated genes at 4 h, we identified an overrepresentation of genes from the fuc operon, i.e., fucI (L-fucose isomerase) and fucA (L-fucose 1-phosphate aldolase) and fucU (L-fucus mutarotase) with an approximately 6-fold increase in expression. Another set of genes involved in exploiting alternative nutrient sources was the putrescine pathway. Like the fucose operon, the Puu-operon (putrescine utilization pathway) genes that degrade putrescine to GABA via γ-glutamylated intermediates were highly activated at 4 h. Genes involved in the putrescine uptake system, i.e., potI and ydcU showed the same pattern. Another gene included in this cluster was tnaA, which encodes tryptophanase, responsible for indole production from l-tryptophan. Its expression was activated 20-fold at 4 h compared with the mid-log growth phase at 3 h and reduced 6-fold at 5 h (Table S4). In addition, increased expression of the fad-operon (fadH, fadM) and the dpp-operon in charge of dipeptide transport were evident.

Downregulation of genes in cluster II included exbB and fhuA involved in siderophore-mediated iron transport, yncE, a DNA binding protein involved in iron metabolism, and bfd bacterioferritin-associated ferredoxin; hence, cluster II suggests that a rapid downregulation of genes involved in iron metabolism occurs transiently before entry into early stationary phase.

Expression of ETEC virulence factors is slightly higher during exponential growth.

Both strains analyzed in this study expressed genes encoding the enterotoxins: eltAB (LT) and estA (STh), and colonization factor operons: csfABCDEF (CS5) and cssABCD (CS6) (18). The transcriptome data showed that eltAB expression was higher than estA at all three time points. In comparison with estA expression, eltAB was 3- to 4-fold higher at 3 h and 2- to 4-fold and 2-fold at 4 h and 5 h, respectively. Expression of eltAB and the CS5 encoding cfs-operon, estA, and the CS6 had a trend of gradual downregulation over time (Table S5). The expression patterns were similar between the strains. Although the expression levels were not significantly changed in this study, the findings are confirmed by our previous work (8). Genome analysis showed that both strains also express additional virulence genes: cexE, clyA, eatA, ecpA, and fimH (18). The expression levels of these genes were not significantly changed.

Global view of the intracellular and secreted metabolome of ETEC growth phases.

Because our transcriptomic data provided a framework of how metabolic pathways were altered during the transition from freely available nutrients to a more restricted environment, we wanted to characterize the impact of nutrient depletion on the bacterial metabolome. We collected pellets (intracellular metabolites) and supernatants (secreted metabolites) of E2265 and E1777 bacteria at 3 h, 4 h, and 5 h of growth in LB broth medium. Each bacterial growth curve was independently repeated 3 times, yielding a total of 36 samples that were used to perform a nontargeted GC-MS-based metabolomics approach. Approximately 2,000 putative intracellular and secreted metabolites were detected in both isolates, of which 288 metabolites were successfully identified (Table S7). Principal-component analysis (PCA) (Fig. 3a) of all samples revealed a profound clustering of samples from the intracellular and secreted metabolomes. The analysis also indicated that the intracellular metabolomes of both E2265 and E1777 are more similar than their respective secreted metabolome. The PCA did not show any variation between metabolomes of the strains per time point (Fig. 3a); however, the hierarchical clustering heatmap showed changes in the abundance of several metabolites along with the time points (Fig. S2).

FIG 3

Intracellular and secreted profile of the metabolic response to ETEC during bacterial growth transition from mid-exponential to early stationary phase. (a) PCA plot generated from all metabolites of different samples. (b) Venn diagrams of number of significant metabolites (padj < 0.05 −2 > Log2Foldchange <2) found in the bacteria (intracellular) or medium (secreted). (c) Heatmap representation of the 35 differentially changed metabolites at any time point. The metabolite abundance per time point was used to generate the heatmap, and Z-score was scaled for visualization. (d) Pathway enrichment analysis of the intracellular and secreted significant metabolites. Non-significant pathways were colored gray. (e) Metabolite-metabolite correlation analysis shows positive correlations in dark purple and negative correlations in yellow. Secreted and intracellular metabolites were marked with orange and gray dots, respectively.

Intracellular and secreted metabolome show unique metabolomic shifts.

Because we were interested in studying the differences in the metabolome during bacterial growth transition mid-exponential to early stationary phase, we used multiple t-tests to compare the metabolite abundance between two time points, i.e., 3 h versus 4 h or 3 h versus 5 h (Table S8, Table S9). The threshold for significant changes was log2 fold change > 2 with a padj < 0.05. As shown in Table 1, Fig. 3b, and Fig. S3, 35 metabolites significantly changed their intracellular and/or extracellular concentration at any time point. Specifically, 20 intracellular and 15 secreted metabolites were significantly altered.

TABLE 1

Summary of significantly altered intracellular and secreted metabolites of ETEC

Metabolite	Classification	HMDB ID	Intracellular				Secreted				Pathways
			E1777		E2265		E1777		E2265		Aminoacyl-tRNA biosynthesisb	Butanoate metabolism	Ala, Asp, and Glumetabolism	Arg and Pro metabolism	Glycolys is/Gluconeogenesis	Glutathione metabolism	Purine metabolism	Gly, Ser, and Thr metabolism	Methano metabolism	Nicotine/Nicotinamide metabolism	Vitamin B6 metabolism	Thiamine metabolism	phenylalanine metabolism	Pyrimidine metabolism	Muropeptide degradaation	Peptidoglycan biosynthesis	NF
			3 h vs 4 h	3 h vs 5 h	3 h vs 4 h	3 h vs 5 h	3 h vs 4 h	3 h vs 5 h	3 h vs 4 h	3 h vs 5 h
L-Glutamic acid	Amino acids and peptides	HMDB0000148	−2.24	−0.73	−1.76	−1.15	−1.59	−2.19	−0.87	−2.46	X	X	X	X
L-Tyrosine	Amino acids and peptides	HMDB0000158	−2.47	−0.11	−1.58	−0.99	−0.05	−0.15	0.21	0.25	X
Succinic acid	TCA acids	HMDB0000254	−3.37	−0.36	−2.01	−2.26	−0.51	−0.43	−2.43	−2.47		X	X
L-Arginine	Amino acids and peptides	HMDB0000517	−3.09	−0.48	−2.37	−1.79	−0.26	−0.03	0.15	0.33	X			X
L-Leucine	Amino acids and peptides	HMDB0000687	−2.44	−0.10	−0.41	0.27	0.19	0.21	0.43	0.45	X
L-Phenylalanine	Amino acids and peptides	HMDB0000159	−2.13	−0.03	−1.27	−0.79	0.04	−0.02	0.01	0.09	X
L-Tryptophan	Amino acids and peptides	HMDB0000929	−2.76	−0.04	−1.48	−1.12	0.11	−0.04	−0.10	−0.06	X							X
N-Acetyl-L-aspartic acid	Amino acids and peptides	HMDB0000812	−3.04	−3.30	0.59	−0.17	0.01	0.12	0.02	0.43			X
S-Adenosylmethionine	Glycosyl compounds	HMDB0001185	1.46	2.50	−0.05	0.49	0.00	0.00	0.00	0.00				X
D-Glucose	Monosaccharides	HMDB0000122	−3.27	−0.35	−1.94	−1.40	−0.03	0.01	0.07	0.02					X
Glutathione	Amino acids and peptides	HMDB0000125	3.02	2.32	3.56	4.27	0.00	0.00	0.00	0.00						X
Oxidized glutathione	Amino acids and peptides	HMDB0003337	0.87	0.71	1.69	2.39	0.40	0.48	1.17	1.29						X
GMP	Purines	HMDB0001397	0.23	1.11	1.60	2.13	0.00	0.00	0.00	0.00							X
Xanthine	Purines	HMDB0000292	−2.86	−0.03	1.14	0.76	0.06	0.11	0	0							X
N-methyl-L-glutamic Acid	Amino acids and peptides	HMDB0062660	2.02	2.48	0.00	0.00	0.63	0.58	0.00	0.00									X
NAD	Nicotinamides	HMDB0000902	0.48	−0.05	0.69	0.62	0.01	1.14	3.17	2.59										X
Pyridoxal 5′-phosphate	Pyridine carboxaldehydes	HMDB0001491	2.21	1.06	2.71	3.33	0.00	0.00	0.00	0.00											X
Thiamine	Pyrimidines	HMDB0000235	−2.46	−0.23	−1.29	−0.99	0.07	0.06	0.09	0.08												X
trans-Cinnamic acid	Cinnamic acids	HMDB0000930	−0.52	−0.72	0.15	2.19	0.16	0.05	−0.03	0.28													X
Uridine	Pyrimidines	HMDB0000296	−2.02	−0.48	−0.91	−2.35	0.03	0.02	0.00	0.00														X
Leucyl-Alanine	Amino acids and peptides	HMDB0028922	−3.14	0.02	−1.21	−1.27	0.57	0.61	0	0																	X
N-Acetylcadaverine	Carboxylic acids	HMDB0002284	0.30	2.63	1.25	3.11	1.68	0.76	0.00	0.00																	X
Desaminotyrosine	Phenylpropanoids	HMDB0002199	−0.49	−0.64	0.31	2.50	0.29	0.17	0.01	0.39																	X
5-Aminopentanoic acid	Fatty acids	HMDB0003355	0	0	0.71	1.10	0	0	−0.85	−2.35				X
Adenosine	Purines	HMDB0000050	−0.49	−0.03	−0.31	−0.82	−0.99	0.03	−2.51	−3.14							X
Alanyl-Phenylalanine	Amino acids and peptides	HMDB0028694	−4.19	−0.06	−2.01	−2.00	−2.57	−1.49	−3.32	−4.60																	X
Betaine	Amino acids and peptides	HMDB0000043	0	0	0.17	0.09	0	0	−1.28	−2.20								X
L-Homoserine	Amino acids and peptides	HMDB0000719	0	0	−2.70	−1.72	−0.29	0.00	−0.69	−2.20								X
Alanyl-Glutamine	Amino acids and peptides	HMDB0028685	−3.62	−0.26	−2.45	−2.43	−2.14	−1.12	−1.13	−3.02															X
D-Alanine	Amino acids and peptides	HMDB0001310	−2.79	−0.40	−1.24	−1.40	−0.32	−0.02	−1.19	−3.05																X
Spermidine	Amines	HMDB0001257	0	0	0	0	0	0	−0.81	−2.21				X
GABA	Amino acids and peptides	HMDB0000112	−2.35	−0.24	−1.77	−1.31	−2.90	−1.35	−0.82	−2.29		X	X	X
Maleic acid	Fatty acids	HMDB0000176	−2.52	−0.71	−1.20	−0.94	0.34	0.52	−3.32	−3.40		X	X
N-Acetyl-L-alanine	Amino acids and peptides	HMDB0000766	0	0	0	0	−0.07	0.02	−2.71	−2.69			X
NADH	Nicotinamides	HMDB0001487	0.86	1.66	−0.67	2.23	0.00	0.00	0.00	0.00										X

Metabolites were classified based on their biochemical structure and role in metabolic pathways. Metabolite abundance is presented in fold change against 3 h. Values in bold represent −2 > Log2Foldchange > 2 and underlined were statistically significant Padj < 0.05.

X, associated metabolic pathways.

To characterize the chemical diversity of the intracellular and secreted ETEC metabolomes and further explore the metabolic changes between the intracellular and secreted ETEC metabolomes across the time points, the significant metabolites were classified according to chemical classes (Human Metabolome DB; www.hmdb.ca) and a hierarchical clustering heatmap based on the metabolite relative abundances was generated. Ten and five metabolite classes were included in the intracellular and secreted metabolome, respectively. The most common class of metabolites were amino acids and peptides, as well as downstream catabolism metabolites. The metabolic profile changes illustrated in the heatmap of Fig. 3c indicated remarkable differences in the metabolite abundance within ETEC metabolomes and some strain-level differences in the metabolic profile. For instance, betaine and the fatty acid 5-aminopectanoic acid were only detected in E2265, and the amino acid N-methyl-l-glutamic acid only in E1777 (Fig. 3c).

Two main clusters of metabolite abundance patterns were identified: the first cluster included a diverse set of metabolites with higher intracellular concentrations or absence of secretion. For instance, S-adenosylmethionine, N-acetylcadaverine, GMP, NADH, glutathione, and pyridoxal 5′-phosphate (PLP, known as the catalytically active form of vitamin B6) were not secreted at any time point by the two strains.

In the second cluster with higher levels of secreted metabolites than intracellular levels, we found most of the amino acids identified in this data set and glucose and succinic acid. The amine spermidine was only detected in the E2265 supernatant, whereas the amino acid L-homoserine was secreted at higher levels by E2265 than E1777, where the metabolite was absent intracellularly (Fig. 3c). Again, these data indicated that although both ETEC strains are closely related genetically, their physiology could vary. The clustering patterns of samples per time point also confirmed significant differences in metabolites abundance over time.

The pathway enrichment analysis of significantly altered metabolites discovered general metabolic differences between both fractions. As is shown in Fig. 3d, the three KEGG pathways—aminoacyl-tRNA biosynthesis; butanoate metabolism; and alanine, aspartate, and glutamate metabolism—were differentially enriched (FDR < 0.05), and several other pathways were associated with the significant metabolites (Table 1). The intracellular metabolites with significant differences over the time points were enriched in the aminoacyl-tRNA biosynthesis pathway, including l-tyrosine, l-tryptophan, l-phenylalanine, l-leucine, and l-glutamic acid. These metabolites showed significantly decreased levels at 4 h which were restored in E1777 at 5 h but remained low in E2265. The intracellular glutathione levels increased 10-fold at 4 h compared with 3 h and remained high at 5 h. On the other hand, the pathways responsible for butanoate metabolism, alanine, aspartate, and glutamate metabolism were significantly enriched among secreted metabolites (Fig. 3c and Table S10).

The butanoate metabolism pathway included the bacterial neurotransmitter GABA, l-glutamate, succinic acid, and maleic acid, and they were mainly secreted at 3 h. These metabolites were found in lower concentrations than metabolites of the aminoacyl-tRNA biosynthesis pathway.

In order to investigate whether there were metabolite-metabolite correlations across the metabolomes, we calculated the Pearson correlation coefficients. We wanted to identify the strongest correlations, and a total of 110 pairwise differentially significant correlations (P < 0.001) were found, of which 74 were positive and 36 negative metabolite-metabolite correlations. Most of these correlations were consistent with the heatmap clustering; however, four sets of tightly correlated metabolite-metabolite interaction were identified. As seen in Fig. 3e (Table S11), all four grouped pairwise metabolites positively correlated (r < 0.8). Most of the amino acid metabolites were clustered in groups A and B, while the C and D groups included more diverse classes of metabolites. The positive correlations emphasized their close biochemical relatedness or overlapped roles among catabolic pathways. In contrast, metabolite-metabolite correlations from A and D as well as xanthine, thiamine, and uridine from the nucleotide metabolisms and A correlated negatively (r < −0.6).

Overall, the metabolome analysis of ETEC revealed distinctive differences in the metabolic composition of the intracellular and external environment of the bacteria, mainly characterized by the secretion or availability of large amounts of essential amino acids, intermediates, and derived molecules such as GABA. In contrast, the intracellular environment was characterized by a transient drop of amino acids at 4 h of growth, and glutathione redox-regulation was implicated in the adaptation to stationary phase as reduced L-glutathione was increased 10-fold intracellularly.

Transcriptome and metabolome integration.

Next, we combined the significant data from the transcriptomic analysis and metabolomics and mapped it into super pathways, including l-arginine, 4-aminobutanoate, and putrescine degradation (Fig. 4 and Table S12). The putrescine degradation pathway uses l-arginine as a carbon source and degrades putrescine to succinate for the TCA cycle. Specifically, at 4 h, the concentrations of l-arginine, GABA, the intermediate molecule in the polyamine putrescine degradation, and succinate had a significant drop by 5-fold in both intracellular and extracellular environments. They remained slightly lower toward the entry of the stationary phase. This was consistent with a significant transient upregulation of all puu operon genes at 4 h. The expression levels of puuR, the main repressor that regulates the intracellular putrescine concentrations by repressing several genes of putrescine utilization and transport, had the highest expression at 4 h. Glutamate and GABA are involved in the acid resistance system 2 (AR2), which enteric bacteria use to survive acidic conditions (19). A significant increase in the expression of gadA and the glutamate decarboxylase at 4 h suggests that bacteria initially convert l-glutamate to GABA, which is then subsequently degraded because levels of glutamate, GABA, and succinate all rapidly decrease intracellularly at 4 h compared with 3 h.

FIG 4

Transcriptomic and metabolic modulation of 4-butanoate degradation, related pathways, and Aminoacyl-tRNA biosynthesis of ETEC during growth. Circles are metabolites, and arrows are reactions. Verified metabolites are labeled in bold (Table S13). The color of the thick arrows indicates the levels of gene expression represented in the fold change. The first and second arrows represent the fold change between 3 h versus 4 h and 3 h versus 5 h, respectively. Concentrations of significant metabolites are represented in the serial charts. The metabolite concentration is adjusted to 3 h, and the data is presented in foldchange ±SD.

Putrescine degradation can also form other polyamines such as spermidine and spermine through l-methionine degradation. Spermine was found in large concentrations at 3 h and drastically decreased over time. The S-Adenosyl-l-methionine, synthesized from the essential amino acid l-methionine was found in low concentration at 4 h and 5 h (only in E1777) and subsequently converted by the spermidine synthase (encoded by speE) to spermidine and later spermine, which was found in lower concentration in the extracellular environment. Other acetylated polyamines derived from spermidine such as N1-acetylspermidine and N1, N8 diacetylspermidine were detected but not significantly changed over time. They displayed mixed patterns with increasing concentrations in both supernatant and intracellular fractions. Acetylation converts the polyamines to a physiologically inert form to protect against polyamine toxicity and is mainly excreted from the cell. On the other hand, spermidine can conjugate with thiol glutathione which plays a role in detoxifying xenobiotics and reactive oxygen species. This thiol was detected at increased concentrations over time intracellularly.

DISCUSSION

Growth in LB broth is commonly used to analyze bacterial properties. Although most natural habitats of prokaryotes do not resemble the nutrients in LB, the transition to more nutrient-depleted conditions, i.e., the stationary phase when growth ceases, is a common phenomenon in nature. We analyzed the transcriptome of two ETEC strains during the transition from late exponential phase/early stationary phase to stationary phase to get a comprehensive view of the regulatory and metabolic pathways involved in ETEC and E. coli growth in LB. Our results indicated several distinct steps during the transition into the stationary phase in support of previous studies (20). Specifically, we identified an interesting transient phase at OD600 = 3 at the immediate onset of stationary phase characterized by decreased gene expression of genes involved in iron uptake and more than 10-fold upregulation of operons involved in, e.g., dipeptide transport, fucose, and putrescine utilization and indole production. Furthermore, integration of the transcriptome with metabolome analyses highlighted the l-arginine, 4-aminobutanoate, and putrescine degradation pathways forcefully induced at the onset of the early stationary phase.

To our knowledge, few studies have focused on transient transcriptomic and/or metabolic responses across different growth phases. A similar study (21) of the E. coli O157:H7 (EHEC) transcriptome during growth phases in MOPS (morpholinepropanesulfonic acid) minimal medium also identified a later transient expression at the stationary phase transition point (4.6 h to 5 h and 5 h to 5.5 h). The transiently changed pathways were involved in nutrient scavenging, particularly alternative carbon sources, as a consequence of a glucose-limited condition and stringent response. These findings were also observed in other studies of E. coli K-12 growth in LB (22) and minimal media (23) as a short-term adaptation to a carbon-limited growth.

Other transient alterations in the transcriptome and/or metabolome of E. coli species have also been reported. Kim et al. (24) investigated the effect of heat (42°C) in E. coli strain BS21 and reported a heat-induced transient accumulation of alanine and adenosine-5´diphosphate during the first 10 min.

The exponential phase metabolomes (25) (OD600 ~ 1.0) of seven industrial strains of E. coli during aerobic growth on glucose did not report transient changes but instead revealed major differences in central metabolic pathways such as amino acids, TCA cycle, and redox. Similarly, transcriptional responses of EPEC prototype strains that grew in LB and DMEM media exhibited distinct transcriptional signatures across early, late exponential, and stationary phases (26). In summary, these findings suggest that although modulation of the transcriptomic and metabolic response can vary among closely related species and with carbon sources, the transition from active exponential division to completely ceased growth at stationary phase consists of several steps of different gene expression and metabolite fluxes. We identified an early transition phase at 4 h in this study, and the results after 5 h represent yet another phase that is different from the gene expression patterns found at the stationary phase.

Our results suggested that the transient phase at 4 h when ETEC/E. coli is preparing for stationary phase is characterized by a temporary reduction in iron uptake. Interestingly, a similar response has been reported in the transition phase between the log and stationary phase for Helicobacter pylori (27), suggesting this phase alteration and its characteristics occur in several bacterial species. Metal ions such as iron are essential for bacteria but highly toxic at the same time. Iron ions in oxidation state II (ferrous (Fe2+) iron are more bioavailable than ferric (Fe3+) ions in oxidation state III but more toxic since they may form hydroxyl radicals through the Fenton reaction. Our results demonstrated that, e.g., fecA and fhuA related to iron uptake and transport of ferrous ion and ferric ion were transiently downregulated during entry into stationary phase. This may have implications in ensuring iron homeostasis in an iron-deficient environment such as the gastrointestinal tract (28–30).

Transcriptome analysis also indicated increased fucose metabolism at the onset of stationary phase at 4 h followed by downregulated at 5 h. Fucose is an abundant mucus-derived metabolite in the intestine generated by commensal bacteria (31). The most common commensal bacterium in the gut is Bacteroides thetaiotaomicron, and some E. coli can utilize fucose as a carbon source (32). We confirmed that E2265 used in this study could metabolize L-fucose in a phenotypic assay, but interestingly this trait was not conserved over all ETEC lineages. Several studies have shown the important role of fucose in virulence. For example, in Salmonella Typhimurium, fucI was significantly upregulated 1 day after infection in germfree mice colonized by Bacteroides thetatiotaomicron, and the respective mutant had decreased competitiveness in vivo (33). Another example of fucose modulating virulence is EHEC, which harbors a pathogenicity island LEE containing a two-component system (TCS) capable of sensing fucose and transcribing the pathogenicity island LEE (34). In EHEC, it has also been demonstrated that fucose is important for colonization, suggesting that bacterial pathogens take advantage of using unexploited sugars by commensal bacteria (35). Hence, the induction of fucose utilization might promote colonization of the mucosa of bacteria in the early stationary growth phase. The fact that E2265 belongs to the commonly isolated and globally spread clonal lineage 5 (4) might implicate that fucose metabolism is important for virulence in certain ETEC lineages.

Pathway enrichment analysis identified two significant pathways, aminobutyrate (GABA) degradation, and 5-hydroxytryptamine (serotonin) degradation, to be significantly enriched in cluster III characterized by genes that progressively increased during entry into stationary phase (Fig. 2). Genes involved, e.g., puuE, gabD, aldAB, feaB, and prr increased their expression levels up to 10-fold at 4 h compared with 3 h and increased further in the last sampling point at 5 h.

The aminobutyrate degradation pathway (Fig. 4) is linked to l-arginine and putrescine degradation pathways. Putrescine and its downstream metabolite spermidine are polyamines present in the gut and introduced by food, microbial, or intestinal cell metabolism (36). They regulate cellular function in both prokaryotic and eukaryotic cells. The transient phase was characterized by upregulation of the puu-regulon. Putrescine is generated by decarboxylation of ornithine or decarboxylation of arginine into agmatine. E. coli uses specific importers such as PotFGHI and YdcU, members of the ATP-binding cassette (ABC) transporter family, and PuuP to take up putrescine across the cell membrane. Once putrescine is imported, PuuA γ-glutamylates putrescine resulting in a γ-glutamyl-γ-aminobutyraldehyde, which is oxidized by PuuB and subsequently dehydrogenated by PuuC. Then, PuuD hydrolyzes the resultant γ-glutamyl group, generating and releasing simultaneously γ-aminobutyrate (GABA) and glutamate. PuuE deaminates GABA to succinic semialdehyde, which is oxidated by GabD (which does not belong to the Puu operon) to produce succinic acid that goes to the TCA cycle (37–39).

In terms of virulence, polyamines, such as putrescine, have been shown to be actively involved. For instance, a mutant of potD in Streptococcus pneumoniae displayed attenuated virulence. In Vibrio cholerae, defective biofilm formation was observed after deletion of the homologous genes of E. coli potD. On the contrary, polyamines might prevent the colonization of the small intestine by V. cholerae because high concentrations of polyamines disrupt pili-pili interaction during autoaggregation (40). In Salmonella enterica serovar Typhimurium, exogenous putrescine and spermidine are sensed to prime intracellular survival and induce virulence (41). Dual-RNA seq of ETEC-V. cholerae cocultures revealed upregulated ETEC puuB and puuR genes as well as putrescine transport. In the same study, the authors used genome-scale reconstruction (GEM) approximations to predict that ETEC can grow on putrescine as sole nitrogen or carbon source while V. cholerae—a common pathogen found in co-infections with ETEC—does not grow on putrescine. These results indicate that during co-infection, ETEC can cross-feed to V. cholerae and support its growth (42). In light of existing literature, the increased capability of ETEC to grow and metabolize alternative nutrient sources in contrast to V. cholerae suggest that ETEC could exacerbate V. cholerae infection during co-infection in a polyamine-rich environment such as the human gut.

Interestingly, the downstream metabolite of putrescine, γ-aminobutyric acid (GABA), is a well-known neurotransmitter used as food supplementation in animal husbandry to reduce aggressive behavior and stress (43). Supplementation of GABA also induces sIgA secretion and increases IL-4 and IL-17 in piglets challenged with porcine ETEC (44). The AraC-like transcription factor GadX is an activator of the glutamate decarboxylase gadAB that converts l-glutamate to GABA. GadX is also a repressor of the transcription factor that activates the expression of adhesion factor bundle forming pilus (bfp) and intimin in enteropathogenic E. coli (45, 46). Hence, factors that increase GABA might also downregulate virulence factors, and the corresponding rapid degradation of GABA at the transient phase could promote virulence and downregulated immune responses to ETEC. We could, however, not see the corresponding pattern in this study because adhesion factors CS5 and CS6 were not changed. Previous studies on transcriptomes in ETEC have indicated that toxin and CF genes are downregulated upon binding to cells as well as influenced by bile salts present in the small intestine (11, 47, 48). We have also previously shown that eltAB expression decreases from exponential to early onset of stationary growth phase (8).

Tryptophan is cleaved to indole, pyruvate, and NH4+ by the tryptophanase (TnaA) enzyme expressed by certain bacteria, including E. coli, Bacteroides, and Lactobacilli. Indole is an extracellular signaling molecule well known for affecting different aspects of bacterial physiology, including biofilm formation, in a concentration-dependent manner (49–52). Indole is also an interspecies signaling molecule. In E. coli, the “indole peak” has been described in several studies as a short window of time at entry into the stationary phase where intracellular levels of indole rapidly peak and then decrease again (49). We found strong induction of tnaA expression at 4 h compared with 3 h and 5 h and concomitant rapid decrease of intracellular l-tryptophan at 4 h, supporting that the transient phase identified in this study is the same as the indole peak.

Tryptophan metabolism in the gut is important since tryptophan metabolites include indole and neurotransmitters and immunomodulators, including serotonin, tryptamine, and kynurenine. The synthesis of these latter molecules is performed in gut cells like enterochromaffin cells from diet-derived tryptophan and constitutes an important part of the gut-brain axis of neurotransmitters (53). In this study, metabolomics identified L-5-hydroxytryptophan (5-HTP), the precursor of serotonin (5-hydroxytryptamine), and we also found L-kynurenine (Table S6). To our knowledge, only one other study reported that E. coli could produce the neurotransmitter serotonin (54). We performed an additional verification analysis (Table S13) to search for serotonin in ETEC but were not able to detect it in our samples.

This study supports that commensal and pathogenic bacteria can both degrade and produce neurotransmitters and their intermediate molecules and hence ETEC infection might interact with gut cell signaling and influence the gut-brain axis to a larger extent than previously thought.

Besides the transient transcriptomic responses at 4 h, our time-series measurements also detected relevant, consistent increase or decrease in gene expression changes toward the early stationary phase that influenced a large portion of the ETEC transcriptome. In corroboration with our results, many downregulated genes (cluster I, Fig. 2) involved in nucleotide synthesis, translation, and transcription were also observed in E. coli K12 (23, 55) and EHEC O157:H7 (21) transcriptomes during the shift to stationary phase. These mechanisms of downregulation respond to a requirement to conserve energy, and it is also required to cope with stress response mechanisms such as the stringent response. Numerous studies indicated that genes involved in the TCA cycle are transcriptionally downregulated in stationary phase (23, 21, 51) while fully active at the early onset of the stationary phase (21). Considering these findings, our transcriptome-based observations show that the gene expression at 5 h in the present study reflects another phase just before entry into the stationary phase. Cluster III that enclose genes with gradually increased expression was significantly enriched with genes involved in the TCA cycle and aerobic respiration. In addition, we show that super pathways l-arginine, 4-aminobutanoate degradation, and putrescine degradation characterize this transition phase.

In summary, our data suggest that at the end of logarithmic growth, there is a distinct growth phase that might pose E. coli/ETEC into a stage of increased survival, virulence, and host competitiveness due to lack of need to sequester iron, retained virulence gene expression, and capacity to compete with the commensal flora for host-derived carbon and nitrogen sources such as fucose and putrescine. This study provides a framework for further studies on E. coli/ETEC gene regulation and comprehensive characterization of transcriptional responses during the transition to the stationary phase that also applies to other bacteria.

MATERIALS AND METHODS

Strains, growth conditions, and bacterial enumeration.

The ETEC strains E1777 and E2265 (LT STh/CS5+CS6), both isolated from adult patients with watery diarrhea in Dhaka Bangladesh in 2005 and 2006, respectively, were used in this study (12, 56). The whole-genome sequences of the two strains are available (18), including a complete assembled chromosome and two plasmids of 142 and 78 kbp, respectively, for E2265 (11, 56). Bacteria from frozen stock vials were grown on blood plates, and 10 colonies were picked and grown under shaking conditions in 10 mL of LB medium to OD600 = 0.8 (109 bacteria/mL) to be used as a starting culture. The starting culture was diluted 100-fold in 20 mL LB medium in a 250 mL Erlenmeyer flask and grown aerated at 150 rpm rotation at 37°C. Samples for optical density, CFU, and RNA extraction and metabolomics were withdrawn after 3 h, 4 h, and 5 h. The track dilution method was performed for bacterial enumeration as described previously (57). In brief, a 20 μL of bacterial culture was collected every hour up to 6 h and overnight time point and subjected to 10-fold serial dilution in a 96-well plate filled with 180 μL of phosphate-buffered saline (PBS) 1X. 10 μL from each dilution were spotted in a column onto LB agar plates, and the plate was tipped onto its side to allow migration of the spots across the agar surface. This step was performed in duplicate. LB plates were incubated overnight at 37°C, and CFU per mL was quantified by multiplying the number of colonies of each tract by their respective dilution factor and inoculated volume (0.01 mL).

For RNA extraction, bacterial samples from 3 h, 4 h, and 5 h time points were immediately mixed with 2 x volume of RNAProtect (Qiagen) using the manufacturer’s protocol. Samples were stored at −80°C until extraction.

RNA preparation.

Total RNA was prepared from lysozyme and proteinase K lysed bacteria using the RNeasy mini-kit (Qiagen, Hilden, Germany) and the instructions provided by the manufacturer for RNA extraction from Gram-negative bacteria. An extra step to remove contaminating DNA on-column was included using the RNase-Free DNase Set (Qiagen). The integrity of the RNA and absence of contaminating DNA was checked by agarose gel electrophoresis, and the RNA concentration was measured spectrophotometrically using a NanoDrop® ND-1000 (NanoDrop Technologies, Wilmington, DE). The RNA was carefully precipitated, washed, and shipped under 99.5% EtOH. The integrity, quality, and concentrations of the RNA were rechecked upon arrival at the sequencing facility at Beijing Genome Institute (BGI), Shenzhen, China, using an Agilent. The RIN values were above 9.8 for all samples.

RNA-seq.

One RNA sample per time point and isolate was depleted from rRNA by RiboZero, and Illumina libraries were generated using the TruSeq protocol described by the manufacturer. The libraries were sequenced using the Hi-Seq 200 using a read length of 100 bp. Reads were assembled using the software SOAPdenovo (65). CAP3 assembled all the unigenes from different samples to form a single set of non-redundant unigenes. All unigene sequences were blasted against protein databases using blastx (e-value < 0.00001) in the following order: Nr Swiss-Prot:KEGG:COG. Unigene sequences with hits in the first or second database did not go to the next search round against later databases. Then blast results were used to extract CDS from Unigene sequences and translate them into peptide sequences. BLAST results information was also used to train ESTScan (58). CDS of unigenes with no-hit in the blast were predicted by ESTScan and then translated into peptide sequences.

Functional annotations of Unigenes, including protein sequence similarity, KEGG pathway, COG, and GO, were performed. First, all-Unigene sequences were searched against protein databases (Nr Swiss-Prot KEGG COG) using blastx (e-value < 0.00001). Then, the Blast2GO program (59) was used to get GO annotations of the Unigenes. After getting GO annotation for every Unigene, we used WEGO software (60) to do GO functional classification for all Unigenes.

Each time point’s differential expression was determined by combining the E2265 and E1777 data sets as replicates (2 replicates per time point) and calculated in DESeq2 package (v1.22.1) and R-3.6.0 (61). The counts were normalized, and the fold change and the log2 of the fold change were calculated based on the following comparisons: control (3 h) versus 4 h, and control versus 5 h. Significant genes from each comparison and each strain were filtrated using the following threshold: Padj < 0.05, log2Foldchange < −2 or log2foldchange > 2. PCA and sample distance heatmaps were plotted to visualize the cluster of groups and outliers. Heatmaps of differential genes were generated using the R package pheatmap (62). For temporal gene expression pattern analysis, the Short time-series Expression miner (STEM) was applied as described by Ernst, Nau (17). For GO enrichment analysis for biological processes and metabolic pathways, PATHER Overrepresentation Test (http://www.pantherdb.org/) (63) was used with an FDR correction applied to all reported P values for the statistical tests.

Untargeted metabolomics.

(i) Bacterial sampling. The bacterial sampling for metabolomics was performed as described previously (64) with some modifications. A total of 4 mL of bacterial culture in triplicates were collected, and each culture was split into two sterile Eppendorf tubes to sample intracellular and secreted metabolites at every time point. In addition, another 100 μL were collected to measure the optical density and pH. For extracellular metabolites, 2 mL of bacterial culture was pelleted by centrifugation at 12,000 × g for 3 min a tabletop centrifuge, and the supernatants were carefully removed and transferred to a new sterile Eppendorf tube for snap-frizzing in liquid nitrogen. Snap frozen samples were stored at −80°C. The fast filtration method was applied for intracellular metabolites using a 3-place EZ-FitTM Manifold (Millipore) connected to a single vacuum that supports simultaneous filtration of three samples. A sterile 22-mm diameter MF-Millipore membrane filter with 0.45-μm pore size was placed onto each manifold, pre-washed with pre-warmed LB medium, and the vacuum set to 50 mbar. Then, 2 mL of the bacterial culture were pipetted in the middle of the filter and subsequently perfused with 5 mL of pre-warmed washing buffer (M9 medium [Sigma-Aldrich] adjusted to pH 7.3) was perfused. Immediately after, the cell-loaded filter was removed and transferred to an Eppendorf tube for snap-freezing in liquid nitrogen. Next, 2 mL aliquots of LB and M9 minimal medium (Sigma-Aldrich) were collected and span-frizzed and used as negative controls. All tubes were kept at −80°C and shipped in dry ice to the Science for Life Laboratory at Uppsala University for metabolite extraction and UPLC-MS analysis.

(ii) Metabolite extraction. Supernatants were extracted by adding 4 mL of 60:40 ethanol: water solution. The mixture was kept at 78°C for 3 min with vigorous mixing every minute. The filters containing the bacterial pellet were transferred to a 4 mL solution and kept at 78°C for 3 min, with vigorous mixing every minute. The samples were transferred to Eppendorf tubes on ice and centrifuged at 13,500 rpm for 5 min, at 4°C. The supernatant was collected and dried under vacuum on a Speedvac concentrator. The pellet was re-dissolved and injected onto the UPLC-MS system.

(iii) UPLC-MS Analysis. Ultra-high-performance liquid chromatography coupled to a mass spectrometer (UPLC-MS) was used to identify metabolites, which differ between the three different time points. Mass spectrometric analysis was performed on an Acquity UPLC system connected to a Synapt G2 Q-TOF mass spectrometer, both from Waters Corporation (Milford, MA, USA). The system was controlled using the MassLynx software package v 4.1, also from Waters. The separation was performed on an Acquity UPLC HSS T3 column (1.8 μm, 100 × 2.1 mm) from Waters Corporation. The mobile phase consisted of 0.1% formic acid in MilliQ water (A) and 0.1% formic acid in LC-MS grade methanol (B). The column temperature was 40°C and the mobile phase gradient applied was as follows: 0 to 2 min, 0% B; 2 to 15 min, 0% to 100% B; 15 to 18 min, 100% B; 18 to 20 min, 100% to 0% B; 20 to 25 min, 0% B, with a flow rate of 0.3 mL/min.

The samples were introduced into the q-TOF using positive electrospray ionization. The capillary voltage was set to 2.50 kV and the cone voltage was 40 V. The source temperature was 100°C, the cone gas flow 50 l/min, and the desolvation gas flow 600 l/h. The instrument was operated in MSE mode, the scan range was m/z = 50 to 1,200, and the scan time was 0.3 s. A solution of sodium formate (0.5 mM in 2-propanol: water, 90:10, vol/vol) was used to calibrate the instrument, and a solution of leucine-encephalin (2 ng/μL in acetonitrile: 0.1% formic acid in the water, 50:50, vol/vol) was used for the lock mass correction at an injection rate of 30 s.

Data analysis.

The obtained UPLC-MS data comparing the different time points were analyzed using the XCMS software package under R (version 3.3.0) to perform peak detection, alignment, peak filling, and integration. The peaks were annotated by comparing their m/z values to the exact molecular masses of all the online platform MetaboAnalyst for pathway analysis. The E. coli K-12 MG1655 from the KEGG database was used for metabolite identification. Confirmed metabolites were co-injected with the bacterial samples for the highest level of confirmation. The structures for the significantly altered metabolites were validated with authentic internal standards, as detailed in Fig. S12. PCA and heatmaps were performed in R software. The metabolite-metabolite correlations were computed using the R function rcorr from the Hmisc package (https://cran.r-project.org/web/packages/Hmisc/index.html). A multiple parametric statistic t test was performed in GraphPad to compare the means of two paired groups, i.e., 3 h versus 4 h, and multiple comparison corrections were applied using the Holm-Šídák method. P-value < 0.05 was set as the threshold for significance.

Data availability.

The metabolomic data used in this study has been submitted to MetaboLights under the accession numbers MTBLS3540. The RNA-Seq raw data are unavailable, but Table S1 contains the raw read counts per gene and their functional annotation.