2 Introduction

The Long Term Evolution Experiment (LTEE) is an experimental evolution system where 12 initially identical but now independent populations of E. coli are propagated in a carbon-limited medium with a 24-hour serial transfer regime [1]. The experiment was started in 1988 and has thus far accumulated over 75,000 generations of growth, making it the longest running experiment of its kind. Over time, all 12 populations have continued to adapt to this environment. For example, relative fitness (measured as growth rate relative to the ancestral strain) has increased in parallel across the 12 populations [2]. Cellular size in each population has also increased [3, 4]. Parallel genomic changes in the LTEE are well characterized, and while some genes are commonly mutated across the replicate populations, most mutations are unique to each line [5–8]. Despite the variability at the genetic level, the evolved populations display similar gene expression profiles [9, 10]. In principle, mutations and changes to gene expression could affect the proteome by altering the function of a protein, the amount of protein produced, or the conditions under which it is translated. Because proteins are the biological catalysts for cellular chemistry, mutations and gene expression changes may affect the metabolome. However, aside from select examples, the relationship between mutations, expression changes, metabolic changes, and how these may impact fitness in the LTEE is unknown.

Metabolic changes have been implicated in adaptive processes in a variety of organisms. An analysis of metabolic differences between species of Drosophila found that these differences may be linked to lifespan and sex [11]. Many organisms produce or engage in symbiotic relationships with organisms that produce biological toxins through specialized metabolic pathways. Examples include the production of formic acid by certain insects [12] and the many marine organisms, such as the Blue-ringed Octopus (genus Hapalochlaena), that house bacterial symbionts that produce tetrodotoxin [13]. In prokaryotes, studies suggest that the genes involved in secondary metabolite production evolve quickly and produce large varieties of phenotypes [14, 15].

In the LTEE, several examples of metabolic changes have been documented. The Ara-3 population has evolved the ability to metabolize citrate under aerobic conditions [16, 17] and Ara-2 has produced unique growth-phase specific ecotypes [18–20]. These ecotypes exhibit growth phase-dependent frequencies and differ in their ability to metabolize acetate, with one ecotype (the L ecotype) growing while glucose is abundant during the exponential phase and the other (the S ecotype) growing during the stationary phase when acetate is abundant. The transposable element induced loss of the rbs operon occurred early and prevents the evolved lines from growing solely on ribose [21]. Reduced ability to grow on maltose has occurred through other mechanisms [22–24]. Despite glucose being the only source of carbon provided, significant changes in the ability to grow on other carbon sources have occurred [25]. In particular, the mutator lines exhibit reduced performance on carbon sources other than glucose. The evolved lines have increased rates of glucose uptake as well as altered flux through some central carbon metabolism pathways [26].

While specific examples of how genetic changes result in metabolic changes in the LTEE (citrate, acetate, maltose, ribose) exist, a more general survey of the metabolome is lacking. Such a survey would allow a better understanding of the effects of genetic and expression changes on the metabolome and provide the ability to test hypotheses from previous data. By integrating genomic [5], expression [10], and metabolic datasets, we can relate mutations to expression changes and their effects at the metabolic level. We provide ample evidence for hypotheses that specific genetic changes may exert their fitness-altering effects by changing gene expression and ultimately, the metabolome.

3 Results

3.1 Survey of metabolic changes in the LTEE

We used liquid chromatography coupled mass-spectrometry (LC/MS) to scan a broad mass range using both positive and negative ionization modes to survey the metabolomes of both ancestral strains and clones from each of the 12 evolved lines of the LTEE at 50,000 generations. We sampled two time points, 2 and 24 hours, to represent exponential and stationary phase samples, each with two biological replicates. Metabolite abundances were estimated using normalized peak areas relative to standards as described in the Data processing and description. Many metabolites are detected in both ionization modes, making it unclear which mode more accurately represents the abundance of the compound in a sample. As a result, we used the combination of a compound and the ionization mode it was detected in as a feature of the data.

We begin by surveying the data in a metabolome-wide manner. Because E. coli exhibits large physiological differences across growth phases [27–29], we expected intraphase comparisons of samples to be more similar than interphase comparisons. Indeed, the distribution of intraphase correlations of metabolite abundances is higher than the distribution of interphase correlations (Fig. S1A). Because growth-specific differences have the potential to drown out signals of adaptation, we performed all subsequent analyses in a growth phase-specific manner.

