High-throughput laboratory evolution of Escherichia coli under multiple stress environments

Bacterial cells have a remarkable capacity to adapt and to evolve to environmental changes. Although many mutations contributing to adaptive evolution have been identified, the relationship between the mutations and the phenotypic changes responsible for fitness gain has yet to be fully elucidated. For a better understanding of phenotype-genotype relationship in evolutionary dynamics, we performed high-throughput laboratory evolution of Escherichia coli under various stress conditions using an automated culture system. One measure of phenotype, transcriptome analysis, revealed that the expression changes which occurred during the evolution were generally similar among the strains evolved in the same stress environment. We also found several genes and gene functions for which mutations were commonly fixed in the strains resistant to the same stress, and whose effects on resistance were verified experimentally. We demonstrated that the integration of transcriptome and genome data enables us to extract the mechanisms for stress resistance. Author summary Understanding the relationship between phenotypic and genetic changes is a fundamental goal in evolutionary biology, which can provide insights into the past and future evolutionary trajectories. Evolution of microorganisms in a laboratory has been the primary approach to clarify the mappings of phenotypic and genotypic changes. Here, we performed high-throughput laboratory evolution with bacteria using an automated culture system, to quantify phenotypic and genotypic changes occurred under various stress conditions. We identified various stress-specific gene expression changes and mutations, and contributions of them to fitness gain were validated. These results demonstrated that the integration of phenotypic and genotypic changes makes it possible to extract the mechanisms for stress resistance evolution, which will contribute to bioengineering applications.


Introduction
Laboratory evolution of microorganisms is a powerful approach for elucidating the nature of evolutionary dynamics [1,2]. Recent advances in measurement technologies including high-throughput sequencing have enabled us to quantify phenotypic and genotypic changes during laboratory evolution, which have provided valuable information on the mechanisms and principles of adaptive evolution [3][4][5]. The impact of laboratory evolution has extended beyond the field of evolutionary biology into engineering and medicine. For example, by using laboratory evolution approaches, a number of candidate mutations that can contribute to antibiotic resistance have been identified, which sheds light on how to control the emergence of antibiotic resistant strains [6][7][8][9].
Laboratory evolution has also become a widely used tool for bioengineering applications [10][11][12][13], in order to generate cells with improved growth, production titer, and stress tolerance, which are essential to improved industrial production by microorganisms.
A central problem in laboratory evolution of microorganisms is identification of phenotypic and genotypic changes responsible for an observed fitness gain in an evolved strain.
Such identification is often difficult, even after whole-genome sequencing analysis of the evolved strains. In some cases, a large number of genotypic changes, mutations, are fixed in evolved strains, and only a portion of such mutations contribute to the fitness gain. Furthermore, various different mutations can cause a similar phenotypic change [5,9,14], and one mutation can affect various cellular functions. The many-to-many relationships between phenotypic and genotypic changes make identifying the mechanisms of fitness gain challenging.
One strategy to overcome such complexity in the analysis of laboratory evolution is to integrate changes of genomic sequence with various levels of phenotypic data. Such phenotypic data include those from omics analysis such as transcriptome and metabolome. Based on such phenotypic data, one can extract correlated phenotypic changes to the observed fitness gain, which facilitates identification of mutations responsible for the phenotypic and fitness changes. For example, analysis of laboratory evolution of antibiotic resistant E. coli revealed that the change in antibiotic resistance can be quantitatively predicted by the transcriptomic change during adaptive evolution, which led to identification of genomic mutations contributing the resistance [9]. Furthermore, cross-stress protection, a phenomenon whereby stress resistance is acquired to one stress after adaptive evolution to another stress, also provided valuable information on which phenotypic and genotypic changes related to the stress resistances when analyzed by integrating transcriptome and genome resequencing data [15]. Integration of various levels of phenotypic data with the change of genomic sequence is the future of investigating evolutionary dynamics and will be an essential method for the application of laboratory evolution in medical and engineering fields.
Given the importance of the integration of multiple levels of phenotypic and genotypic data in laboratory evolution, we performed high-throughput laboratory evolution of E. coli cells under 11 different stress conditions and then quantified phenotypic changes, specifically transcriptome expression changes and growth rate, and genotypic mutations of independently evolved strains. Transcriptome analysis revealed that, the expression changes which occurred during the adaptive evolution were generally similar among the strains evolved in the same stress environment. We analyzed the correlation between these expression changes and mutations commonly fixed in the strains resistant to the same stress. This analysis revealed the mechanisms for some of the stress resistances, which were validated by introducing the relevant mutations into the genome of the parent strain. The integration of the multilevel phenotypic and genotypic data enables us to understand a wide spectrum of molecular mechanisms in E. coli stress resistances.

