Evolutionary dynamics determines adaptation to inactivation of an essential gene

Genetic inactivation of essential genes creates an evolutionary scenario distinct from escape from drug inhibition, but the mechanisms of microbe adaptations in such cases remain unknown. Here we inactivate E. coli dihydrofolate reductase (DHFR) by introducing D27G,N,F chromosomal mutations in a key catalytic residue with subsequent adaptation by serial dilutions. The partial reversal G27->C occurred in three evolutionary trajectories. Conversely, in one trajectory for D27G and in all trajectories for D27F,N strains adapted to grow at very low supplement folAmix concentrations but did not escape entirely from supplement auxotrophy. Major global shifts in metabolome and proteome occurred upon DHFR inactivation, which were partially reversed in adapted strains. Loss of function mutations in two genes, thyA and deoB, ensured adaptation to low folAmix by rerouting the 2-Deoxy-D-ribose-phosphate metabolism from glycolysis towards synthesis of dTMP. Multiple evolutionary pathways of adaptation to low folAmix converge to highly accessible yet suboptimal fitness peak.


Introduction
The ability of microorganisms to adapt to different environments is crucial for their long-term survival and particularly important in the development of antibiotic resistance. Bacterial metabolism is often equipped with alternative pathways specialized in the utilization of diverse nutrients, which provides microbes an exceptional ability to thrive in various environments and under stress from antibiotic treatment. In addition, multi-layers of regulation that global control of cellular functions ensures both robustness, against fast fluctuating conditions, and flexibility, through slower transcriptional responses to alternating nutrients (Bennett et al., 2008;Buescher et al., 2012;Millard et al., 2017). When important cellular functions are inactivated, e.g. by genetic mutations, a severe disruption of cellular networks can occur, which poses a major adaptive challenge to cells. Recent studies investigated evolution of E. coli upon inactivation of nonessential enzymes of carbon metabolism that lead to re-wiring through less efficient pathways (Krusemann et al., 2018;Long et al., 2018;McCloskey et al., 2018a, b, c, d). These studies highlight a crucial interplay between regulatory responses and imbalances in metabolite concentrations resulting from gene knockouts, which can be subsequently corrected by mutations elsewhere. Importantly all these studies involved gene-knockout so that adaptation by reversal to WT genotype was not possible.
Less is known about the adaptation mechanisms that follow inactivation of unique cellular processes that are deemed indispensable in microbes. The genes that perform such functions, often classified as essential, are believed to face higher selective pressure and thus evolve slower (Luo et al., 2015). It is thus expected that disruption of essential genes can create a major barrier to the adaptability of bacteria. Nevertheless, a comprehensive study involving knock outs of >1000 genes classified as essential in yeast has shown that a small percentage of mutants could recover viability after laboratory evolution (Liu et al., 2015), revealing that essentiality can, in some cases, be overcome through adaptation. Loss of essential biosynthetic genes have also been observed to occur naturally, under evolutionary conditions where pathway products are provided externally (D'Souza and Kost, 2016), highlighting a context-dependent nature of essentiality. Nevertheless, while adaptation upon inactivation of essential genes in microbes has been demonstrated, its mechanisms remain unknown.
Here we study evolutionary adaptation upon functional inactivation of an essential E. coli enzyme dihydrofolate reductase (DHFR). Past efforts to link fitness effects of chromosomal variation in the folA locus encoding DHFR and their biophysical effects on DHFR protein allowed us to develop an accurate quantitative biophysical model of DHFR fitness landscape (Bershtein et al., 2015a;Bershtein et al., 2013;Bershtein et al., 2012;Bershtein et al., 2015b;Rodrigues et al., 2016). The availability of a clear genotype-phenotype relationship makes DHFR a unique model to study the dynamics and outcomes of evolution to recover from its inactivation by point mutations. In contrast to previous studies that used gene knockouts here we inactivate the chromosomal E. coli DHFR by introducing mutations in the folA locus at a key catalytic residue in position 27, generating strains that express inactive DHFR protein but are viable only with external metabolic compensation, allowing the mutants to adapt to lack of DHFR function by decreasing supplement concentration. The advantage of this approach is that it presents cells with an obvious evolutionary "solution" of reverting the mutant back to WT-variant without massive rewiring that could lead to potentially lesser fit variants. However, an actual outcome may depend on evolutionary dynamics which could revert to WT or other form of active DHFR (higher fitness peak) or converge to a potentially more accessible solution of rewiring to a less efficient metabolic pathways that do not require DHFR function. Therefore, this setup allows us to assess relative roles of height and accessibility of fitness peaks in determining the outcome of evolutionary dynamics.
We show that, depending on starting DHFR variant, partial reversion of DHFR phenotype may indeed occur. However, adaptation to low concentration of external metabolites through metabolic rewiring is the prevalent evolutionary solution due to the availability of a greater number of trajectories leading to consecutive gene inactivation events in two key loci, thyA and deoB. Using omics analysis, we observe global perturbations in metabolites and proteins levels occurring due to DHFR inactivation and upon adaptation, highlighting a key role of regulatory circuits in directing evolution. Finally, we show how adaptation to loss of drug target generates highly resistant strains, and that one important evolutionary solution found here, inactivation of thyA gene, can be also found in clinical isolates of resistant strains of S. aureus (Chatterjee et al., 2008;Kriegeskorte et al., 2014) and H. influenza (Rodriguez-Arce et al., 2017).

Results
DHFR is encoded by the folA gene and is essential for the biosynthesis of purine, pyrimidine and glycine. We sought to inactivate folA in E. coli by introducing mutations in the key catalytic residue D27 ( Figure 1A). - However, such approach poses a methodological challenge since, by definition, the knockout of essential genes results in lethality or infertility. However, folA knockout mutants are conditionally lethal if grown in a culture medium supplemented with a mixture of metabolites comprising thymine, glycine, methionine, inosine and adenine (folAmix) (Howell et al., 1988;Kwon et al., 2010). We envisaged that inactivated DHFR mutants growing in the presence of folAmix could be progressively challenged to adapt to ever decreasing concentrations of this supplement in the growth medium ( Figure 1B). The next step was to develop a protocol that allowed controllable decrease of supplement concentration in growth medium in response to the fitness status of the cells. This was achieved by implementing a fully automated serial passaging scheme as depicted in Figure 1C-D using a Tecan liquid handling robot. In this setup, the growth rate of replicate cultures is monitored by periodic OD readings, and the concentration of folAmix is adjusted downward in each serial passage when cultures exceed a defined growth rate threshold. This feedback control loop ensures mutant strains are continuously challenged to grow at sub-optimal conditions, sustaining a strong selective pressure on the loss of DHFR function. This approach combines the medium-throughput capabilities of plate-based serial dilution methods (currently, up to 32 independent replicate trajectories can be evolved in parallel), and both real-time monitoring of fitness status and control of growth conditions featured in morbidostat setups (Toprak et al., 2012;Toprak et al., 2013), which allows cells to continuously experience exponential growth and sustained selection pressure.

D27 mutations confer growth defects
We selected single or double nucleotide base mutations to replace a key catalytic residue Asp 27 either for Asn, Gly or Phe; the presence of a carboxylate side chain in position 27 is strictly conserved among DHFR orthologues. Purified D27 mutant proteins were characterized (Table 1)  1.3 a Kinetic properties for dihydrofolate reductase catalytic activity (kcat = enzymatic turnover number, KM = Michaelis-Menten constant, kcat/KM = catalytic efficiency, Ki = inhibition constant for trimethoprim) and protein stability properties (ΔTm = difference in melting temperature of folding with respect to wild type, bis-ANS = relative fluorescence from binding of bis-ANS to molten-globule intermediates with respect to wild type protein.) b from Ref (Tian et al., 2015b) and the catalytic activity measurements confirm the lack of significant DHFR function in these variants; catalytic efficiencies (kcat/KM) are several orders of magnitude lower than WT. However, thermal denaturation data obtained for the D27 mutant proteins indicates that their stability is mostly unaffected (or even significantly increased in the case of D27F mutant (Tian et al., 2015b)) showing that, despite the lack of catalytic activity, these proteins retain the ability to fully fold. Importantly, this provides a pathway for restoration of DHFR function over the course of evolution, either by revertant mutations at the D27 locus or, potentially, compensatory mutations elsewhere in the protein. D27 mutant strains were generated by lambda red recombination (Bershtein et al., 2012) and plated in folAmix-containing agar media, however, growth was only visible after 48h and the colonies formed were miniscule. When DHFR function was assayed in cell lysates of D27 mutant strains, no DHFR activity could be detected, confirming that these mutations inactivate DHFR. As expected, growth experiments showed that D27 mutants grow only at high concentrations of folAmix ( Figure 2). Interestingly, the D27G mutant shows slightly better growth at lower concentrations of folAmix than D27F and D27N. Given that the catalytic efficiency of D27G and D27N mutants is comparable, and extremely low, it is very likely that the difference in growth originates from a potential acquisition of a slightly advantageous mutation elsewhere upon genetic manipulation to introduce D27G mutation in the folA locus (see Discussion). We tested D27 mutants for growth in the presence of individual components of folAmix and their different combinations and found that only thymine was essential, although growth with thymine alone is slower compared to growth with all folAmix components (Supplementary Figure S1). The growth rates measured for all D27 mutants at high folAmix fall in the range 70-80% of WT and lag times were typically 1.5-2 fold longer (Figure 2 B). Overall, weaker growth and small colony phenotype show that D27 mutants are severely compromised even in presence of folAmix. Using the evolutionary scheme described previously, we allowed six replicates of each mutant strain to evolve in parallel for about 50 days.

WT D27F
A B

Evolution of D27G mutant
Cultures of D27G mutant strains were first grown in the presence of folAmix diluted 0.1 fold with respect to the initially defined folAmix composition. Throughout the course of the experiment the concentration of supplement mixture further decreased, as imposed by the feed-back control loop, in response to increasing fitness of the mutant strains ( Figure 3A).  (Rodrigues et al., 2016). The resulting prediction (30% of the wild type value) is in fair agreement with the experimental value measured for a BW25113 strain in which D27C mutation was introduced in the wild type folA gene (D27C reconstituted). D) Comparison of growth rates of wild type, evolved D27G and reconstituted D27C strain. E) Dose-dependent growth inhibition by trimethoprim determined for wild type and reconstituted D27C mutant. Solid lines are fits with logistic equation, from which IC50 was determined to be 0.7 ± 0.1 µg/mL for wild type and 77 ± 6 µg/mL for D27C mutant.