To identify large scale differences in metabolite abundances, we performed a principal component analysis (PCA) of metabolite abundances in a phase-specific manner. In the exponential phase, the PCA strongly separates Ara-3 from the other samples (Fig. 1A). This is likely due to the unique ability of Ara-3 to metabolize citrate. Citrate is present in the medium originally as a chelating agent. However, mutations in Ara-3 have allowed it to metabolize citrate under aerobic conditions, giving it access to additional carbon that the other lines do not have [30]. Half of the evolved lines have a mutator phenotype and have at least an order of magnitude higher number of fixed point mutations than the non-mutators [5]. Despite the large differences in mutation accumulations, we do not observe a strong effect of mutator status on the separation between evolved lines in the PCA. This is consistent with previous observations of a minimal effect of mutator status on changes in gene-expression patterns [10]. The first two principal components are dominated by contributions of nucleoside monophosphates, amino acids, and compounds involved in carbon metabolism, such as glucose and succinate (Fig. S3A-B).

Figure 1.

Similar to the exponential phase, Ara-3 separates from the rest of the evolved lines in a PCA of metabolites in the stationary phase (Fig. 1B). The primary driver of variation in the first principal component (PC1) is the relative abundances of compounds in Ara-3 (Fig. S4A). In addition to Ara-3, Ara-2 also shows a higher degree of separation from the other lines. Though this experiment was performed with single clones, the 12 flasks in the LTEE are not isoclonal, but instead consist of competing subclones that sometimes coexist for extended periods [7]. Ara-2 is the strongest example of this, consisting of two major subclones called L and S [18–20]. The observed separation might be because the Ara-2 clone in this experiment comes from the L ecotype, which is more strongly adapted to the exponential phase and has poorer growth relative to it’s S ecotype counterpart in the stationary phase. In contrast to exponential phase PCA, mutator status appears to be a driving factor of PC2 in the stationary phase (Fig. S4A-B) with mutator lines having lower abundances of nucleoside monophosphates relative to non-mutators during stationary phase.

A common theme in studies of LTEE is the high degree of parallelism observed at the genetic [5], gene expression [10], and fitness [31]. For example, genes involved in flagellar biosynthesis and amino acid metabolism are commonly mutated across the evolved lines [5] with associated changes in expression levels [10]. To examine the extent of parallel changes in metabolites, we calculated the ratio of the peak area in an evolved line to the peak area in the ancestor for each compound after averaging biological replicates. Pairwise Spearman’s correlations using these ratios are significantly higher than expected when comparing the observed distribution to a distribution generated by randomizing fold-changes within a sample prior to calculating correlations (Fig. 1C, p < 0.001, t-test). In the exponential phase, the lowest correlations were mainly between a pair involving one of Ara-2, Ara-4, or Ara-6, and another evolved line (Fig. 1D). The correlations of metabolite abundances across evolved lines are more similar in the exponential phase (Fig. 1D) than in the stationary phase (Fig. 1E) (t-test comparing the observed exponential and stationary phase correlations in Fig. 1C, p < 0.001). It may be that the strategies of metabolite usage and abundance in the rapidly growing exponential phase are more similar across lines than the strategies of survival are in the stationary phase.

To further quantify the extent of shared changes, we considered the difference between the number of shared changes that were observed and expected to be observed by chance. We approximated an expected distribution of shared changes using a Sum of Independent Non-Identical Binomial Random Variables (SINIB) method [32], essentially asking what the chance of repeatedly observing the same change in an increasing number of evolved lines is. For more details on this method, see section 5.1.5. For metabolic features (the combination of a compound and the ionization mode it was detected in) that shared alterations in only a few lines (generally six or less), we cannot suggest selection on these changes because we see fewer or as many changes as expected (Fig. S5, Fig. S6). These may represent metabolites that negatively affect fitness when altered. On the other hand, we observe more changes shared in a higher number of evolved lines than expected. These shared changes may be beneficial, but this would need to be clarified by additional experiments. By combining genomic, expression, and metabolic datasets, we can better explore the nature of these changes and how they may impact fitness.