Laboratory evolution under 11 stress conditions and analysis of cross-protection
We selected 11 stresses that cover a wide range of environmental stresses (Table 1). E. coli MDS42 cells [16] were cultured in 200 μL of M9 synthetic medium [17] with constant concentration of stressors. The concentrations of these stressors were set to levels which initially decreased the specific growth rate by approximately one half of the non-stress growth condition. Every 6 hours, a fraction of the cells was transferred into fresh medium with the stressor. The transfer volume was adjusted to keep the final cell concentration below a threshold and to maintain the culture in exponential growth phase. The specific growth rate was quantified by initial and final cell concentration, which was used as the measure of fitness under these stress environments (Fig 1a). To evaluate the reproducibility of the evolutionary pathways, for each stress, five independent culture lines were propagated in parallel. For this high-throughput laboratory evolution, we used a custom-developed automated system for laboratory evolution [18], by which we can maintain hundreds of independent culture series in a fully automated manner. After 906 hours of propagation, we observed significant increases in specific growth rate in all 11 stress environments (Fig 1b to 1m).
In addition to the cultures with stressors, we maintained the cells in the synthetic medium without adding stressors (Fig 1n). In these cultures, without stress, slight increases of specific growth rate were also observed, although these increases were significantly smaller than those in the stress conditions. For all resistant strains, we confirmed maintenance of the stress resistances by reintroduction after cultivation in the absence of the stressors for at least 60 generations, indicating that the phenotypes of stress resistance were stable.
To explore evolutionary cross-protection, for each obtained resistant strain, we measured the specific growth rate under all 11 stress conditions. In Fig 2, the relative growth rates of the resistant strains to the parent strain are presented (All the growth rates are presented in S1 Table).
The results demonstrated that the stress resistant strains often show cross-protection and hyper-sensitivity to multiple stresses. A clear symmetric cross-protection was observed among NaCl and KCl stress resistant strains, i.e., NaCl resistant strains exhibited significant KCl resistance, and vice versa. This fact suggested that, at least a part of the mechanisms of the resistances were shared among NaCl and KCl resistant strains. For other cases, the cross-protections were asymmetric, which suggested that the mechanisms of these resistance and hyper-sensitivity were independent, or there are hierarchical relationships between the resistance/hyper-sensitivity mechanisms. For example, methylglyoxal (MG) resistant strains exhibited resistance to lactate (Lac), while Lac resistant strains did not significantly affect MG resistance. This asymmetric cross-protection might be explained by a hierarchical structure of the resistance mechanisms. To degrade MG in metabolic pathways, one possible intermediate is Lac [19], and thus MG resistance strains might need to have Lac resistance at the same time. In contrast, MG is not involved in Lac the degradation pathway, and thus Lac resistance would not necessarily be accompanied with MG resistance.