Figure 3-
The change in folA mix concentration over time reflects the effect of adaptation. Three trajectories (1,3 and 6) stand out by reaching the ability to grow in the complete absence of folAmix. Trajectories 2 and 5 reached a point where growth dropped below detection limit, and ultimately could not be recovered, whereas trajectory 4 had adapted to grow at low concentrations of folAmix. We hypothesized that mutations in folA locus could have restored DHFR activity in the cultures where reversion of phenotype was observed. Accordingly, Sanger sequencing analysis of this region revealed that in all these three trajectories residue Gly27 had mutated to cysteine. This finding was surprising because alignment of multiple known DHFR sequences shows that cysteine does not occur naturally in position 27; only aspartate or glutamate are observed in this locus. To assess whether Cys27 mutant is functional, this variant was purified and characterized in vitro, and its catalytic properties are compared with wild type and other mutants in Table 1. Although stability-wise D27C mutant is similar to wild type protein, it shows much weaker catalytic efficiency, mostly in terms of retaining only about 0.1% kcat/KM of wild type. To assess if this residual catalytic function is sufficient to explain the reversion of phenotype in the evolved strain we first predicted fitness based on a previously established model (Rodrigues et al., 2016). Taking as input the vitro biophysical properties of purified D27C mutant, namely kcat/KM and bis-ANS fluorescence values, the model predicts a growth rate of about 30% of the wild type strain. To check the validity of this prediction, the mutation D27C was reconstituted in the wild type background by lambda red recombination and cells were plated in folAmix-containing agar plates to lift any selective pressure for DHFR activity. This strain was able to grow in the absence of folAmix supplementation, and the measured growth rate was 85% of the wild type, which is in fair agreement with the in vitro-based prediction which does not take into account a possibility of upregulation in response to DHFR deficiency (Bershtein et al., 2015a) (see Supporting Information for discussion) ( Figure 3B-D). This analysis shows that D27C mutation alone explains most fitness recovery of the evolved strains. Of note is also the fact that significant loss in catalytic efficiency in D27C mutant pays-off with a dramatic 4 orders of magnitude increase in the inhibition constant for trimethoprim. ( Figure 3E) This is in line with previous observations that active site mutations compromise catalytic function but also disrupt pocket interactions with the drug, resulting in high levels of resistance (Bader et al., 2006;Rodrigues et al., 2016). Evolutionary trajectories of D27N and D27F mutants represented in Figure 4 A-B show a marked decrease in folAmix necessary to sustain growth, reaching concentrations nearly three orders of magnitude lower than initially required for naïve strains. Adaptation to grow at low folAmix concentrations was verified by showing that growth curves of individual colonies isolated from evolved strains are similar to the ones obtained in evolutionary dynamics ( Figure 4C). However, these strains cannot grow in the complete absence of folAmix. This result indicates that, unlike some trajectories in D27G evolution, adaptation of D27N and D27F strains did not involve reversion of the folAphenotype. Accordingly, no mutations were found in folA locus in evolved D27N and D27F strains. Other mechanisms must therefore be responsible for the ability to grow at low folAmix concentrations. The growth rate of the evolved strain measured at high concentrations of folAmix reached values comparable to wild type, showing a clear fitness increase with respect to naïve strain. However, at lower concentrations of folAmix the growth rate of evolved strain becomes markedly reduced to nearly half of the wild type value and appears to plateau in an intermediate range to folAmix concentration. We observed a similar plateau when growth was measured at different concentrations of thymine alone (Supplementary Figure S1), from which we concluded that at low folAmix concentrations only thymine is being essential to sustain growth whereas other components are probably too diluted to have a positive impact in fitness. We also evaluated sensitivity to trimethoprim for both naïve and evolved D27F strains. Not surprisingly, in the absence of a functional DHFR and in the presence of folAmix D27F mutant strains are extremely resistant to trimethoprim (Supplementary Figure S2A). We noted, however, that while no decrease in growth was evident in naïve D27F strain up to 1000 µg/mL, the growth rate of the evolved strain dropped abruptly above 500 µg/mL trimethoprim, suggesting that the drug is acting on an unknown target in the cell which is essential in the evolved but not the parent strain. We then tested how the resistance of evolved D27F strain would compare with wild type at very low concentrations folAmix (Supplementary Figure S2B). We found that in those conditions IC50 for the wild type is similar to that measured in the absence of supplement, yet the evolved D27F still shows an extremely high resistance to trimethoprim (IC50=656±78 µg/mL).