3.2 Increased NAD may facilitate higher energy demands during adap-tation

In the LTEE, mutations to the nadR gene have occurred in all of the evolved lines and appear to have been under strong selection to improve fitness [5, 33]. NadR is a dual-function protein that negatively regulates genes involved in the synthesis of nicotinamide adenine dinucleotide (NAD) and has its own kinase activity [34, 35].

NAD and its phosphorylated derivative NADP function as electron carriers and are necessary for many catabolic and anabolic pathways in the cell. While untested, hypotheses for how nadR mutations might increase fitness involve derepression of its target genes and increases to intracellular NAD; this may aid in shortening the lag phase and modulating redox chemistry [33, 37, 38]. The abundance of NAD and related compounds in a cell can modulate reaction rates and other processes [35, 39]. We had previously shown that mutations in nadR in the evolved lines lead to consistent upregulation of its target genes [10] (Fig. 2). Here we test if there are indeed higher abundances of NAD and NAD-related metabolites in the evolved lines consistent with changes observed at the genetic and gene expression levels.

Figure 2.

The genes regulated by NadR participate in specific reactions along NAD synthesis pathways. We detected some of the compounds in these pathways as well as both oxidized and reduced forms of NAD and NADP (Fig. 2). Compared to the ancestors, both redox states of NAD and NADP are almost universally increased in the evolved lines during the exponential phase (Fig. 2, Fig. S7), with median fold-changes of 3.32 and 4.65 for NAD and NADH, respectively. NADP is generated through phosphorylation of NAD by NadK [35, 40] and saw median fold-changes of 2.34 and 3.54 for NADP and NADPH, respectively. Because NADP is generated from NAD, and each exists in two oxidation states, we expected that all four compounds would see similar increases within an evolved line. Indeed, increases in the various redox and phosphorylation states were consistent within an evolved line (0.77 < R < 0.97, Fig. S7B). Aspartate, the starting point for NAD synthesis, was also increased in most of the evolved lines, with a median fold-change of 1.94. Nicotinamide mononucleotide, which can be used to make NAD, was consistently increased, with a median fold-change of 2.51. In the stationary phase, most of these patterns remain unchanged except for aspartate, which is generally lower in abundance compared to the ancestor (Fig. S7A, S3C). These results suggest that the end product of the mutations in nadR is increased intracellular NAD.

3.3 Derepression of Arginine biosynthesis at genetic, gene expression, and metabolic levels

Similar to nadR, another consistent target of mutations in the LTEE is a transcriptional repressor of genes involved in arginine biosynthesis argR [5]. ArgR, along with arginine, negatively autoregulates the arginine biosynthesis pathway by downregulating the participating genes when arginine is abundant [41]. As an amino acid, the primary role of arginine is in protein synthesis. However, various reactions produce and consume arginine; both its synthesis and the regulation of the enzymes involved are complex (see [42] and the EcoCyc [36] pathway “superpathway of arginine and polyamine biosynthesis”). For example, arginine can be a source of nitrogen and carbon or be used for ATP synthesis [43]. Because of these complications, if mutations in argR and the subsequent changes in gene expression do affect this pathway, it is unclear at which step one might observe a change or if that change might be in a different pathway.

Mutations in argR are localized to specific regions, potentially leading to disruption of the DNA binding domain, the arginine sensing domain, or interaction domain, which is required for hexermerization of the protein [41]. These mutations in argR could reduce its ability to repress its target genes, resulting in their increased expression. Consistent with these expectations, we had previously observed parallel increases in the expression of genes repressed by ArgR (Fig. 3) [10].

Figure 3.

We were able to detect 11 compounds involved in the superpathway of arginine and polyamine biosynthesis (Fig. 3, Fig. S8A). Compared to the ancestor, arginine was increased in nine (negative mode) or ten (positive mode) evolved lines and experienced a median fold-change of 1.9 across both ionization modes in the exponential phase. This pattern also persisted in the stationary phase (Fig. S8A). While arginine is increased in both growth phases in most of the evolved lines, other compounds in this pathway show highly variable changes (Fig. S8A). Despite the consistent relationship between mutations, expression changes, and changes in the abundance of arginine, how this might affect fitness is not obvious. Higher intracellular abundances of amino acids could facilitate higher translation rates and promote faster growth. Consistent with this hypothesis, we observe that most amino acids show increases in their abundances in most evolved lines in the exponential phase (Fig. S8).