Changes in expression profiles in stress resistant strains
To elucidate the mechanisms of the various stress resistance phenotypes, we performed transcriptome analysis of all the resistant strains we created by using microarray experiments (all transcriptome data are presented in S2 Table). The results demonstrated that the expression changes which occurred during the long-term cultivation were generally most similar among the strains evolved in the same stress environment, as shown in the hierarchical clustering analysis transcriptome data (Fig 3a). To characterize the similarity in transcriptomic changes in these resistant strains, the expression changes of representative transcriptional factors (TFs) are presented in Fig 3b. As shown, these TFs exhibited similar expression patterns within each stress resistant clade. For example, expression levels of marR encoding a repressor protein for the marRAB operon were commonly down-regulated in BuOH and Na 2 CO 3 resistant strains (S1a Fig). The marRAB operon is involved in resistance to several stresses including antibiotics and organic solvents [20], and the inactivation of marR has been shown to increase the resistance of E. coli to multiple antibiotics and organic solvents by up-regulation of the AcrAB-TolC multidrug efflux pump [21,22]. Our results also showed a negative correlation between marR and acrAB expression levels (S1b Fig), which suggests up-regulation of acrAB contributed to the observed BuOH and Na 2 CO 3 resistances. In contrast, marR expression was up-regulated in Crotonate (Cro) resistant strains. Although this characteristic up-regulation might be involved the resistance, the mechanism remains unclear. For another example, hns encoding a histone-like protein H-NS were up-regulated in all CoC, Lac, Mal, MG, and CPC resistant strains (S2a Fig). H-NS is a global regulator involved in the adaptation to various stress environments, whose change in expression level affects hundreds of genes. Among these target genes, we found that rcsD expression levels in the resistant strains clearly showed a negative correlation to hns expression level (S2b Fig). rcsD is known to be involved in acid stress resistance [23], and thus the correlation suggests the contribution of hns regulation to stress resistance through rcsD expression changes. The observed similarity of expression changes in the resistant strains, as represented by the TFs expression patterns in Fig 3b, strongly suggests that such expression changes contributed to the acquisition of the resistances.
To further analyze the correlations between two elements of phenotype, transcriptome changes and stress resistance, we used a simple, previously published [9] mathematical model to predict the resistances using the obtained gene expression profiles (see Materials and Methods for details). Briefly, we assumed the changes in the specific growth rates under stresses are represented by a linear combination of log-transformed expression changes during the long-term cultivation.
Then, we sought the optimal number and combination of genes which have the highest predictive accuracy by using cross-validation and genetic algorithm. As result, we found that the combinations successfully screened genes whose expression changes are highly correlated to resistance acquisition, which can contribute to highly accurate descriptions of complex evolutionary dynamics.

Mutations fixed in stress resistant strains
The mutations identified in the resistant strains are shown in Fig 4, and the detailed information about the mutations is presented in S3 Table. Less than 10 mutations were fixed in each of the resistant strains. We also sequenced two strains maintained without addition of stress (Fig 1n), and found that relatively small number of mutations were fixed in these control strains, as shown in Table S3.
We found several genes and gene functions for which mutations were commonly fixed in the resistant strains, suggesting contributions by these mutations to the resistances. To verify the possible contribution of these mutations to the stress resistance, we introduced some of such commonly identified mutations (highlighted by a red character in Fig 4) into the genome of the parent strain, and quantified the change of growth rate under corresponding stresses ( Fig 5). Below, we discuss some examples of the relationship between resistance acquisition and phenotype/genotype changes. A full description of the discussions is presented in S1 Text.
The first example of the common mutations is that of to the proU ABC transporter which was identified in all NaCl resistant strains and 4 KCl resistant strains (Fig 4a and 4b). All of these mutations caused frame-shifts, suggesting that they disrupted the activity of the proU transport system. The proU operon (proVWX) encodes a binding protein-dependent transport system which is essential for the uptake of osmoprotectants such as glycine betaine and is known to be up-regulated in response to osmotic stress [24]. Our expression analysis showed that, the expression of proU operon was significantly up-regulated (> 50 times) in response to the initial NaCl stress addition, however after long-term cultivation the expression levels decreased when the cells acquired NaCl resistance (S5 Fig). These results suggested that, the programmed up-regulation of proU operon in response to osmotic stress is not beneficial in an environment without osmoprotectants such as that which we used. Instead, the up-regulation of proU operon suppresses the cellular growth presumably due to its energy consumption. Thus, the disruption of proU transport activity can contribute to the active cellular growth under NaCl stress. To support this hypothesis, we introduced a frame-shift mutation identified in one of the NaCl resistance strains into the parent strain using site-directed mutagenesis, and confirmed that this mutation increased the growth under NaCl stress ( Fig 5). ydhM expression controls the expression of gloA, which is involved in the degradation pathway of MG to lactate and is known to contribute to MG resistance [19]. There is a clear correlation between The detailed descriptions for the other common mutations found in the resistant strains are presented in S1 Text. We demonstrated that the identified common mutations generally contributed to resistance acquisition, as shown in Fig 5. It should be noted that, to understand the mechanism how such mutations affected the stress resistance phenotype in addition to the identification of genomic mutations, the information on expression changes, transcriptome phenotype, during the adaptive evolution was also highly valuable, as in the cases of proU and ydhM mentioned above.
It is worth noting that, although in some cases the data suggested simple causal relationships between genomic mutations and phenotypic changes, the relationships between them were not always a simple one-to-one correspondence. As shown in Fig