D27 mutation causes major metabolic and proteomic changes.
The previous observation that D27 mutant strains show severe growth defects suggests that DHFR inactivation imposes considerable homeostatic imbalance. We focused on the strain D27F to carry out a detailed characterization of the systems-level effects of DHFR inactivation. To that end we carried out high throughput proteomics analysis of D27F mutant strains using LC/MS TMT approach as described earlier (Bershtein et al., 2015a). The method based on differential labeling provides abundances of proteins in the proteome relative to a reference strain which in our case was WT (Bershtein et al., 2015a). We computed Z-scores of the log of relative (to wild type as reference) abundance according to the following equation: , eq. 1 and are the gene i abundances obtained for the mutant and reference (wild type) strains, respectively, denotes a quantity averaged over all genes for a given strain, and is the standard deviation of . Then, we performed a Student t-test to determine which groups of genes had statistically significant variation of protein levels in naïve D27F strain, with respect to wild type, following the functional and regulatory classification of genes into groups by Sangurdekar et al (Sangurdekar et al., 2011). Numerous processes were significantly altered, reflecting broad genome-wide effects of DHFR inactivation (Supplementary Figure S3). Several processes were downregulated in naïve D27F strain, from energy metabolism (aerobic respiration, TCA cycle), to the metabolism of nitrogen, several amino acids, pyrimidines and lipopolysaccharides. On the other end, we found significant upregulation of stress responses, peptidoglycan recycling and salvage of guanine and xanthine. We then performed both targeted and untargeted metabolomic analysis of naïve D27F mutant strain and wild type (see experimental details) to characterize significant changes at the level of metabolites. Likewise, Z-scores for metabolite levels were computed for naïve D27F mutant, with respect to wild type. The metabolites with the highest absolute Z-scores (>1.96) were selected for an enrichment test using MBRole online software (Lopez-Ibanez et al., 2016) to identify pathways in which altered metabolites are overrepresented. The analysis revealed that the metabolism of purines, pyrimidines, beta-alanine, histidine and sulfur were the most significantly changed in naïve D27F mutant comparatively to wild type (Supplementary Figure S4). A detailed scheme representing the changes in metabolites and proteins levels of nucleotide synthesis pathways is shown in Figure 5A.  Not surprisingly, we observe build-up of metabolites upstream the reactions that require reduced folate cofactors (purT, purH, and thyA), whereas metabolites downstream of those reactions are found to be strongly depleted.