3.4 Functional changes in the central carbon metabolism

Low carbon availability is the key feature of the minimal medium used in the LTEE, and selection for better utilization of carbon is a driving force of adaptation in the system. Hence, how mutations and expression changes might affect fitness by altering central carbon metabolism pathways is of interest. The evolved lines have seen significant positive and negative changes in their ability to grow on different carbon substrates. Rather than metabolic specialization, this was found to be due to mutation accumulation [25]. Moreover, low glucose has driven innovations such as citrate metabolism in Ara-3, granting it access to extra carbon from an unused ecological niche [16, 17]. The evolved lines retain more carbon in their biomass compared to the ancestors [30], and increased glucose uptake rates have been demonstrated [26]. Despite the abundance of molecular data associated with the LTEE, the highly networked nature of the genes and metabolites involved in central carbon metabolism makes relating the various molecular data to each other challenging. In particular, mass-spectrometry data alone will not allow us to differentiate between sources for compounds that are generated by many reactions. We might overcome this limitation by combining genomic, expression, and metabolic datasets.

We again sought instances of parallel sets of mutations and expression changes that may exert their effects at the metabolic level. Three genes that are key components in the glyoxylate cycle, aceB, aceA, and aceK were commonly upregulated in the evolved lines [10]. Their transcriptional regulators icIR and arcB are heavily mutated [5], and this is known to be beneficial in the LTEE [17]. The glyoxylate cycle allows the use of acetate as a carbon source by converting it to succinate, which can later be fed into the citric acid cycle [36, 44]. Because these genes had increased expression, the compounds they produce (succinate, glyoxylate, and malate) may be elevated in the evolved lines relative to the ancestor. Succinate and malate were reliably identified in our data, but glyoxylate was not. Unfortunately, acetate was part of the mobile phase of our liquid chromatography setup, thus preventing its accurate quantification. Likewise, most of the molecules in the glyoxylate cycle, including succinate and malate, are also components of the citric acid cycle and can be produced or consumed by other enzymes. Because we cannot distinguish where changes in these metabolites come from, we chose to look at compounds from central carbon metabolism in general, considering compounds from the glyoxylate cycle, glycolysis, and the citric acid cycle.

Overall, 18 compounds from these pathways are detected in our LC/MS assay. Regarding the glyoxylate shunt, malate, aconitate, and succinate were generally elevated in the evolved lines (median fold-change of 5.38, 1.66, and 2.30 respectively (Fig. 4)). Interestingly, the evolved lines appear to have similar or lower amounts of glucose despite having been shown to have increased glucose uptake rates [26]. This discrepancy is likely due to the fact that in [26], glucose uptake rates were calculated by measuring the depletion of glucose in the medium, whereas we measured glucose from inside the cells. Once inside the cell, we suspect that the glucose is quickly utilized. This is supported by the fact that other downstream glycolysis compounds, like phosphoenolpyruvate, are generally elevated (median fold-change of 6.72, Fig. 4).

Figure 4.

The citric acid cycle is a major metabolic pathway that extracts energy via the reduction of electron carriers like NAD. These electrons are shuttled to the electron transport chain by NAD, where they are used to generate ATP. NAD can be reduced at two points during the citric acid cycle, the conversion of alpha-ketoglutarate to succinate and the conversion of malate to oxaloacetate. If the citric acid cycle is running faster in the evolved lines due to their higher growth rates and hence higher energy demands, then more NAD may be required to efficiently shuttle electrons from the citric acid cycle to the electron transport chain. Interestingly, we measured a median fold-change of 3.13 for alpha-ketoglutarate and 5.38 for malate Fig. 4. Additionally, increases in these compounds are correlated with an increase in NAD compounds within an evolved line (Fig. S9). It may be that increases in NAD allow faster operation of the citric acid cycle. However, experiments that utilize radioactive tracers to study flux through these pathways are necessary to confirm this hypothesis.