Conclusion
Laboratory evolution experiments combined with organism-wide phenotype and genotype analyses enable us to understand mechanisms of adaptive evolution. In this study, we demonstrated that the high-throughput system for laboratory evolution we developed can create resistant strains of E. coli to various stress environments by long-term cultivation. Using this system, we can cultivate up to 44 microplates (96 or 384 wells) simultaneously, and thus more than 10,000 independent culture series can be maintained in a fully automated manner. This system allows us to trace evolutionary dynamics under various environmental conditions and initial conditions (e.g., all E. coli strains in single-gene knockout library[30]), with numerous replicate experiments.
In this system for laboratory evolution, we maintain cells under relatively small culture scale (i.e., 200 μL). One might question whether such small culture scale might result in small population size, which causes fixation of random drifts and accumulation of neutral mutations. Our data suggested that this is not the case. We confirmed that the population size was maintained more than 10 5 cells. The ratio of synonymous to non-synonymous substitutions over all resistant strains was relatively small in comparison with that of neutral mutations, suggesting that the fixation of a majority of substitutions were driven by evolutionary selection pressure. The fact that, for many stresses, resistant strains to a common stress shared common mutations also supports the idea of evolutionary selection.
The genome-wide expression and resequencing analyses showed common phenotypic and genotypic changes in the resistant strains to an identical stress. It is worth noting that, our results demonstrated that combined use of transcriptome and genome resequencing analysis greatly accelerated our interpretation about how E. coli strains acquired the stress resistance. For example, mutations related to the proU operon found in all NaCl resistant strains were thought to cancel the deleterious regulatory program. Using only sets of fixed mutations in resistant strains often makes it difficult to establish the causal relationship with stress resistance. Integrating both gene expression and genotype data can overcome this problem and it allows a highly accurate description of the mechanisms responsible for stress resistance.
The adaptive evolution to environmental changes is a phenomenon that involves changes in genome, transcriptome, metabolome, and so on, meaning a complex interaction network is involved. One possible strategy to understanding such complex dynamics is to analyze large-scale data for each hierarchical layer and then to integrate the analyses to extract the essential components for the phenotypic and genotypic changes. The present study is one of the first to achieve genome-wide analysis of evolved strains obtained by high-throughput laboratory evolution, and we succeeded in extracting phenotypic/genotypic changes responsible for stress resistance. Such large-scale data of laboratory evolution will, we believe, help us to unveil principles of evolutionary dynamics and provide valuable information on rational design of industrially useful microbial strains.

