Bacterial populations have a remarkable capacity to cope with extreme environmental fluctuations in their natural environments. In certain cases, adaptation to one stressful environment provides a fitness advantage when cells are exposed to a second stressor, a phenomenon that has been coined as cross‐stress protection. A tantalizing question in bacterial physiology is how the cross‐stress behavior emerges during evolutionary adaptation and what the genetic basis of acquired stress resistance is. To address these questions, we evolved Escherichia coli cells over 500 generations in five environments that include four abiotic stressors. Through growth profiling and competition assays, we identified several cases of positive and negative cross‐stress behavior that span all strain–stress combinations. Resequencing the genomes of the evolved strains resulted in the identification of several mutations and gene amplifications, whose fitness effect was further assessed by mutation reversal and competition assays. Transcriptional profiling of all strains under a specific stress, NaCl‐induced osmotic stress, and integration with resequencing data further elucidated the regulatory responses and genes that are involved in this phenomenon. Our results suggest that cross‐stress dependencies are ubiquitous, highly interconnected, and can emerge within short timeframes. The high adaptive potential that we observed argues that bacterial populations occupy a genotypic space that enables a high phenotypic plasticity during adaptation in fluctuating environments.
Escherichia coli cells were evolved over 500 generations and profiled in four abiotic stressors to observe several cases of emerging cross‐stress behavior whereby adaptation to one stressful environment provided fitness advantage when exposed to a second stressor.
Cross‐stress dependencies were found to be ubiquitous, highly interconnected and can emerge within short timeframes.
Several targets were implicated in adaptation and cross‐stress protection, including genes related to iron transport and flagella.
Adaptation in a first stress can lead to higher fitness to a second stress when compared with cells adapted only in the latter environment.
Adaptation to any specific stress and the growth media was found to be generally independent.
Outside the laboratory, microbes rarely live in conditions of optimum growth as environmental parameters constantly fluctuate. ‘Stress’ is the term used when these fluctuations vary severely from the optimum, and the ability of microbes to sense and respond to environmental stress is critical to their survival. Due to their importance to medical, agricultural and industrial fields, the underlying response mechanisms have been studied extensively and as a result we now understand their key components (Storz and Hengge‐Aronis, 2000).
Previous work on the adaptive potential of the bacterium Escherichia coli to fluctuations in temperature (Riehle et al, 2003; Bennett and Lenski, 2007; Sleight et al, 2008), UV irradiation (Alcántara‐Díaz et al, 2004), high ethanol (Goodarzi et al, 2010), isobutanol (Atsumi et al, 2008), and various carbon sources (Silver and Mateles, 1969; Dekel and Alon, 2005; Blount et al, 2008) has provided initial evidence regarding the genetic basis of microbial long‐term adaptation to abiotic environmental factors. However, little is known about how stress selection affects the capacity of an organism to be protected under other stresses. ‘Cross‐stress’ protection, that is, acquired stress resistance to a specific stress after pretreatment with another stressor, has been documented for various stress combinations and it has been observed in many species across the tree of life, ranging from microbes (Tesone et al, 1981; Jenkins et al, 1988, 1990; Leyer and Johnson, 1993; Mattimore and Battista, 1996; Lou and Yousef, 1997; Fletcher and Csonka, 1998; Garren et al, 1998; Ryu and Beuchat, 1998; Cánovas et al, 2001; Begley et al, 2002; Sghaier et al, 2007; Berry and Gasch, 2008; Machielsen et al, 2010), to plants (Schenk et al, 2000; Durrant and Dong, 2004) and humans (Raffaghello et al, 2008; Zhao, 2009). For example in E. coli, pre‐adaptation to glucose or nitrogen limitation increases survival rates after a heat shock or oxidative stress (Jenkins et al, 1988), and a possible link was suggested to be RpoH, a heat shock regulator, which was found to have an important role in protein synthesis under carbon starvation (Jenkins et al, 1991). Similarly, increased survival was observed after osmotic shock, when cells were pre‐treated with a mild increase in osmolarity (Jenkins et al, 1990). Transcriptional profiling in E. coli revealed a high degree of overlap between different stresses (Supplementary Results; Supplementary Figure S1), such as starvation, osmotic and acidic stress (∼140 genes) (Weber et al, 2005), and it was later found that high osmolarity or high temperature induces the oxidative stress regulons (SoxRS and OxyR) that can at least partially explain the cross‐stress protection between these stresses (Gunasekera et al, 2008; White‐Ziegler et al, 2008). Responses to n‐butanol were recently found to share the same high overlap with those in heat shock, oxidative and acidic stress (Rutherford et al, 2010). Interestingly, recent work with gene‐deletion libraries in yeast shows that acquired hydrogen peroxide (H2O2) tolerance after three distinct pretreatments occurs by different mechanisms depending on prior cellular experiences (Berry et al, 2011).
Given the importance of cross‐stress protection in bacterial physiology, we investigate the extent to which cross‐stress protection emerges spontaneously during microbial evolution, the evolutionary trade‐offs in stress selection, and the genetic basis of both acquired resistance and evolutionary cross‐stress dependencies. Towards this goal, we performed laboratory evolution of E. coli cells over 500 generations in minimal media and four stresses: osmotic, acidic, oxidative, and n‐butanol (Figure 1). Whole genome resequencing (Araya et al, 2010; Cooper and Lenski, 2010; Conrad et al, 2011) and transcriptional profiling (RNA‐Seq) revealed several genetic and transcriptional differences between the control and evolved strains. To understand the genotypic–phenotypic dependencies and genetic basis of acquired stress resistance during laboratory selection, mutations were individually reversed and the fitness of the corresponding repaired mutants was assessed. Systematic analysis of growth curves and competition assays among the evolved strains reveals a wide spectrum of evolutionary cross‐stress protection and evolutionary trade‐off cases for various strain–environment combinations.
Evolutionary potential and cross‐stress behavior
We adapted E. coli strains for 500 generations in five environments (osmotic, acidic, oxidative, n‐butanol, and control; four biological replicates per environment) with minimal M9 salts media and glucose as the sole carbon source (Figure 1). Prior to evolutionary adaptation in these environments we tested the MG1655 and ΔlacZ strains that were used in order to enable direct competition assays of evolved populations, in all growth environments to ensure neutrality of the lacZ deletion marker in all conditions (Supplementary Results; Supplementary Table S‐I; Supplementary Figures S2–S4). Growth curves and maximum growth rates of all strain–environment combinations were initially acquired to assess the growth differences in the respective populations (Supplementary Table S‐II and S‐III; Supplementary Figures S7–S9). Competition assays between adapted strains in all environments were then performed at a population level to obtain a relative fitness profile for each strain (four biological replicates per environment, competition over population averages, Supplementary Tables S‐IV and S‐V; Supplementary Figures S10–S15). We observed significant phenotypic heterogeneity in individual populations (Supplementary Figure S6), which is in agreement with what has been reported previously in short‐term evolution experiments (Wang et al, 2010). Notably, there was only a small evolutionary trade‐off between adaptation to the medium and to the respective stress, as all strains had similar fitness (Darwinian fitness w=0.96–1.04) in the control (non‐stress) condition (Figure 2). In contrast, medium‐only evolved cells (G500) had significantly lower fitness under other stressors (w=0.51–0.75 relative to the stress‐adapted strains), which argues that adaptation to the medium and to the specific stress are independent processes.
Source data for Figure 2A and 2B [msb201276-sup-0001-SourceData-S1.txt]
Evolutionary cross‐stress protection was observed in various stress pairs. The n‐butanol‐adapted strains (B500) had the highest fitness under hyper‐osmotic conditions, which was similar to the fitness of the osmotic‐adapted strains O500 (w=1.51±0.05 and w=1.39±0.02, respectively). A similar profile was observed under n‐butanol stress where the O500 strains had an increased fitness (w=1.10±0.02), although still lesser than the B500 strains (w=1.58±0.03). Remarkably, all B500 and O500 strains had higher fitness in the acidic environment, outcompeting even the acidic‐evolved strain (P500), a surprising result that clearly argues that similar or higher fitness neighborhoods can be reached through adaptive evolution under different environmental contexts. This result is in agreement with previous reports that adaptation in pH‐related stress has been known to confer fitness effects of low magnitude (Hughes et al, 2007), and lower viability in acidic stress was observed in ethanol‐adapted (Goodarzi et al, 2010), but not isobutanol‐adapted (Atsumi et al, 2010) cell lines. Evolutionary cross‐stress protection was also observed in the case of cells evolved under oxidative stress (H500), which were found to have the same Darwinian fitness as acidic‐evolved cells (P500) under osmotic stress. A negative cross‐stress relationship of lesser magnitude was observed in the case of butanol‐evolved (B500) and osmotic‐evolved (O500) cells under oxidative stress. Interestingly, in competition assays between the B500‐O500 strains in the oxidative environment, the O500 clearly showed increased fitness relative to the B500 strains (1.436±0.030, see Supplementary Text and Supplementary Table S‐V), which further substantiates the negative evolutionary cross‐stress dependency between selection in n‐butanol and oxidative stress.
To investigate whether the phenotypic similarity of the B500 and O500 strains holds for higher stress levels and for different medium composition, we evolved them for another 500 generations and compared their growth in different levels of osmotic and n‐butanol stress (Supplementary Figures S19 and S20). We found that both sets of strains demonstrated almost linear correlation between stress tolerance and number of generations in the stressful environment. Interestingly, osmotic‐adapted cells exhibited a striking trade‐off (sixfold lower density compared with wild‐type) when the medium was switched to LB under the same stress conditions, in contrast to the butanol‐adapted cell lines.
Identifying the genetic basis of adaptation
To elucidate the genetic basis of the observed behavior, we used high‐throughput sequencing to resequence the genomes of selected clones from the adapted cell lines (Figure 3; Supplementary Tables S‐VI–S‐VIII; Supplementary Figure S16). Although we did not resequence the ancestral genome, its reconstruction from the sequences of the five adapted cell lines revealed seven novel genetic variations to the MG1655 genome (Hayashi et al, 2006), six of which were confirmed independently (Freddolino et al, 2012) (Supplementary Table S‐IX). The overall number of fixed mutations and type of observed genomic changes is consistent with our expectations from previous laboratory evolution experiments (see Supplementary Results and Supplementary Table S‐X).
The glucose‐adapted strain G500 showed four SNP mutations in the ynfL, ybcS, pykF and rpoB genes (Figure 3; Table 1). From these, the ynfL mutation is outside the coding region, and the mutation in the ybcS gene may be neutral or inactivating as the latter encodes a prophage lysozyme that is implicated in cell lysis (Srividhya and Krishnaswamy, 2007). Mutations in the pykF gene are known to be beneficial for M9‐glucose minimal media adaptation, as a previous study has shown that an IS insertion into the pykF gene has a significant positive effect in glucose‐limited media (Barrick et al, 2010). Similarly, mutations on the β subunit of RNA polymerase (rpoB gene) have been shown to confer higher growth rates and efficiency in carbon source metabolism than the wild‐type strain (Conrad et al, 2010). In fact, we identified rpoB mutations in five out of the six strains (with P500 being the exception). Previous results highlight the importance of RNA polymerase subunits as targets during microbial laboratory selection experiments. For example, the ratio of the housekeeping sigma factor σ70 and the stress‐related sigma factor σS is crucial for shifting the balance from nutritional competence to stress protection (Patten et al, 2004; Ferenci, 2008; Klein‐Marcuschamer et al, 2009). In the current study, we identified only mutations in the rpoB gene and no other RNA polymerase subunits. This may indeed indicate that selection towards increased performance in a selected growth medium and an additional stressor can be independent from each other. Sequence analysis in the osmotic‐adapted strain O500 revealed a 12‐fold amplification of a 8.5‐kbp region that includes full copies of the enterobactin‐related proteins fepA, entD, fes, hokE, and ybdZ genes (Supplementary Figure S17A). The three first genes are involved in enterobactin‐iron transport and catalysis (Sansom, 1999) and there is evidence that fepA is upregulated under adversarial environments (Wachi et al, 1999; Muela et al, 2008), although its implication in stress protection is not clear. Additionally, we identified a premature stop codon mutation in proV, a subunit of the proline ABC‐transporter which is essential for the uptake of potent osmoprotectants, such as glycine betaine (Perroud and Le Rudulier, 1985; Lucht and Bremer, 1994) that E. coli does not efficiently synthesize de novo (von Weymarn et al, 2001). The promoter region of yijO, a putative transcriptional regulator with unknown function also bore a SNP mutation, which suggests that this putative regulator may be implicated in stress response.
Interestingly, the butanol‐adapted strain B500 has undertaken an alternative path to alter fepA expression by introducing a mutation in its promoter region (Table I), which affects both the fepA and entD genes. We also found a mutation on the acrA gene that is involved in the acr multidrug efflux system for n‐ and iso‐butanol resistance (Atsumi et al, 2010). A SNP mutation was identified in the marC gene, which was recently identified to also play a role in butanol tolerance (McDermott et al, 2008; Atsumi et al, 2010).
In the H2O2‐adapted strain (H500) we find mutations in rsxD, ccmD, katG and rpoB. The rsxD gene is a member of the rsxABCDGE cluster, encoding the soxR reducing system. The deletion of any member of this gene cluster results in constitutive activation of soxS transcription (Koo et al, 2003), which was found to participate in superoxide removal (Demple, 1996), a finding that suggests an overlap of the hydrogen peroxide and superoxide stress response of E. coli. The ccmD gene is a membrane protein and involved in heme transport and cytochrome c maturation by its interaction with ccmC and ccmE (Ahuja and Thony‐Meyer, 2005). Since H2O2 can destroy heme prosthetic groups and therefore inactivate cytochrome c (Villegas et al, 2000), ccmD changes can influence cellular cytochrome c production and thereby compensate increased cytochrome c turnover in oxidative stress conditions. We also identified a mutation in the promoter region of katG which is an oxyR‐dependent gene that was found to be differentially expressed under H2O2‐induced oxidative stress (Zheng et al, 2001). As such, mutations in its promoter region can influence intracellular catalase levels and further detoxification of hydrogen peroxide.
Finally, in the acidic‐adapted strain P500, one of the insertions affected the evgS membrane sensor kinase, which together with the evgA response regulator has already been described to function in acid and drug resistance in E. coli (Itou et al, 2009). The strain also contains a duplication of a 147‐kb fragment, which includes the acid fitness island region involved in the acid response (Supplementary Figure S17B), and is enclosed between homologous rhsB and rhsA genes of the same transcriptional orientation. Additionally, a deletion in rph, a ribonuclease involved in tRNA processing (Kelly et al, 1992) and an insertion in yjjM, essential for the growth on l‐galactonate (Reed et al, 2006) was found, but the importance of these mutations is unclear. Surprisingly, rpoB mutations were not found in the P500 strain, indicating that adaptation to low pH and to the medium may be under antagonistic pleiotropy, something that is further supported by the fact that the P500 strains are unfit when compared with most of the other strains on M9 glucose media and in the absence of any stressor (Figure 2).
Differential gene expression under cross‐stress protection
To further elucidate the genes and pathways that are implicated in cross‐stress protection, we performed transcriptional profiling (RNA‐Seq, verification by qRT–PCR, see Supplementary Materials) of the five strains under the osmotic stress, an environment where cross‐stress protection was clearly observed. Populations evolved under butanol and osmotic stress showed almost the same fitness in this environment, while the fitness of populations evolved under acidic and oxidative stress was increased under osmotic stress when compared to the control populations. We analyzed the data by measuring statistically significant differential expression (DE) of genes between each stress‐evolved strain and the G500 strain (Figure 4). In the H500, O500, P500 and B500 strains we identified 98, 79, 150, and 8 DE genes, respectively (P‐value <0.05 after correction for multiple hypothesis testing). Among all DE genes in at least one of the stress‐evolved strains (267 genes), several GO terms are statistically over‐represented and include functions related to the ribosome, membrane, protein binding, flagella, peroxidase, RNA biosynthesis, and glycolysis (Figure 5A). We observed a significantly larger number of genes simultaneously DE in multiple strains versus what is expected by a pure chance: 49 versus 4 in 2 strains simultaneously, 8 versus 0.06 in 3 strains, and 1 versus 0.001 in 4 strains (see also Supplementary Results). However, despite the significant number of differentially expressed genes in at least one strain, only 9 DE genes (related to protein, membrane, and RNA binding; Figures 4 and 5B) where shared in at least three out of the four strains. It is a remarkably small overlap (by the absolute gene count) that is consistent with recent observations on acquired stress response in yeast (Berry et al, 2011). The majority of DE genes (53% of all DE genes and 83% of genes that do not belong in an amplified DNA region) are under‐expressed relative to the control in at least one of the experimental strains, which is apparently an attempt to shut down non‐essential pathways in stressful environments. There is a strong overlap between DE genes in H500 and O500 strains (52 genes, P‐value <0.05), all of which are downregulated and they encode for parts of the transcription/translation machinery or are involved in carbon metabolism (Figure 5E). In the P500 strain, DE genes include protein binding, flagella and membrane‐related genes (Figure 5G and H). The strongly upregulated flagella cluster was unexpected, and it suggests that flagella genes are implicated in acid response, a link that has been hypothesized before (Maurer et al, 2005). Interestingly, B500 expression profile is very similar to that of the G500 control strain with only eight genes being differentially expressed, which is a stark contrast to their fitness dissimilarity in the various stress environments.
Source data for Figure 4A and 4B [msb201276-sup-0002-SourceData-S2.txt]
Several of the differential transcript expression levels relative to the G500 strain can be explained by the resequencing results. The SNP mutation at the promoter of katG in H500 results in its 3.6‐fold over‐expression in that strain (P‐value=2 × 10−8), as does the SNP mutation in the fepA promoter in B500 (2.3‐fold upregulation, P‐value=2.5 × 10−15). An average 2.2‐fold upregulation (P‐value <10−3), of Acid Fitness Island (AFI) genes in the P500 strain is consistent with the average 1.6‐fold DNA AFI amplification that we identified. Similarly, we observe a 10‐fold upregulation of the fepA gene (P‐value <10−130) in the osmotic stress‐evolved clone as expected from the 12‐fold DNA amplification of several genes from the enterobactin cassette that includes fepA. The proVWX operon was found to have a premature stop codon mutation in the proV gene of the O500 strain, and the operon was found to be significantly downregulated (2.7–4.4‐fold, P‐value <10−30), although we have no explanation on what the role of the mutation is and why the operon is downregulated in the O500 strain. Notably, our analysis did not reveal any major regulators other than rpoB, which has been implicated with increased evolvability in mutant cells (Barrick et al, 2010).
Effect of individual mutations on cross stress protection
To associate the mutations identified with the phenotype observed, we reversed individual mutations in the evolved B500 sequenced clone (see Materials and methods) and again performed direct competition assays (Supplementary Figure S22–S25). Since a single B500 clone was used for repair of the mutations, we first assessed its phenotypic profile through competition assays and we confirmed its similarity to the population averages that were used in the previous analysis (Figure 6; Supplementary Figure S27). One striking result was that mutations to the ferric transporter fepA were found to be adaptive in all environments and confer its upregulation (Figure 4; Supplementary Table S‐XII). This is partially explained by the fact that fepA is responsible for medium adaptation, as fepA reverse mutants had a significant disadvantage against the G500 population (w=0.850±0.005) in M9 glucose competition, while the B500 clone mutants without the repair outcompeted the G500 strains in the same environment (w=1.130±0.019, Supplementary Figure S26). However, in osmotic and acidic stress, the effect of fepA extends above whatever is expected from medium adaptation only, as the reverse mutants were significantly less fit than the original B500 clones (w=0.72±0.03 and w=0.70±0.0006, respectively). The marC mutation was validated to be the dominant driver for the increased fitness of B500 cells under n‐butanol stress, with acrA and rpoB mutations contributing significantly lower. Interestingly, both the marC and acrA mutations were advantageous in the case of cross‐stress protection of B500 in the acidic environment (w=0.85±0.03 and w=0.89±0.03, respectively). Finally, the rpoB mutation was found to have a negative impact to the B500 fitness under oxidative stress and is the main reason of the poor performance of these strains under that environment. Similarly, the lack of an rpoB mutation in the P500 strains can explain the increase fitness of the latter when compared with both O500 and B500 strains.
Through integration of sequencing and transcriptional profiles, functional gene analysis, and measured Darwinian fitness, we are able to get a glimpse of the evolutionary trade‐offs and dependencies that underlie microbial stress resistance. Cross‐stress protection emerged at various degrees in all evolved strains, suggesting a ubiquitous phenomenon that may stem from high pleiotropic dependencies in the underlying genes and pathways. The rapidly emerging and reproducible adaptation patterns argue that bacterial populations occupy a genotypic space that allows high phenotypic plasticity within short timeframes. Adaptation to any specific stress and the growth media was found to be generally independent, as only a few trade‐offs were observed between the control and stress‐evolved strains. Surprisingly, our results argue that long‐term adaptation to an environment does not necessarily lead to the fittest phenotype for that particular environment, as it is the case of cross‐stress protection in acidic stress. This result suggests that the environment in which a cell evolves plays a significant role in the fixation rate of advantageous mutations, possibly by altering the selection dynamics.
We found several different types of mutations, such as SNPs in the coding and non‐coding regions that can potentially affect different biological processes. Although rare, such point mutations are rapidly fixed in a population when they confer a significant fitness advantage. Insertion of transposable elements and the amplification of several genetic loci were also observed in various clones. Transposable elements and their random insertion represent a way for a given population to rapidly explore the genetic network space and develop new phenotypes and can be used to decouple evolutionary conserved responses even on a laboratory time scale (Tagkopoulos et al, 2008). In the current analysis, gene amplifications result in increased gene dosage and have an impact to the organism's fitness in a stressful environment but our data also indicate minor fitness trade‐offs upon removal of the stressor and transfer to a less selective environment. Our observations are in accordance with recent theoretical analyses on the role of gene duplication in adaptation and evolution of novel functions (Kondrashov and Kondrashov, 2006; Bergthorsson et al, 2007). Furthermore, our data on low pH evolved populations indicate that the emergence of optimal fitness can be slow under certain circumstances and that improved phenotypes are hardly accessible from a given point in network space. As such, evolutionary adaptation to low pH may represent a suitable model system to further analyze the importance of historical contingency and population size in molecular evolution and neutral network theory (Wagner, 2008).
Two of the strains, O500 and B500, were found to have similar profiles and they exhibited the most profound cases of cross‐stress behavior. A common denominator in both these strains is the upregulation of genes that are related to iron acquisition, either through mutation in the respective regulatory region (B500) or through DNA amplification (O500), which can partially explain the observed phenotypic similarity in the osmotic and acidic environments. Indeed, the fepA mutation was found to be the dominant factor behind the high B500 fitness in both these environments, as reversal of this mutation significantly decreased the fitness of the B500 mutants. This link that we establish in this work is in support of indirect prior evidence that iron‐related pathways are generally over‐expressed under osmotic stress (Snyder et al, 2004; Hagan and Mobley, 2007; Muela et al, 2008; Nielubowicz et al, 2008) and the high correlation of iron‐regulation and acid response pathways in Helicobacter pylori (Bijlsma et al, 2002); Salmonella typhimurium (Hall and Foster, 1996), Shigella flexneri (Oglesby et al, 2005), Corynebacterium glutamicum (Jakob et al, 2007) and E. coli (Zhu et al, 2002). The role of iron transport in n‐butanol and more generally in stress tolerance is not clear although it has been shown that iron is critical to the control of genes implicated in lipopolysaccharide biosynthesis, which in turn is the most prominent mechanism of toxic tolerance (Hagiwara et al, 2004; Reyes et al, 2012). Our results show that mutations in iron‐related genes are a common genetic factor that drives bacterial tolerance across multiple stresses.
Both osmotic‐evolved and butanol‐evolved cells performed poorly under oxidative stress, which is in contrast to their superior performance in all other stresses. One hypothesis may be that the upregulation of iron‐related genes in both strains may lead to an excess of iron inside the cell that can be toxic under oxidative stress (Touati et al, 1995; Andrews et al, 2003) as it facilitates the creation of reactive hydroxyl radicals that damage nucleic acids. While our mutation reversal experiments support this hypothesis, the effect of the iron transporter fepA was found to be secondary to the effect of the rpoB mutation, which was maladaptive in the case of oxidative stress. In addition, since the two strains have very different mechanisms to alter iron‐related cells (i.e., mutation of the fepA promoter in B500 versus 12 × amplification in O500), the conservation of the regulatory control in the O500 can attribute to the increased fitness of the O500 to B500 strains when compared in oxidative stress. This is also supported from the transcriptional analysis: fepA was found to be significantly downregulated in G500 strains under oxidative stress (qRT–PCR, 3.8 × downregulation), upregulated in O500 cells (8.2 × upregulation), and strongly upregulated in B500 strains (16 × upregulation). Another interesting finding of the transcriptional analysis is the strong upregulation of the flagella cluster in the acidic environment. It would be interesting to further investigate the importance of this finding in the organism's physiology within its ecological context, for example as a mechanism to maintain pH homeostasis during its transition through the mammalian GI tract.
Given its abundance, it is important to further investigate what are the underlying causes of evolutionary cross‐stress protection and trade‐offs for different stress combinations. If specific stresses have similar effects at a cellular level, then a high degree of cross‐stress protection can be expected since responses to these stresses have to protect the cell versus common perturbations. Alternatively, the environmental context under which cells have evolved may be responsible for cross‐stress protection or trade‐off, even in the absence of common underlying mechanism. Indeed, E. coli and baker's yeast were recently shown to exhibit cases of cross‐stress protection when stresses were applied in the order that exist in their natural environment, while this coupling was lost quite rapidly when adapted in environments with the inverse dynamics (Tagkopoulos et al, 2008; Mitchell et al, 2009). Concomitantly, the absence of cross‐stress protection or strong trade‐offs may imply that the respective stresses did not co‐exist in the organism's natural habitat, nor was the latter exposed to rapid fluctuations across these dimensions. The expansion of laboratory evolution experiments beyond what we perceive to be their natural habitat parameters is relevant to other areas, including the evolution and treatment of food‐borne and nosocomial pathogens. Due to their fast rate of evolutionary adaptation, bacteria pose a major problem in those areas and a better understanding of their adaptive potential in natural and non‐natural stress conditions can lead to techniques that aim to control their proliferation. Furthermore, besides its significance in the development of sterilization techniques, analyses similar to what presented here may also greatly impact the development of microbial strains in industrial biotechnology. The biotechnological production of biofuels and other organic substances (including non‐natural metabolites) depends on microbial cells as biocatalysts. In most cases such systems depend on bacterial growth and microbial cells can suffer stress during cultivation, including variations in pH, aeration and temperature (e.g., during cultivation in large bioreactors due to improper mixing of the culture broth). Additionally, cells are often engineered to produce chemical compounds that occur only in very small quantities or not at all in their natural habitats. As such, the interplay and compatibility of ubiquitous stress parameters and artificial stresses is of great interest for the optimization of bacterial growth in complex industrial environments. As our knowledge of bacterial physiology and ecology continues to expand, so will our understanding of the underlying causes of cross‐stress potential and behavior. Until then, evolutionary experiments such as the one presented here will provide a unique perspective on the emergence of complex microbial behavior.
Materials and methods
Bacterial strains and growth media and conditions
E. coli MG1655 and MG1655 ΔlacZ strains were used in all experiments of this study. The inclusion of ΔlacZ mutants allowed us to perform competition assays using X‐Gal and IPTG staining. The neutrality of the ΔlacZ mutation was confirmed in all environments prior to adaptation (Supplementary Table S‐I; Supplementary Figure S4). Minimal M9 salt medium with 0.4% (w/v) glucose as carbon source was used for evolutionary adaptation, growth tests, and competition assays. All cultivations were performed in 10 ml volume in 125 ml shake flasks at 37°C on an orbital shaker at 150 r.p.m.
Evolutionary adaptation by serial dilution experiments
Adaptive evolution experiments were performed by daily serial dilutions after 24 h of growth in fourfold replication (2 × E. coli MG1655 and 2 × E. coli MG1655 ΔlacZ) for each growth condition. For each of the daily transfers, the OD600 was monitored and cultures were diluted (1:500) in fresh growth medium yielding ∼7–9 generations per day. Serial dilutions were performed for a total of 500 generations.
E. coli strains were evolved in M9 salts medium with glucose as the sole carbon source and the following stressors: osmotic stress (0.3 M NaCl, O500 strains), acidic stress (pH 5.5, P500 strains), oxidative stress (100 μM H2O2, H500 strains), n‐butanol stress (0.6% n‐butanol, B500 strains), and control (no‐stress, G500 strains). The OD600 of each culture was measured each day before the daily transfers to ensure that the estimated nine generations per day were reached (Supplementary Material).
Growth curves of the ancestral and evolved cell lines were performed in the growth media described above. For growth curves, the starter cultures of all strains were grown and therefore adapted (∼7–9 generations) to the medium to be tested, for 18 h. Cultures were started at OD600=0.1 or 0.05, depending on the growth condition. Samples were taken in regular intervals and OD600 was measured on a biophotometer (Bio‐Rad). Two independent replicate growth tests were performed for each strain and population.
For direct competition assays, the two competing populations (E. coli MG1655 and E. coli MG1655 ΔlacZ) were grown under adaptive conditions overnight and inoculated at an approximate 1:1 ratio at a starting OD600=0.004. After 24 h of growth, cultures were diluted 1:500 in fresh medium and competed for additionally 24 h. Samples were taken at 0, 8, 24, and 48 h. Two independent replicates were performed for each individual competition. Cell counts were determined on LB agar plates containing 0.25 mM IPTG (Isopropyl β‐d‐1‐thiogalactopyranoside) and 40 μg/ml X‐gal (bromo‐chloro‐indolyl‐galactopyranoside). Plates were incubated at 37°C overnight. Darwinian fitness (W) was calculated as described in the Supplementary Methods section.
Maximum stressor concentration
Starter cultures of evolved and ancestral populations were grown in M9 medium at 37°C overnight and a main culture (M9 and LB medium with different stressor concentrations) was inoculated at OD600=0.1. Growth potential was monitored by measuring optical density after 6 and 24 h. Maximum stressor concentration is referred to as the concentration of NaCl and butanol in the growth medium, which exceeded the concentrations that were used for the initial evolutionary adaptation experiment to test whether the maximum level of abiotic stressor was increased as compared with the ancestral strains.
Genome resequencing of evolved clones
Selected clones were grown on LB medium overnight and genomic DNA (gDNA) was isolated using a Wizard Genomic DNA Purification Kit (Promega). gDNA was quantified photometrically and 7 μg were used for fragmentation in a bioruptor (Diagenode) sonication apparatus. DNA was blunt ended using T4 DNA polymerase and Klenow fragment. A tailing was performed using Klenow fragment (exo). After A tailing DNA linkers were ligated using fast ligase. PCR enrichment and DNA size selection was performed on 2% agarose gels. DNA purification steps were performed using a NucleoSpin Extract II purification kit (Machery and Nagel). DNA libraries were sequenced by a paired end run with 85 bp read length on an Illumina Genome Analyzer GAII. The shotgun sequence data was assembled, using the IDBA assembler (Peng et al, 2010), after preliminary filtration and quality‐trimming of the Illumina reads using SGA (Dinh and Rajasekaran, 2011). Scaffolding was performed with SSPACE script (Boetzer et al, 2011). Parameters for the assembly and scaffolding were selected based on the estimations provided by an automated pipeline (Darling et al, unpublished data). Scaffolds were reordered using Mauve Contig Mover (Rissman et al, 2009) and aligned to the reference E. coli strain K12 substrain MG1655 genome (Hayashi et al, 2006) using progressive Mauve genome aligner (Darling et al, 2010). SNPs and variances between sequenced strains and the reference E. coli genome were obtained for each strain by analyzing both: (i) shotgun reads mapped to the reference genome using BWA (Li and Durbin, 2010) and SAM tools (Li et al, 2009) and (ii) the de novo assembly followed by alignment to the reference genome. Variances that are found in all six strains are attributed to the ancestral strain; unique variances are the result of the independent evolution of the strains.
Gene expression analysis
For transcriptional profiling, the resequenced clones were grown in M9‐gluc medium with 0.3 M NaCl (osmotic stress medium) to mid‐log phase (OD600 approx. 0.8). In all, 1 ml aliquots were harvested and mixed with 0.5 ml Phenol/ethanol (5% (v/v) Phenol in EtOHabs). Cells were collected by centrifugation and stored at −80°C until use. RNA was extracted using an RNeasy kit (Qiagen) and after first‐ and second‐strand cDNA synthesis, linker ligation, and size selection subjected to shotgun sequencing using single read 40 bp read length runs on an Illumina Genome Analyzer GAII. llumina reads for RNA‐Seq libraries were mapped to the reference E. coli strain K12 substrain MG1655 genome (GenBank accession no. U00096.2) by BWA toolkit (Li and Durbin, 2010). Differentially expressed genes were identified by processing raw counts of mapped reads with edgeR library. Gene ontology analysis was performed by goseq R package. Differentially expressed (DE) genes in all stress‐evolved strains were found relative to the expression levels in the reference G500 strain using the edgeR R package. Genes with BH (Benjamini and Hochberg) adjusted P‐values below 0.05 threshold were selected as DE genes. Genes from the amplified regions in O500 (12‐fold) and P500 (twofold) strains are significantly over‐expressed relative to the reference: average log Fold Change for concentrations is +3.6 and +0.8 (both with P‐values <10−20), respectively. Full list of genes DE in at least one of the strains and sorted by the expression patterns is shown in Supplementary Table S‐XI. The RNA‐seq data set can be obtained from the GEO repository, record no. GSE39926.
Mutation repairs in B500‐evolved strain
To perform single mutation repairs, DNA fragments from wild‐type E. coli MG1655 were used to replace each of the mutations in the B500‐evolved strain by P1 transduction (Miller, 1992) using a chloramphenicol resistance gene as a selectable marker. PCR‐amplified chloramphenicol resistance cassettes were introduced into wild‐type E. coli chromosome in sites corresponding to genomic regions close to the mutation sites present in the B500 strain (described in Table I) by linear DNA transformation (Poteete and Fenton, 1984, 2000; Murphy, 1998; El Karoui et al, 1999; Datsenko and Wanner, 2000). The newly transformed strains carrying one chloramphenicol cassette linked to one of the mutations were used as donors to replace the B500 strain mutated sequences. New B500 mutants repaired a single mutation site, while still carrying three of the original four mutations. The repaired fragments were verified by Sanger sequencing.
Fitness effect of individual mutations
Repaired B500 mutants were competed against the original B500‐evolved strain on M9 salts medium with glucose under several conditions. The two competing strains were grown overnight in 2 ml M9 salts medium with glucose, washed and inoculated at an approximate 1:1 ratio for direct competition under adaptive conditions. After 24 h, they were diluted 1:50 into fresh medium and competed for an additional 24 h. Two independent replicates were performed for each individual competition. Culture samples at 0, 24, and 48 h were diluted and plated in LB medium and LB medium with 20 μg/ml chloramphenicol, and incubated at 37°C overnight. The number of colonies present on the LB medium represents the total number of all cell types in the population. The number of colonies on the LB chloramphenicol medium represents the total number of chloramphenicol resistant cells in the population (repaired mutant). The difference between colony numbers from the two plate types will show the total number of chloramphenicol susceptible cells in the population (original B500 strain). From this experiment, we determined the relative fitness of each repair mutant under different growth conditions relative to the original B500 mutant. E. coli MG1655 strain with an added chloramphenicol resistance gene was competed to wild‐type E. coli MG1655 (chloramphenicol sensitive) as a control.
Verification of mutations and transcriptional abundance
All mutations were independently verified by sequencing. To test gene expression fold changes, evolved E. coli MG1655 cells were grown in M9 salt glucose medium. The overnight cultures were diluted 1:10 in fresh M9 salt glucose medium with different stress (no stress, 0.3 M NaCl, 0.6% n‐butanol, or 100 μM H2O2) and grown at 37°C with shaking (150 r.p.m.) to the log phase. Total RNA was prepared according to the manufacturer's instructions (RNeasy mini kit, Qiagen). cDNA was synthesized from 1 μg of total RNA by using the iScript cDNA synthesis kit (Bio‐Rad). Cells were grown in Luria Broth (LB) media overnight, and genomic DNA was isolated via the Wizard Genomic DNA Purification Kit (Promega). Samples were then analyzed by Chromo4 instrument (Bio‐Rad). Primers were designed using the Primer Quest (IDT DNA). Primers were made for the target genes and σ70 gene, which was used as reference. qPCR was performed in 50 μl mix containing 1 × SsoAdvanced SYBR Green Supermix (Bio‐Rad), 40 ng genomic DNA or cDNA, and 5 pmol primers. Experiments were done in triplicate and qPCR data were analyzed by Opticon Monitor 3 (Bio‐Rad).
We thank Aaron Darling, Jonathan Eisen, Marc Facciotti, Artyom Kopp, John Roth, and Saeed Tavazoie for the manuscript comments and helpful discussions. This work was supported by the UC Davis Opportunity Fund and DoD award W911NF1210231‐0.
Author contributions: MD performed the strain evolution, strain characterization, and library preparation for genome sequencing and RNA‐Seq. VM performed the data analysis, the genome assembly and RNA‐Seq. SQS repaired the mutations and performed competition assays for validating the fitness effect of the mutations. JP performed qRT–PCR validation experiments. MD, VM, SQS and IT evaluated data and prepared the manuscript. IT conceived of the study and supervised all aspects of the project.
Conflict of Interest
The authors declare no conflict of interest.
Supplementary Tables SI–XII, Supplementary Figures S1–27 [msb201276-sup-0001.pdf]
Supplementary Table S‐VIII
De novo assembly breaks analysis file [msb201276-sup-0002.xlsx]
Supplementary Table S‐XI
Log‐fold change and p‐values for transcriptional profiling (RNA‐Seq) data [msb201276-sup-0003.xlsx]
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2013 EMBO and Macmillan Publishers Limited