Evolved D27 strains partially revert the omics effects of DHFR inactivation
Next, we carried out same LC/MS TMT proteomics analysis of evolved strains (D27F, D27N and D27G) to help identify systems-level changes emerging from adaptation to lack of DHFR function.
To get a glimpse of global proteomics changes we applied PCA analysis which revealed that both wild type and evolved D27G occupy the same quadrant in the space of two principal components, whereas D27 mutants that are folAmix dependent cluster more closely in another quadrant (Supplementary Figure S5). We then computed the Z-scores for the proteomic changes in evolved strains with respect to naïve D27F (Z D27_evo/D27F_naïve ), to assess the effect of evolution on the proteomic levels, and plotted these values against Z D27F_naïve/WT to compare with the initial changes caused by D27F mutation. A clear anticorrelation was observed for all mutants ( Figure 5B), being the strongest for evolved D27G strain. These results indicate that adaptation leads to proteomic changes that generally oppose the immediate effects caused by DHFR inactivation. In the case of D27G, the strong global proteomic shift towards wild type levels is somewhat expected since in the evolved strain, the DHFR activity is partly restored by G27->C mutation. To better characterize the proteomic changes that were specific to evolution of folAmix dependence we performed a functional and regulatory classification analysis to identify which groups of genes were significantly altered in both evolved D27F and D27N, with respect to naïve D27F strains. We found an upregulation of genes involved in glycolysis, fermentation, sugar alcohol degradation and response to low pH, suggesting a metabolic shift towards mixed acid fermentation as a result of adaptation. On the other hand, the processes of DNA restriction/methylation and D-ribose uptake were significantly down regulated with respect to naïve D27F in both evolved D27F and D27N strains. Next, we focused on evolved D27F strain to perform a detailed metabolomic analysis and characterize the changes occurring at the level of metabolites. We computed Z-scores for metabolite changes with respect to naïve D27F (Z D27_evo/D27F_naïve ) and we found significant anticorrelation with Z-scores obtained for naïve D27F (Z D27F_naïve/WT ) ( Figure 5C), indicating that metabolite concentrations in evolved strain generally change to partially recover wild type levels, as observed with proteomics. We found that the pathways significantly enriched in metabolites with the highest Z D27_evo/D27F_naïve scores (>1.96) in evolved D27F were the same as those that had the most significant changes in naïve D27F, with respect to wild type (Supplementary Figure S6). Overall these results show that adaptation to low folAmix concentrations is accompanied by an overall shift in the abundance of metabolites and proteins to partially revert the system-level changes that are caused by DHFR inactivation.