Laboratory evolution.
The IS elements-free Escherichia coli strain MDS42 [16] was purchased from Scarab Genomics and used as the parent strain of the laboratory evolution. The use of the IS elements-free strain can make the resequencing analysis reliable since the determination of the precise position of IS element insertions is often difficult by using high-throughput sequencers. For serial transfer culture, we used 200 μL modified M9 medium [17] including 5g/L glucose as a carbon source. The cells were cultured in 96-well microplates (Corning Inc. 3595) with shaking at 300 rotations/min at 34 °C. All cultivations were performed by the automated culture system [18] (Fig 1a) which consists of a Biomek® NX span8 laboratory automation workstation (Beckman Coulter, Tokyo, JP) in a clean booth connected to a microplate reader (FilterMax F5; Molecular devices, Tokyo, JP), shaker incubator (STX44; Liconic, Mauren, LI), and microplate hotel (LPX220, Liconic, Mauren, LI). The movie of this automated culture system for laboratory evolution is available at youtube (https://www.youtube.com/watch?v=4k6qCN7ppsk). We diluted the cells into a fresh medium every 6 hours. The cells were maintained in the exponential growth phase by adjusting the initial cell concentration of each dilution to keep a final cell concentration of less than 10 7 cells per well as measured by optical density at 620 nm (OD 620 ). Before laboratory evolution under the stress conditions, we cultivated the cells without adding stressors for 96 hours (approximately 90 generations) to adapt them to M9 medium. The specific growth rate was calculated based on the initial and final cell concentrations of the each dilution. The cells after the evolution experiments were single-cloned and were stored as glycerol stocks at −80 °C and used for further analysis. The quantifications of specific growth rates of evolved or genetically manipulated strains (Fig 2 and Fig   5) were performed after 60 hours preculture (approximately 30 to 60 generations) while maintaining exponential growth phase. For the preculture, the same culture conditions as the laboratory evolution with the corresponding stressor were used.
Transcriptome analysis by microarray.
The cells were precultured for 60 hours under the same condition with the laboratory with or without the corresponding stressors. Then, 5×10 7 cells in the exponential growth phase were killed immediately by the addition of an equal volume of ice-cold ethanol containing 10% (w/v) phenol.
After that, the cells were harvested by centrifugation and stored at −80 °C before RNA extraction.
Total RNA was isolated and purified from cells using an RNeasy mini kit with on-column DNA digestion (Qiagen, Hilden, Germany) in accordance with the manufacturer's instructions. The quality of the purified RNA was evaluated using Agilent 2100 Bioanalyzer with an RNA 6000 Nano kit (Agilent Technologies). The purified RNAs were stored at −80 °C prior to transcriptome analysis.
Microarray experiments were performed using custom designed the Agilent 8×60K array for E. coli  Table. In that table, we also present the biological triplicate data to check the reproducibility of the analysis, i.e., expression data obtained from different cultures of the parent strain without addition of stress. So as to use only quantitatively reliable data, genes with low expression levels (less than 100 a.u. in any strain) were excluded from the following analysis (approximately 60 % of genes remain). We confirmed that after the removal of low expression genes, more than 99 % of expression ratios between the biological repetitive data were within the range 1/1.3 to 1.3.

Predicting stress resistance by gene expression levels.
To examine the contribution of the gene expression changes to the stress resistances, we constructed a simple mathematical model to predict the resistances using the obtained gene expression profiles [9]. Here, we assumed that the stress resistances, quantified by the change of growth rate after stress addition, are determined as a function of gene expression levels and neglected any direct effect of the mutations on the resistance. Furthermore, for simplification, we neglected non-linear effects and cross terms of the gene expression changes. Thus, we assumed the following simple linear model to predict the growth rate change by the expression levels of N genes: ∆ is the growth rate changes in the j th strain for the k th stress, is the log 10 -transformed expression level of the i th gene in the j th strain after standardization to zero mean and unit variance, and and are fitting parameters. Since the number of genes used in this analysis is much larger than the number of fitness data, the use of all transcriptome data for the fitting resulted in overfitting leading a meaningless prediction of the fitness. In order to avoid such overfitting and to investigate the number of genes with the highest predictive accuracy, we used the cross-validation method, in which the fitness data were separated into training data used for the parameter fitting and test data used to verify the prediction accuracy. When N was large, the prediction accuracy for the test data became small due to overfitting, while the accuracy became small when N was small since the linear combination of expressions was insufficient to represent changes of the fitness data. In this analysis, we used a 4-fold cross validation method. That is, the resistant strains were randomly partitioned into 4 equally sized subgroups; 1 subgroup was used as the test data set for validation and the remaining 3 subgroups were used for the fitting.
Expression levels of genes in the same operon are generally correlated, which can cause problems in the gene selection procedure. Thus, in each operon, we selected a gene having the highest average expression level over all samples and use it for the fitting, while we discarded the data of other genes. Furthermore, we excluded data with relatively small expression changes among the resistant and parent strains, since the expression changes of such relatively unchanged genes dominated the experimental errors. These selection criteria left 413 genes for the analysis.
N genes used for the fitting were selected by a genetic algorithm (GA), in which the correlation coefficient between the predicted and observed growth rate changes of the training data sets was used as the fitness function. As an initial population, 1000 sets of N genes were randomly selected, and the fitness of the sets was calculated. Then, gene sets within the top 5% of the highest fitness were selected as the parent sets of the next generation from which mutant sets were generated by randomly replacing a single gene. We iterated 300 cycles of the mutant sets and selected gene sets with the highest fitness to obtain sets of small number of genes whose expression levels can represent changes in drug resistances and susceptibility. We repeated the selection of gene sets using 10,000 different training data sets prepared by randomly partitioning the total data set to obtain the frequency of genes selected after the GA, as shown in S3b The quality of the sequence data was first assessed with FastX Toolkit 0.0.13.2 (http://hannonlab.cshl.edu/fastx_toolkit), and the raw reads were trimmed using PRINSEQ[31], in which both ends with quality scores lower than Q20 were trimmed. The potential nucleotide differences were also validated with BRESEQ version 0.28 [3]. For structural variations, we used the same method as above [5] to detect large ins/del.
Only test data not used for the fitting are plotted. The error bars in the y-axis were obtained by predicted growth rate calculated from 10,000 different sets of test data and training data.
S4 Figure. Expression changes of (a) tbpA and (b) ydiH in the resistant strains.  Table. Summary of relative growth rates of the resistant strains in stress conditions.

S2 Table.
Transcriptome changes observed in the resistant strains.

S3 Table.
List of mutations identified in the resistant strains. Supporting Text S1: Detailed description of common mutations identified in resistant strains.
The number of mutations identified in the resistant strains is shown in Fig. 4, and detailed information about the mutations is presented in Supplementary Table S2. Below, for each stress we discuss the relationship between resistance acquisition, genome mutations and gene expression changes.

Sodium chloride (NaCl) resistant strains
See Main Text.

Potassium chloride (KCl) resistant strains
In KCl resistant strains, in addition to mutations in proU operon as described in Main Text, we found three resistant strains have mutations in the coding region of ftsI, which is involved in the cell division process. Since ftsI is an essential gene, we did not evaluate the effect of the identified ftsI mutations by introducing them into the parent strain. The role of ftsI mutations in KCl resistance remains unclear.

Cobalt chloride (CoCl 2 ) resistant strains
All resistant strains to CoCl 2 stress had mutations in corA 1 , which encodes a transporter mediating influx of Mg 2+ , Ni 2+ , and Co 2+ . The identified mutations included frame-shift and deletion of ORF, suggesting that the disruption of corA activity contributed to CoCl 2 resistance. This hypothesis was supported by the fact that the introduction of the corA mutation found in a CoCl 2 resistant strain into the parent strain significantly increased the growth under CoCl 2 stress (Fig. 5). Four out of five CoCl 2 resistant strains also had mutations in feoA/feoB encoding ferrous iron transporters 2 . Even though the role of FeoA/FeoB in cobalt uptake has not been demonstrated 3 , these identified feoA/feoB mutations suggested that FeoA/FeoB transporters are involved in Co 2+ transport.

Sodium carbonate (Na 2 CO 3 ) resistant strains
Two out of five Na 2 CO 3 resistant strains had mutations in sapABCDF operon, which encode Sap (Sensitive to antimicrobial peptides) ABC importer 4 . To evaluate the effect of the mutations, we introduced sapA mutation in its coding region (found in Na 2 CO 3 -2 resistant strain) into the parent strain, and confirmed that it significantly increased the growth under Na 2 CO 3 stress. Sap transporter is known to contribute alkali stress resistance in Sinorhizobium meliloti, root nodule bacteria 5 , and thus our result might suggest that Sap transporter is involved in alkali stress resistance also in E. coli.

L-Lactate (Lac) resistant strains
Four Lac resistant strains (Lac-1, 2, 3, and 4) had mutations either purT or purU genes, which are involved in purine biosynthesis 6 . In Lactobacillus lactis, it was suggested that perturbing purine biosynthesis is related to multi-stress resistance by changing 3',5'-bispyrophosphate (ppGpp) concentration 7 , which is known as a master regulator of the stringent response to amino acid starvation. Thus, the identified mutations in pur genes might be related to the control of the stringent response, although the details still remain unclear. We confirmed that the introduction of the mutation in the purT coding region significantly increased the growth under Lac stress (Fig. 5).
The expression levels of yojL in these four strains significantly increased (Fig. S7), probably due to the mutations in the promoter region. yojL encodes a periplasmic lipoprotein which is involved in thiamine biosynthesis 8 , and also it is thought to play a role in assembly or maintenance of iron-sulfur clusters 9 . These mutations commonly fixed in yojL promoter region suggests their contribution to Malate resistance by up-regulating its expression level. However, we could not observe a significant growth rate increase under Mal stress when one of these mutations was introduced in the parent strain (Fig. 5).

Methacrylate (MCL) resistant strains
All resistant strains to MCL had mutations in pykF, which encodes a pyruvate kinase that catalyzes the conversion of phosphoenolpyruvate (PEP) into pyruvate in the central metabolic pathway. All these mutations in pykF were nonsynonymous single nucleotide substitutions, which might suggest that the change of enzymatic activity contribute to the methacrylate resistance. To verify this, we introduced the mutation found in a MCL-3 resistant strain into the parent strain, and confirmed that it significantly increased the growth rate under methacrylate stress (Fig. 5). Although the mechanism how pykF mutations contribute to MCL resistance remains unclear, controlling intracellular accumulation of pyruvate caused by MCL stress might be related to the mechanism of resistance.
MCL is known to inhibit the enzymatic activity of pyruvate formate-lyase (Pfl) that converts pyruvate into acetyl-CoA and formate 10 , and thus MCL stress can cause an increase in intra-cellular pyruvate. Since the PEP/pyruvate concentration ratio is an important parameter to control various metabolic pathways, including glucose uptake via sugar phosphotransferase system (PTS) 11 , the decrease of PykF activity by the mutations might play a role in rebalancing the disrupted PEP/pyruvate ratio by MCL stress.

Crotonate (Cro) resistant strains
No common mutation was identified in Cro resistant strains.

Methylglyoxal (MG) resistant strains
See Main Text.

n-butanol (BuOH) resistant strains
Three BuOH resistant strains (BuOH-3, 4, and 5) had mutations in the coding region of cspC, and we confirmed that one of these mutations can slightly increase the BuOH stress resistance (Fig. 5).
CspC is a constitutively produced member of the CspA family of RNA-binding proteins, which increases RpoS expression level presumably by stabilizing its mRNA 12 . RpoS is a central regulator of the general stress response whose up-regulation leads to growth arrest 12 . We found that the rpoS mRNA expression increased in response to BuOH stress addition, while it relaxed to the original levels in BuOH resistant strains (Fig. S8). This result suggests that the growth rate increase under BuOH stress can be partially explained by the down-regulation of rpoS caused by the cspC mutations. Furthermore, all BuOH resistant strains had deletion of region including yneK, ydeA, and marC genes. Interestingly, disruption of marC was also identified in isobutanol resistant E. coli strains obtained by laboratory evolution 13,14 . Although marC is a poorly characterized gene whose function is unknown, these results suggest a positive effect of marC disruption on butanol resistance.

Cetylpyridinium chloride (CPC) resistant strains
Three out of five CPC resistant strains (CPC-2, 4, and 5) had mutations in dcm, which encodes a DNA cytosine methyltransferase. To verify the effect of dcm mutations, we introduced the nonsynonymous single nucleotide substitutions found in CPC-5 strain into the parent strain. This mutant strain exhibited a slight increase of growth rate under the CPC stress (Fig. 5), although the difference was not statistically significant. This result might suggest that the mutation in dcm has a weak fitness gain effect in the CPC stress environments. The mechanism for CPC resistance by these dcm mutations remains unclear. A recent study demonstrated that the deletion of dcm contributes to antibiotic resistance through up-regulation of sugE encoding multidrug efflux transporters 15 . In three of five CPC resistant strains (CPC-3, 4, and 5), mutations are found in the promoter region of sugE, which suggests that these resistant strains acquired CPC resistance by activating the SugE efflux transporter. The fact that the expression levels of sugE significantly increased in CPC-3, 4, and 5 strains and weakly increased in CPC-2 strain supported this hypothesis (Fig. S9).