4 Discussion

The metabolome of a cell is the integration point of an organism’s environment, genetics, and gene expression patterns. Metabolism has been shown to contribute to generating adaptive phenotypes [11–14]. However, how mutations affect metabolism through changes in gene expression and ultimately affect fitness is less clear. Specific examples of this have been studied in the LTEE, such as the acquisition of citrate metabolism in Ara-3 [16, 17]. We used metabolic mass-spectrometry to study other aspects of the LTEE, relating mutations and expression changes to metabolic changes. This allowed us to lend support to previous hypotheses of how certain mutations affect fitness in the system. For instance, we show how mutations in nadR alter the expression of its target genes and increase the intracellular abundance of NAD. It also enables new areas of investigation, such as what role mutations to argR and subsequent expression and metabolic changes might mean for the system.

A key limitation of the data presented here is that it contains measurements of 196 metabolites out of the 3755 metabolites currently annotated in the E. coli metabolome database [45]. Nonetheless, analysis of this subset of data can reveal some global patterns of metabolite changes. For example, distinct clustering of Ara-3 and Ara-2 in the PCA’s suggests that even with a limited sampling of the metabolome, the unique characteristics of these lineages are observable at the metabolic level. Interestingly, while the mutator lines do not form a well-defined group in the exponential phase, their reduced abundances of nucleoside monophosphates in stationary phase groups them together (Fig. 1B). Due to deletions of key biosynthetic genes, the mutators are less flexible than the non-mutators in their ability to grow on different carbon sources [25]. Ara-6, in particular, has one of the higher mutational burdens and is the least flexible in its ability to grow on different carbon substrates. Pathways with missing or broken genes may cause the shunting of those compounds to different pathways while also starving downstream reactions of the now missing reactants.

The analysis of the evolved metabolomes in the LTEE presented here has focused on specific pathways. However, due to the networked nature of the metabolome, the exact consequences of any particular perturbation are hard to predict. Experiments in other systems have shown that even altering the dosage of a single gene can have cascading and convoluted effects [46]. In LTEE, this phenomenon is further compounded by the fact that the evolved genomes are highly variable and typically lack tens to hundreds of genes. The work presented here lays the groundwork for comprehensive whole-cell modeling of expression and metabolism [47–49] to tease apart the relative contributions of individual gene expression and metabolic changes to fitness. Our strategy of relating genetic, expression, and metabolic changes bring us one step closer to understanding the mechanisms by which mutations lead to fitness increases in the LTEE. Our study shows the benefit of combining multiple datasets to study the impact of mutations in an evolutionary context.

Acknowledgements

We thank Richard Lenski for generously providing clones from the LTEE. Premal Shah is supported by NIH/NIGMS grant R35 GM124976 and start-up funds from the Human Genetics Institute of New Jersey at Rutgers University. Samhita S. Yadavalli is supported by NIGMS R35 GM147566 and institutional start-up funds. Mass-spectrometry data were generated by the Rutgers Cancer Institute of New Jersey Metabolomics Shared Resource, supported in part with funding from NCI-CCSG P30CA072720-5923.

Competing interests

Premal Shah is a member of the Scientific Advisory Board of Trestle Biosciences and is a Director at Ananke Therapeutics. Srujana S Yadavalli consults and collaborates with Designs for Vision Inc. Kyle S. Skalenko is a scientist at Specialty Assays Inc.

5 Supplementary text

5.1 Methods

5.1.1 Cell culture

The Ara-3 clone used in this experiment is capable of metabolizing citrate, and the clone of Ara-2 is the L variant [10]. Frozen bacterial stocks were revived by growing them in 5mL LB broth for 24 hours. Following this, 1% of each culture was transferred to 5mL standard LTEE medium for another 24 hours. After 24 hours, these cultures were used to initiate the final experimental cultures. Each clone was grown in 250mL of medium in a 1L flask following standard LTEE protocols. After 2 hours of growth, 160mL of the sample was removed for exponential phase mass-spectrometry samples. After 24 hours of growth, 40mL of the sample was removed for stationary phase samples. Each line had two independent biological replicates.

5.1.2 Mass-spectrometry sample collection