Regulatory responses are altered in evolved D27F strain
Metabolites of purine and pyrimidine biosynthetic pathways were among the most strongly depleted in naïve D27F strain and showed highest increase upon evolution (Supplementary Figure  S4). We reasoned that these depleted metabolites could be limiting growth of the D27 mutants, and that supplementing the culture medium with these metabolites would improve growth. Surprisingly, supplementation with purines strongly inhibited growth of the naïve D27F strain, whereas pyrimidines addition had a stimulatory effect ( Figure 5D, and Supplementary Figure S7). The detrimental effect of purine supplementation suggests that drop in purine metabolites in naïve D27F was a regulatory response to high levels of the purines supplied by the culture medium. In the case of evolved D27F strain, the effect of supplementation of individual nucleotides on growth was comparably weaker or non-existent, indicating that the inhibitory/stimulatory effects were relaxed upon adaptation. It is possible that the profound metabolic and proteomic changes upon DHFR inactivation amount to cells switching their proteomes and metabolomes to a novel stable state. In these conditions, the "programmed" responses that are triggered by high external concentrations of metabolites might actually worsen the fitness of the cell, instead of providing the beneficial advantage that they have evolved for. On the other hand, after evolution, even though metabolites and proteins levels are still quite distinct from wild type values, the regulatory network appears to be better adapted to that new metabolic state. Adaptation can either result from changes in the regulatory network itself, due to mutations in regulatory genes, or from re-wiring of metabolic pathways that establish a different metabolic state for which existing regulatory networks provide more optimal responses.

Whole genome sequencing
Whole genome sequence (WGS) analysis was performed to identify the genetic basis for the systems-wide changes observed in evolved strains. In addition, to characterize the dynamics of mutation fixation, individual colonies from intermediate evolutionary time points of a selected D27F trajectory were also sequenced.  Figure S6). However, deoB inactivation saves 2-deoxy-D-ribose-1-phospahte that is required in dTMP production, instead of being utilized for energy production in glycolysis, which is significantly upregulated in evolved D27F, with respect to naïve D27F (average Z-score= 0.593, P= 0.015 Figure 6 A-B summarizes the WGS results obtained for evolved and naïve strains from each D27 mutant. We found that D27N and D27G strains, but not D27F, had additional background mutations that must have been acquired prior to starting the evolution experiments; although all D27 strains were constructed from the same parent wild type strain, we could not prevent the appearance of mutations at any given stage prior to the evolution experiments. It is possible that these mutations may provide some fitness advantage with respect to "pure" D27F strain (see below and discussion). Comparison between evolved and naïve D27G strains reveals that no other mutation was fixed upon evolution besides the one nucleotide change at the D27 locus of the folA gene, corresponding to the Gly->Cys substitution described earlier. This is fully consistent with the earlier conclusion that D27C mutation alone explains the fitness recovery observed in D27G evolution. From the analysis of evolved D27F and D27N we can identify two common hotspots for mutations, the genes thyA and deoB, encoding thymidylate synthase and phosphopentomutase, respectively. For that reason, the sequence of these two loci were determined by Sanger analysis for all other trajectories that evolved to grow at low folAmix concentrations, including trajectory 4 in D27G evolution. Strikingly, in a total of 12 trajectories all were found to have mutations in thyA and 11 had deoB mutated ( Figure 6C). We noted as well that evolved D27F strains also carried mutations in other genes, crr (7bp deletion), rpe and yhdH (4bp deletion), however, these were not observed among other trajectories indicating that these are most likely either random passenger mutations or provide only marginal advantage on a specific genetic background. On the other hand, the occurrence of deletions in both thyA and deoB strongly implies that functional inactivation of these gene products arises as a main mechanism of adaptation to the lack of DHFR function. We reasoned that complementation with either thymidylate synthase or phosphopentomutase should create a significant fitness disadvantage to the evolved cells, which can be reverted if these enzymes become inactivated. To verify this prediction, we transformed evolved D27F strains with plasmids expressing either wild type thyA or deoB and compared the ensuing growth rates of these strains with transformants carrying the same plasmids expressing inactivated thyA and deoB, respectively. Expression of functional ThyA and especially DeoB were found to be highly deleterious at low folAmix concentrations, as shown Figure 6 D-E, whereas no significant change in growth was observed upon expressing inactive mutants. The gene deoB is part of an operon responsible for deoxyribose degradation, that includes deoA, deoC, and deoD, which is under the control of 2-deoxyribose-5-phosphate-inducible deoR repressor (Hammer-Jespersen and Munch-Ptersen, 1975). DeoB catalyzes the interconversion of 2-Deoxy-D-ribose-1-phosphate and 2-Deoxy-D-ribose-5-phosphate, and, while the former is necessary for synthesis of dTMP from thymine by DeoA, the latter can be further degraded in glycolysis, through conversion into acetaldehyde and D-glyceraldehyde 3-phosphate by DeoC ( Figure 6F). Repression of deo operon might be deleterious because expression of DeoA is essential for thymine uptake, however coexpressing high levels of DeoB and DeoC can also have a negative effect by directing 2-deoxy-Dribose-1-phosphate for degradation. Therefore, inactivating deoB or deoC is an "economy" solution preventing key metabolite in thymine uptake to be wasted in energy production, allowing cells to efficiently use small amounts of thymine available from media supplement for nucleotide synthesis (Dale and Greenberg, 1972). It is also clear that, since we found no mutations in regulatory genes, the origin of the attenuation of the nucleotide-mediated inhibitory/stimulatory effects observed in evolved D27F strain are most likely due to global changes in protein/metabolite concentrations.

Figure 6 -Loss of function mutations in two genes, thyA and deoB, lead to adaptation to low folAmix concentrations. A) Mutations identified in naïve and evolved D27 mutant strains. Vertical lines identified with NGS represent mutations identified by whole genome sequencing whereas the remaining cases were identified by Sanger sequencing. Details are presented for trajectory 2 of D27F evolution and include the mutations identified at various passages and the folAmix profile obtained for that trajectory. B) Mutations identified in the most relevant genes. C)
Overall, these results show that two key loss of function mutations alone are responsible for the partial adaptation of D27F to loss of DHFR activity and result in significant metabolic, proteomic and regulatory shifts observed in the evolved strains.