Cells were collected via vacuum filtration using Millipore Omnipore 0.2 um PTFE filters (JGWP04700 Sterile plastic petri dishes were placed onto a metal tray on dry ice, and 1.2mL of extraction solvent (40:40:20 acetonitrile: methanol: water + 0.5% formic acid) was added to the petri dish to chill. When filtration was complete, the filter was placed cell-side down into the petri dish. The metal tray was moved to wet ice to extract for 20 minutes. After 20 minutes, the acid was neutralized by the addition of 1.07% (final) ammonium bicarbonate, then cell debris was spun down at 14,000 rpm at 4C for 10 minutes, and the clarified extract was transferred to a chilled 2mL tube and stored at −20C. The filter was extracted again by the addition of 0.4mL extraction solvent and sat for 15 minutes. The extract was neutralized, clarified, and consolidated with the first extract and incubated overnight at −20C.

The next day, precipitated protein was spun down at 14,000 rpm at 4C for 10 minutes, and solvents were removed by speed vac for 2 hours at room temperature. Lastly, the sample was concentrated by removing the remaining water by lyophilization. Dried samples were reconstituted in 40uL of extraction solvent and submitted for mass spectrometry.

5.1.3 Mass-spectrometry

Mass spectrometry was performed as described previously in [50] at the Rutgers University Cancer Institute of New Jersey metabolomics core facility using a Thermo Q Exactive PLUS coupled with a Thermo Vanquish UHPLC system.

5.1.4 Data processing and description

The raw mass spectrometry data is deposited at the metabolomics workbench under the study ID ST002431 and the code for the analysis at https://github.com/shahlab/ltee_massspec. Normalization of the mass-spectrometry data was performed in two steps. First, raw peak areas were normalized against internal standards as in [50]. After these values were generated, all other data processing steps were performed using the R programming language [51] and the tidyverse set of packages [52]. Code for this section can be found in the document titled data_processing.Rmd. Peak areas resulting from the first step of normalization were then taken as a proportion of the total peak area for a single sample to normalize against differences in input amounts.

As noted in the text, not all of the 196 compounds were detected in all samples. This varied depending on the compound, line, growth phase, and ionization mode. Overall, 168 compounds were present in all samples in at least one ionization mode, with the remaining 28 compounds being undetected in at least one of the samples in at least one of the ionization modes. We used methods from [53] and the accompanying R package imputeLCMD to evaluate different imputation methods for imputing missing data, settling on the quantile regression imputation of left-censored (QRILC) data method. Imputed values theoretically represent values below the limit of detection rather than a complete absence of compounds. As expected, our imputed values always fell below the detected values for a given compound (Fig. S2). After imputation was performed, correlations between the replicates were high (Fig. S1B-C). Compounds that were detected in both ionization modes show a modest correlation within a sample (Fig. 1D). Distributions of normalized peak areas are similar across replicates and samples (Fig. S1C). This completed dataset was used for further analysis and is available as supplementary table 1. Where appropriate, we show data from both ionization modes, treating the combination of a compound and the ionization mode it was detected in as a feature of the data and compare these across samples.

5.1.5 Theoretical distributions for parallel changes in metabolites

The code for this analysis can be found in the document titled parallelism.Rmd. We used the R package SINIB [32] to calculate theoretical probabilities of finding a shared metabolic change in a particular number of evolved lines. First, we designated metabolic features (the combination of compound and the ionization mode it was detected in) as significant if they experienced an |log₂(fold–change)| ≥ 1. Then, we determined the probability of randomly picking a significantly altered metabolic feature in each evolved line in a manner specific to the growth phase and direction of change. We then parameterized the dsinib function with these probabilities, essentially asking what the chance of repeatedly finding the same change in an increasing number of evolved lines is. We used these probabilities to determine the total number of metabolic features ones might expect to find altered in the same direction. These numbers can be compared to the observed distribution, showing that more parallel changes are observed than expected.

6 Supplementary figures

Supplementary Figure 1.

Supplementary Figure 2.

Supplementary Figure 3.

Supplementary Figure 4.

Supplementary Figure 5.

Supplementary Figure 6.

Supplementary Figure 7.

Supplementary Figure 8.

Supplementary Figure 9.