Discussion
In this study we set to explore the evolutionary mechanisms that follow inactivation of the essential gene that encodes DHFR. Such scenario is akin to antibiotic-mediated target inhibition, but distinct in several aspects. When treated with DHFR inhibitor trimethoprim, bacterial cells recurrently acquire resistance by mutations that decrease drug binding affinity and/or enhance target abundance levels (Toprak et al., 2012). In contrast, genetic inactivation of DHFR at the conditions studied here constrains evolution towards solutions that appear to be mutually exclusive, either restoration of DHFR catalytic activity or adaptation to lack of DHFR function. Interestingly, both evolutionary outcomes were observed. While D27F, D27N and one trajectory in D27G convergently adapted to lack of DHFR function, other trajectories in D27G mutant reached reversion from thymine auxothrophy phenotype. It is important to address the possible reasons for the observed outcome. A major confounding factor could be the presence of background mutations in some of these strains. We note the potential role of the mutation in upp gene present in D27G strain. This gene encodes the enzyme uracil phosphoribosyltransferase that participates in pyrimidine salvage pathway and, quite strikingly, appeared consistently depleted in proteomics analysis of all D27 mutants. The slightly higher fitness previously noted for D27G (Figure 2A) could be due some beneficial effect of the upp and/or ymfD mutation. Conceivably, such small difference in fitness could be a key factor determining the outcome of evolution, but further experiments are needed to test this hypothesis. The codon structure of each mutant studied here and bias in the frequency of each mutation type can also affect the likelihood that phenotypereverting mutations may arise. We note that the nearest accessible high-fitness Asp/Glu codon for D27F mutant is two nucleotide substitutions away, whereas both D27G and D27N could be rescued by a single G->A or A->G transition, respectively. Only D27G was found to revert, however, not to wild type codon, but to a less catalytically efficient cysteine residue caused a G->T transversion. In this respect, it is rather surprising that D27C mutation was fixed in 3 independent trajectories when transversions are generally regarded to be less likely than transitions. Another important factor in determining the fate of D27 mutant evolution appears to be the early fixation of thyA inactivating mutations. Thymidylate synthase is the only known enzyme in E. coli able to synthetize dTMP de novo from dUMP, and inactivation of thyA inevitably commits the cell to thymine auxothropy; reversion of DHFR function on the thyAbackground is not expected to change this phenotype, as it is known that folA + thyAstrains are thyminedependent (Bertino and Stacey, 1966). Competition between DHFR reversion and thyA inactivation thus appears to be decisive in determining the evolutionary solution. However, since inactivation of a gene can be achieved by multiple ways, there is an incomparably greater number of accessible evolutionary trajectories leading to thyA knockout when compared with the very limited available mutational routes that restore DHFR function.
Upon thyA inactivation, the only available route for dTMP synthesis involves DeoA-catalyzed formation of thymidine from thymine and 2-deoxy-D-ribose-1-phosphate ( Figure 6F). However, DeoA is part of the deo operon that is induced by high levels of 2-deoxy-D-ribose-5-phosphate. This regulatory scheme allows cells to recycle excess metabolites into energy production. However, on the background of thyA inactivation, the co-expression of deo operon proteins creates simultaneously a solution and a problem for the uptake of thymine. Now DeoA becomes an essential protein, but its function is opposed by the roles of DeoC and DeoB which divert 2-Deoxy-D-ribose-1-phosphate into degradation via the glycolysis pathway. The evolutionary solution found in D27 strains converged to the inactivation of a single gene of the operon, deoB, an event that is sufficient to block the degradation pathway and yet the regulatory structure is retained. This critical example illustrates the role of regulatory circuits in constraining evolution when disruption of an essential gene forces the system to adopt a metabolic state well beyond normal-operating homeostatic boundaries. Understanding these processes is critical for guidance towards new strategies for antibiotic resistance prevention. Perhaps not surprisingly, we found that pathways of adaptation to genetic inactivation of an antibiotic target provide a route to high levels of resistance. Strikingly, mutations in thyA are also known to confer trimethoprim resistance in bacterial clinical isolates of S. aureus (Chatterjee et al., 2008;Kriegeskorte et al., 2014) and H. influenza (Rodriguez-Arce et al., 2017), which emphasizes the need of laboratory studies of microbial adaptation as useful models to tackle serious global threat. Sangurdekar, D.P., Zhang, Z., and Khodursky, A.B. (2011). The association of DNA damage response and nucleotide level modulation with the antibacterial mechanism of the anti-folate drug trimethoprim. BMC Genomics 12, 583. Scheltema, R.A., Jankevics, A., Jansen, R.C., Swertz, M.A., and Breitling, R. (2011). PeakML/mzMatch: a file format, Java library, R library, and tool-chain for mass spectrometry data analysis.

Contact for Reagent and Resource Sharing
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Eugene I. Shakhnovich (shakhnovich@chemistry.harvard.edu).

Construction of D27 mutants
Inactive D27 mutants were created by lambda-red recombination (Datsenko and Wanner, 2000) according to a previously described procedure (Bershtein et al., 2012) with some modifications. Briefly, a pKD13 plasmid was modified to contain the entire regulatory and coding sequence of folA gene, flanked by two different antibiotic markers (genes encoding kanamycin (kanR) and chloramphenicol (cmR) resistances) and approximately 1kb homologous region of both upstream and downstream chromosomal genes flanking folA gene (kefC and apaH, respectively). The entire cassete was amplified and transformed into BW25113 cells with induced Red helper plasmid pKD46, and cells were recovered in SOC medium containing 1x folAmix (adenine 20 ug/mL, inosine 80 ug/mL, thymine 200 ug/mL, methionine 20 ug/mL and glycine 20 ug/mL). Transformants were plated in agar media containing both antibiotics and folAmix. Correct integration of the desired mutations was confirmed by Sanger sequencing of folA gene. Plasmid pKD46 was removed by plating cells at 37 °C twice in the absence of antibiotic selection.

Automated experimental evolution
A liquid handling instrument Tecan Freedom Evo 150 equipped with Magellan plate reader and Liconic shaker was used in this work. The experiments were done at 30 ˚C and using M9 minimal media supplemented with 2g/L glucose, 34µg/mL chloramphenicol and 50µg/mL kanamycin. The general procedure involves up to four 96-well plates that are used for serial dilutions of bacterial cultures. Each working plate can carry 8 trajectories that are positioned in a single column. In this work we used the first and eight wells in the column with media alone to control for contamination. The experiment starts by placing the 200 µL cultures/control in the first column of the 96-well plate. All the four plates are incubated in a shaker at 30 ˚C and at every 30 min the OD of each plate is determined alternately. The growth rate of each culture is calculated from OD measurements over time. To ensure proper comparison, the growth rate was computed only from OD values within a specified range (0.1-0.25). When the average OD of the six experimental replicates in a plate exceeds a threshold of 0.30, each culture is diluted into the next adjacent column by mixing a calculated volume of both culture and fresh medium and folAmix so that the initial OD is 0.01 in a total of 200 µL. At this point, the remaining portion of the previous culture is mixed with glycerol in an auxiliary plate and subsequent freezed at -80C. The cycle is repeated through the entire experiment. At every passage, the growth rate is compared with a threshold value (0.16 h -1 ) and whenever a culture exceeds this value the concentration of folAmix is halved.

Media and Growth Conditions
Cell cultures were recovered by inoculating fresh M9 media medium supplemented with 0.2% glucose and 1x folAmix with a portion of -80˚C glycerol stocks taken at different passages. After recovery for several hours, the cultures were plated in selective agar media containing 1x folAmix and incubated at 30˚C. Individual colonies were grown in M9-media+folA mix and stored in glycerol at -80˚C for later analysis. Cell cultures grown overnight were diluted in fresh M9 minimal media containing 1x folAmix and antibiotics and were grown for additional 4-6 h. Cultures were then pelleted by centrifugation and washed 3 times in fresh M9 without folAmix. Microplates containing 150 µL of M9 minimal with 0.8g/L glucose, antibiotics, and varying concentrations of folAmix were then inoculated with each culture at a starting OD of 0.0005. Growth measurements were performed in a Infinite® 200 PRO plate reader for 48h at 30˚C with constant shaking. Growth rate values are represented as mean ± SD from at least 3 biological replicates.

Protein purification and characterization
D27 mutant fused to C-terminal (6x) His-tag were overexpressed using pFLAG expression vector under isopropyl β-D-1-thiogalactopyranoside (IPTG) inducible T7 promoter. The recombinant proteins were purified from lysates on Ni-NTA columns (Qiagen) as described previously (Rodrigues et al., 2016). DHFR kinetic parameters were derived from the analysis of progress-curve kinetics of NADPH oxidation in the presence of dihydrofolate using the software Kintek Explorer (Johnson, 2009) as described before (Rodrigues et al., 2016).

Predicting fitness of D27 mutants
The fitness of E. coli DHFR mutants can be predicted from in vitro biophysical parameters using a simple metabolic model described in an earlier work (Rodrigues et al., 2016). The metabolic flux through DHFR reaction in DHFR mutant strains, normalized to wild type, can be approximated to: eq. 2 Where, , and are, respectively, the catalytic efficiency, trimethoprim inhibition constant and intracellular protein abundance of the mutant or of the wild type, and and are, respectively, the concentration of trimethoprim in the culture medium and the ratio of intracellular vs extracellular concentrations of trimethoprim, defined previously to be 0.  (Rodrigues et al., 2016). To connect flux to fitness, we first measure the growth rate of wild type E. coli cells at various concentrations of DHFR inhibitor trimethoprim.
Then fitness is plotted against the flux through DHFR reaction calculated at each concentration of inhibitor using equation 2, and catalytic parameters determined in vitro (Table 1). Protein abundance is determined by measuring the total catalytic activity in cell lysates, as described earlier (Rodrigues et al., 2016). Finally, we use equation 3: eq. 3 to fit the fitness vs V norm data points to determine the parameter B, which represents the normalized flux at which growth rate is reduced by half. From this equation, the in vitro parameters determined for D27 mutants are used in equation 2 and 3 to predict fitness. The protein abundance of D27 mutants, however, cannot be determined from enzymatic assays in cell lysates, as these are inactive variants. Instead, protein abundance was predicted using an inverse relationship with relative bis-ANS fluorescence of those variants (Bershtein et al., 2013;Rodrigues et al., 2016). In qualitative terms, our model accurately predicted that the D27C variant is able to grow in the absence of folAmix supplement. However, the measured fitness was somewhat higher than the value predicted by the calculations. It is possible that some assumptions are not entirely met. One important assumption is that the concentration of DHFR and substrate dihydrofolate are either always constant or change only as a function of the flux calculated by equation S1. The latter implies that, for every given value of flux, the changes in DHFR and dihydrofolate concentration will always be the same, and thus always comparable regardless of what parameters are used as input in equation S1. However, this might not be always true. For example, a tight feedback control regulates the DHFR promoter activity, in response to the metabolic needs of the cell (Bershtein et al., 2015a;Bershtein et al., 2015b). Likewise, a decrease in DHFR catalytic function is also expected to result in a considerable level of substrate build up. In these conditions, the magnitude of substrate KM for each variant might become more relevant. In our calculations, however, substrate saturation is not considered. This may lead to deviations, especially when KM values differ significantly, as we find in this work. How much these dynamic processes may differently affect the flux through DHFR reaction in each mutant is still to be determined.

Effect of deoB and thyA expression
The genes deoB and thyA and corresponding mutants deoB p.309Stop and thyA 755-762del were cloned into a modified pTRC plasmid in which the laqI and promoter region were replaced by the tetR repressor gene and promotor region derived from plasmid Plasmid pEM-Cas9HF1-recA56 (addgene Addgene plasmid # 89962 (Moreb et al., 2017)