Introduction

The dosage of a gene, i.e., relative copy number, directly impacts its level of expression and can alter organismal phenotype and fitness [1–3]. Mutations that increase copy numbers, termed as gene duplications and amplifications (GDAs), are implicated in diverse biological phenomena including cancer progression [4–6], speciation [7–9], reproductive isolation [10–12] and antimicrobial resistance (AMR) [2, 13]. GDAs are important for rapid adaptation since their frequencies of occurrence can be higher by several orders of magnitude than point mutations [2, 13–16]. Following the initial GDA mutation, however, amplified genes could be lost over time, pseudogenized, sub-functionalized or neo-functionalized [17–19]. Which of these possibilities is realized depends on the nature and strength of selection imposed by the environment, as well as the fitness cost associated with maintaining multiple gene copies [17–19].

In bacteria, GDA mutations are frequently encountered during adaptation to antibiotics [13] or heavy metals [20, 21]. They occur by RecA-dependent nonequal recombination as well as RecA-independent processes and often involve repeat sequences such as in IS-elements [22–25]. Amplified genes in drug-resistant bacteria code for effectors of resistance such as drug inactivating enzymes, drug target modifiers or “by-pass” mechanisms [26–28], but may also compensate for resistance-associated fitness costs [29]. More recently, the view that GDAs are less stable than point mutations has emerged. For instance, transient alterations in gene dosage are a common mechanism of heteroresistance, i.e. resistant sub-populations within a drug-sensitive bacterial isolate [20, 30, 31]. Further, it is estimated that as many as 20% of cells in a bacterial population may harbour a duplication of some genomic region [32]. These copy number polymorphisms create within-population heterogeneities in gene expression, serving as a bet-hedging strategy with adaptive value in fluctuating environments [33].

Dihydrofolate reductase (DHFR) is a ubiquitous enzyme crucial for growth and division. It catalyzes the NADPH-dependent reduction of tetrahydrofolate, an important cofactor in nucleotide and amino acid biosynthesis. DHFR is coded by an essential gene in most organisms and its deletion or pharmacological inhibition results in cell death or growth stasis [34]. Not surprisingly, anti-folates are effective therapeutic agents for diverse pathologies [35, 36]. For instance, methotrexate, used in the treatment of cancers and autoimmune conditions [37], pyrimethamine, used as an anti-malarial [38], and trimethoprim, used as an anti-bacterial [39], are all competitive inhibitors of DHFR. Interestingly, copy number amplification and consequent overexpression of DHFR confers resistance in all three systems [40–43]. Thus, the gene dosage of DHFR appears to be tightly regulated across evolutionarily distant organisms, making it a useful system to understand the principles underlying selection and evolutionary fate of GDA mutations.

In Escherichia coli DHFR is coded by folA. Several point mutations that confer trimethoprim resistance have been identified in this enzyme [44–47]. Additionally, mechanisms that transcriptionally up-regulate folA expression are also known in trimethoprim resistant E. coli. These include cis-regulatory mutations in the folA promoter and trans-acting loss of function mutations in mgrB, which upregulate folA transcription in a PhoQP-dependent manner [44–46, 48]. Earlier work from our group identified Lon, a conserved quality control protease in bacteria, as a post-translational regulator of trimethoprim resistance. Lon degraded mutant DHFRs and reduced their half-lives in vivo, directly affecting genotype-phenotype relationships of resistance-conferring mutations [49, 50].

In this study, we first compare early mutational events in wild type and Lon-deficient E. coli challenged with trimethoprim. We show that a large genomic duplication encompassing folA confers trimethoprim resistance in both genetic backgrounds, though its frequency of occurrence is elevated in lon-knockout E. coli. We then use this system to ask what factors determine the fate of duplicated folA gene copies during evolution at different trimethoprim pressures.

Results

A large IS-flanked duplication encompassing folA mediates rapid adaptation to trimethoprim in E. coli

We had previously evolved 5 replicate populations of wild type E. coli K-12 MG1655 in 300 ng/mL of trimethoprim ( ∼0.3 x MIC) for 25 generations. We then isolated and sequenced the genome of a dominant trimethoprim-resistant clone from each population (designated as WTMPR1-5) (Figure 1 A). These sequencing results have been reported earlier [46]. Briefly, all isolates harboured loss-of-function mutations in the mgrB gene (Figure 1 B) that resulted in PhoQP-dependent overexpression of folA [46]. Additionally, isolates WTMPR1, 3 and 5 had mutations at the active site of DHFR, while WTMPR4 harboured a large genomic duplication of ∼0.6 Mb encompassing folA (Figure 1 B, C). We performed similar evolution experiments starting with an isogenic E. coli Δlon strain to ask if early adaptation to trimethoprim was altered by Lon-deficiency (Figure 1 A). Wild type and lon-deficient E. coli have similar MIC for trimethoprim, though the lon-knockout has a higher IC50 [49]. Among trimethoprim-resistant isolates derived from E. coli Δlon, designated as LTMPR1-5, three harboured mutations at the mgrB locus (LTMPR1-3, Figure 1 A, B, Table 1). LTMPR4 and LTMPR5 harboured mutations in pitA, coding for a metal-phosphate symporter, that have been reported earlier in trimethoprim adapted E. coli [51]. LTMPR5 had a deletion encompassing pitA, indicating that loss of PitA was beneficial in trimethoprim. A knockout of pitA indeed showed greater colony formation on trimethoprim-supplemented media demonstrating that mutations in this gene were adaptive (Figure S1), though the underlying molecular mechanisms are unknown at present. Interestingly, none of the LTMPR isolates had mutations in DHFR or its promoter. Instead, LTMPR2-4 had evolved GDA mutations similar to WTMPR4, while LTMPR1 had a duplication in an overlapping but shorter genomic stretch (Figure 1 B, C, Table 1).

Figure 1.

Table 1.

Duplicated genomic regions in WTMPR4 and LTMPR1-4 were flanked by IS-elements (Figure 1 C). Seeking to explain the higher frequency of GDA mutations in the trimethoprim-adapted lon-knockout, we asked whether IS-transposases were sensitive to proteolysis by Lon. Indeed, heterologously expressed InsB and InsL transposases had greater stability in the lon-knockout (Figure S1), suggesting that their accumulation may contribute to higher GDA frequency in E. coli Δlon. It would follow that a higher proportion of GDA mutants would be isolated from E. coli Δlon during adaptation in other environments as well. In agreement, GDA mutations were detected after evolution of E. coli Δlon in 2 of the 5 additional growth conditions tested (Figure 1 D, S2). In all cases duplications were flanked by IS elements (Figure S2). Further, duplicated genomic regions were similar among isolates derived from the same growth condition, but different across environments (Figure 1 B, S2). On the other hand, wild type E. coli did not show GDA mutations in any of the other tested conditions (Figure 1 D).

The duplicated region in trimethoprim-resistant isolates included genes previously associated with resistance to trimethoprim, namely, folA, and acrA and B, which code for components of the AcrAB-TolC efflux pump (Figure 1 C). In one of the isolates, i.e. LTMPR1, the amplified stretch encompassed folA but not acrAB without affecting trimethoprim IC50 values (Figure 1 B, C). We therefore reasoned that duplication of folA was the main effector of adaptation to trimethoprim. Immunoblotting showed that DHFR protein was indeed overproduced in isolates LTMPR1 and WTMPR4 over and above the increase mediated by inactivation of mgrB alone (Figure 2 A).

Figure 2.

In addition to folA, ∼173 and ∼531 other genes were also duplicated in LTMPR1 and WTMPR4, respectively. To understand the global transcriptional burden imposed by this amplification, we analysed the gene expression profiles of WTMPR4 and LTMPR1 using RNA-sequencing (Supplementary File 1). Since both isolates had loss-of-function mutations in mgrB, which cause large changes in global gene expression [51], we also analysed the transcriptomes of E. coli ΔmgrB and E. coli ΔlonΔmgrB strains for comparison (Supplementary File 1). Transcript level of folA was higher in LTMPR1 (2.3-fold) and WTMPR4 (5.7-fold) over respective controls (Supplementary File 1). Other duplicated genes were mildly overexpressed by ∼1.7-fold on average in LTMPR1 (over E. coli ΔlonΔmgrB) and ∼1.3-fold in WTMPR4 (over E. coli ΔmgrB) (Figure 2 B, Supplementary File 1). However, there was large variation in expression level among duplicated genes (Figure 2 B), indicating that increase in gene dosage did not necessarily result in enhanced expression levels. Instead, global gene expression profiles of WTMPR4 and LTMPR1 were correlated with E. coli ΔmgrB and E. coli ΔlonΔmgrB respectively (Figure 2 C, D) and most of the top 50 up/down-regulated genes in WTMPR4 and LTMPR1 could be attributed to loss of mgrB alone (Figure S3). In line with this result, strains with and without genomic amplifications had comparable fitness costs (∼20 %) in drug-free media, similar in magnitude to the cost of mgrB inactivation (Figure 2 E). Further, the fitness cost was alleviated in growth media supplemented with high Mg²⁺, i.e. by inhibiting the PhoQP pathway, confirming the source of this cost to be loss of mgrB rather than the GDA (Figure 2 F). Thus, we concluded that large GDA mutations encompassing folA led to DHFR overexpression and represented an alternate mechanism for trimethoprim-resistance in E. coli with undetectable additional cost.

Maintenance of folA duplication is linked to antibiotic-generated demand for DHFR expression

We next investigated the evolutionary fate of duplicate folA gene copies by evolving the LTMPR1 isolate, which harboured the shortest duplication, in the absence of antibiotic. We also evolved this isolate at sustained drug pressure, i.e. trimethoprim at the same concentration at which it had originally emerged (300 ng/mL). Using three replicate populations in each condition, we asked how DHFR expression and trimethoprim IC50 values changed over ∼252 generations of evolution and used genome sequencing to identify the underlying mutational basis. As controls, we performed similar evolution experiments with WTMPR4 (which harboured a GDA and a functional Lon protease) and WTMPR5 (which harboured a folA point mutation and a functional Lon protease) (Figure 3 A).

Figure 3.

Evolution in trimethoprim-free media

In drug-free conditions, trimethoprim IC50 of LTMPR1 rapidly declined and a concomitant loss of colony formation on trimethoprim-supplemented plates was observed (Figure 3 B, C). Evolved WTMPR4, but not WTMPR5, showed similar loss of colony formation on trimethoprim indicating that resistance mediated by GDA mutations was less stable than point mutations regardless of the presence of Lon (Figure 3 C). For LTMPR1, these trends were explained by a decline in DHFR expression levels over evolution and reversion to a single folA copy (Figure 3 D, E, S4). Mutation in mgrB was, however, consistently maintained till the end of the experiment (Supplementary File 2). In agreement with genomics data, trimethoprim IC50s of evolving populations stabilized at a value similar to E. coli ΔlonΔmgrB (Figure 3 B). Interestingly, a second GDA mutation (“GDA-2”, ∼2.2 Mb to 3.3 Mb) emerged in these populations (Figure 3 E, S4). The GDA-2 mutation was apparently lost and regained multiple times during the experiment (Figure 3 E, S4), suggesting a dynamic maintenance in the population rather than fixation. Further, the region of the genome amplified in GDA-2 was similar to what we had observed when E. coli Δlon had evolved in drug-free media (Figure S2, S4), indicating that GDA-2 was likely a fitness enhancing mutation unrelated to trimethoprim adaptation.

To further investigate the emergence of GDA-2 during evolution of LTMPR1 in drug-free media we sequenced 6 randomly picked colonies from each of the lineages at the final time point (Figure 3 F). All the sequenced isolates had a single copy of folA and a majority of isolates showed the presence of GDA-2 (Figure 3 F, Supplementary File 2). Interestingly, some of the isolates that did not harbour GDA-2 had evolved mutations in rpoS (Figure 3 F, Supplementary File 2). Loss-of-function mutations in RpoS are fixed at late time points during adaptation to trimethoprim and enhance the fitness of E. coli without conferring resistance [46, 51]. We also sequenced a few randomly picked isolates from WTMPR4 and WTMPR5 evolved lines (2 from each of the 3 replicate lineages in each case) to check for the presence of GDA-2 and rpoS mutations. As expected, clones derived from WTMPR4 had reverted to a single copy of folA, while WTMPR5 retained its ancestrally inherited folA mutation. None of them had evolved the GDA-2 mutation, while a few of the isolates from WTMPR4 and WTMPR5 lines had acquired rpoS mutations (Figure 3 G, H).

Consistent with earlier studies, and the role of mutations in rpoS in general fitness enhancement in E. coli, we found that an LTMPR1-derived isolate with mutant rpoS recovered fitness in drug-free media (Figure 4 A). Interestingly, an isolate with the GDA-2 mutation had also recovered fitness (Figure 4 A). Since GDA-2 occurred only in a lon-deficient background, and GDA-2 and rpoS mutations were mutually exclusive, we wondered whether the fitness enhancing effects of GDA-2 were also mediated by reduced RpoS activity. Lon is known to reduce the levels of RpoS during the exponential growth phase of E. coli [52]. Therefore, we reasoned that at least some targets of RpoS would be overexpressed in a lon-knockout. To test this prediction, we analysed the transcriptome of E. coli Δlon and specifically looked at the expression levels of 168 known targets of RpoS (Supplementary File 3). Of these, 41 genes were overexpressed in the lon-knockout by at least 1.5-fold over wild type, including genes induced in the stationary phase and by various stresses such as poxB, aidB and cbpA (Figure 4 B, C). Levels of rpoS mRNA were however unaffected in this strain. Next, we compared the expression levels of these 41 RpoS-targets in LTMPR1 and E. coli Δlon ΔmgrB with E. coli Δlon as the reference and found that most had similar expression levels in these three strains (Figure 4 C). Indeed a few of these genes showed higher expression levels in LTMPR1 and E. coli Δlon ΔmgrB, which was not unexpected since loss of mgrB is known to stabilize RpoS [53]. We checked their expression in isolates derived from LTMPR1 which harboured either GDA-2 or rpoS mutations. In both isolates RpoS target genes were down-regulated indicating a re-setting of RpoS activity (Figure 4 C). Further, in the isolate harbouring GDA-2, rpoS mRNA levels were themselves lowered to 0.49x of E. coli Δlon (Supplementary File 3). Finally, deletion of rpoS from E. coli Δlon enhanced the fitness of the strain (Figure 4 D). Thus, evolution of LTMPR1 in antibiotic-free media was characterized by rapid reversal of GDA-mediated resistance and acquisition of mutations that enhanced fitness by reducing the levels of RpoS-transcribed genes (Figure 4 E).

Figure 4.

Evolution at constant trimethoprim pressure

In the presence of trimethoprim (300 ng/mL), antibiotic resistance level and DHFR expression remained high over the course of evolution of LTMPR1 (Figure 5 A, B, C, D). WTMPR4 and WTMPR5 too remained resistant, and in the case of WTMPR4 an increased colony forming efficiency was observed at a late time point (Figure 5 C). These data suggested that under drug pressure high folA copy number was likely to be maintained. Indeed, population sequencing of LTMPR1 lineages showed that the ancestral GDA mutation was maintained until 196 generations of evolution in all three replicates (Figure 5 E, S5). By 252 generations, however, a majority of bacteria in Lines A and C had reverted to a single folA copy (Figure 5 E, S5) and had acquired a point mutation either in the folA gene (Line C) or its promoter (Lines A) (Figure 5 E, Supplementary File 4). Line B, on the other hand, retained high folA copy number, but with a reduction in the length of the genomic duplication (Figure 5 D, S5). Line B had also acquired a mutation in the folA promoter, however its frequency in the population was low (Figure 5 D, Supplementary File 4).

Figure 5.

In agreement with population sequencing, none of the resistant isolates derived from LTMPR1 Line A at the final time point showed the ancestral GDA. Instead, they all harboured a C_-35T mutation in the folA promoter that is known to up-regulate DHFR expression [51] (Figure 5 F, Supplementary File 4). Among Line B isolates, many had retained the ancestral duplication, without acquiring any additional mutations in folA. One of the isolates from this Line that had reverted to a single folA copy had gained the A_-64T folA promoter mutation (Figure 5 F, Supplementary File 4). In Line C too we found that isolates that had lost the ancestral GDA had evolved a point mutation in the folA promoter (+TG_-30) or a Trp₃₀Arg mutation in DHFR (Figure 5 F, Supplementary File 4). We found similar results in sequenced isolates from WTMPR4 lineages, i.e. loss of GDA and gain of a point mutation either in folA or its promoter (Figure 5 G). WTMPR5-derived isolates did not accumulate additional mutations in folA but did acquire mutations in rpoS (Figure 5 H). Interestingly, two isolates from LTMPR1 Lines B and C had evolved new GDA mutations encompassing folA. In one of the Line B isolates (Isolate 2) the ancestral GDA was replaced by 5x amplification of a shorter genomic region (∼12 kb) encompassing folA (Figure S6). Both ends of this newly amplified region were flanked by IS1 elements (Figure S6). On the other hand, one of the Line C isolates (Isolate 5) had an expanded duplication beyond what was originally found in LTMPR1 (Figure S6). Also, both copies of folA in this isolate coded for DHFR-Trp₃₀Arg (Figure 5 E, S6), suggesting that the original GDA was lost and an extended GDA mutation was gained after fixation of the mutant folA allele.

Based on these results we concluded that even under drug pressure, GDA mutations were reversed in most bacteria. Notably, reversion to a single folA copy was slower than in drug-free media as it first required the fixation of point mutations at the folA locus. However, in a smaller sub-population, more than one folA gene copy was maintained. This maintenance was dynamic and likely involved repeated loss and gain of GDA mutations.

To further investigate the replacement of GDA by point mutants, we checked to see if point mutants could invade a population harbouring GDA mutations under drug pressure. We allowed LTMPR1-derived isolates with point mutations in folA (marked with ΔlacZ::Cat) to compete with a GDA-harbouring isolate from the same time point in evolution and asked how rapidly the former could invade a population of the latter. We performed these competitions with a contemporary GDA-harbouring isolate rather than the ancestor to nullify the effects of other media adaptations such as rpoS mutations that may have been acquired over evolution. Since for this assay we used point mutant strains genetically marked with ΔlacZ::Cat we verified the neutrality of the ΔlacZ::Cat marker by checking if the marked isolate could invade its unmarked parent. In the presence of trimethoprim (300 ng/mL) point mutants were enriched over the GDA mutant in co-culture by about 100-fold (Figure 6 A) recreating the results of our evolution experiment. A possible explanation for this finding may be higher resistance level conferred by folA point mutations than GDA. Indeed, some isolates with folA point mutations had higher IC50 values compared to a contemporary GDA-harbouring isolate (Figure 6 B). However, regardless of IC50 values, all isolates also had similar doubling times in the presence of 300 ng/mL trimethoprim (Figure 6 B). Thus, though there were no large observable differences in growth rate between isolates, point mutants could indeed invade and outcompete GDA mutations, reflecting either subtle differences in fitness between them or their relative stabilities in culture.

Figure 6.

Maintenance of high folA gene copy number under drug pressure is driven by unstable mutants of DHFR

Though sustained drug pressure did extend the lifetime of folA gene duplication, acquisition of trimethoprim-resistant point mutations facilitated reversion to a single gene copy. We therefore asked whether folA duplication could be maintained under conditions of increasing antibiotic pressure, i.e., trimethoprim concentrations starting at 1 μg/mL and ending at 10 μg/mL over 252 generations of evolution (Figure 7 A). Under these conditions, trimethoprim IC50 values increased sharply over the course of evolution eventually reaching ∼40 to ∼160-fold higher than ancestral LTMPR1 (Figure 7 B). Colony formation at high trimethoprim concentrations also increased at later generations of evolution (Figure 7 C). Similarly evolved WTMPR4 and WTMPR5 also showed enhanced colony formation at high trimethoprim concentrations (Figure 7 C).

Figure 7.

Curiously, we observed that the 3 replicate lineages of LTMPR1 showed different trends in DHFR expression. Line A maintained a high level of DHFR protein for the duration of the experiment. On the other hand, Lines B and C showed high DHFR expression until 196 generations, after which it sharply declined (Figure 7 D). Genome sequencing revealed that all 3 lines convergently evolved the Trp₃₀Arg mutation in DHFR early during evolution and had fixed this mutation by 196 generations. Until this time point, all 3 lines retained folA duplication (Figure 7 E, S7). Subsequently, Line A evolved an additional C_-35T folA promoter mutation and had partially lost the ancestral GDA (Figure 7 E, S7, Supplementary File 5). In Lines B and C, a second missense mutation in DHFR, i.e. Tyr₁₅₁Asp, appeared between 196 and 252 generations (Figure 7 D, Supplementary File 5). Line B also showed an Ile₉₄Leu mutation in DHFR at low frequency at generation 252 (Figure 7 D, Supplementary File 5).

Unexpectedly, despite acquiring additional mutations and lower DHFR expression level, LTMPR1 Lines B and C retained two copies of folA in contrast to Line A (Figure 7 E, S7, Supplementary File 5).

Individual resistant clones from each of the lineages provided greater insight into the genetic trajectories of evolution. In LTMPR1 Line A, 2 of the 6 sequenced isolates had lost gene duplication and acquired the Trp₃₀Arg coding and C_-35T promoter mutations in folA (Figure 7 F, Supplementary File 5). Interestingly, the other 4 sequenced isolates from Line A retained the ancestral GDA and harboured the same combination of folA point mutations, but only in 1 of the 2 gene copies (Figure 7 F, Supplementary File 5). Ten of the twelve isolates from LTMPR1 Lines B and C, on the other hand, retained folA duplication. In all ten isolates, both copies of folA had acquired Trp₃₀Arg and Tyr₁₅₁Asp mutations (Figure 7 F, Supplementary File 5). Unlike LTMPR1, all sequenced isolates of WTMPR4 at 252 generations of evolution had lost their GDA mutation and gained 2 folA mutations (Figure 7 G). WTMPR5 too had gained an additional folA mutation (Figure 7 H). Notably, none of the isolates from WTMPR4 and WTMPR5 lineages had the Trp₃₀Arg and Tyr₁₅₁Asp combination of mutations that occurred repeatedly in LTMPR1 lineages.

We next sought to understand why some LTMPR1 bacteria had retained high folA copy number despite acquiring multiple point mutations in the drug target. To address this question, we first asked what the advantage of folA duplication was in the two different mutational trajectories. Isolates from Line A had similar IC50 values for trimethoprim regardless of whether they had one or two copies of folA (Figure 8 A). On the contrary, isolates from Line B that harboured only one copy of mutant folA (coding for DHFR-Trp₃₀Arg-Tyr₁₅₁Asp) displayed ∼3.5-fold lower IC50 than those that harboured two copies (Figure 8 A). Immunoblotting showed that isolates from Line A had similar expression level of DHFR (Figure 8 B). However, DHFR protein level was significantly compromised by the Trp₃₀Arg-Tyr₁₅₁Asp mutation combination and double the number of folA gene copies restored this defect (Figure 8 B). Thus, higher gene dosage of DHFR was dispensable once a promoter up-regulatory mutation was fixed. On the other hand, bacteria with alleles of DHFR harbouring multiple missense mutations continued to benefit from overexpression due to folA duplication.

Figure 8.

Mutations at position Trp₃₀ in DHFR render the protein resistant to inhibition by trimethoprim but are accompanied by lower in vivo proteolytic stability due to perturbation of the hydrophobic core of the protein [54]. This affects the expression level of DHFR, resulting in lower cellular abundance that limits the resistance level conferred by mutations at Trp₃₀ [49, 54]. Interestingly, residue Tyr₁₅₁ lies in the same hydrophobic pocket as Trp₃₀ (Figure 8 C). To test the effect of combination of mutations at these sites on trimethoprim resistance and DHFR expression, we first heterologously expressed wild type DHFR, Trp₃₀Arg or Tyr₁₅₁Asp single mutants and the Trp₃₀Arg-Tyr₁₅₁Asp double mutant in wild type E. coli (Figure 8 D). As expected, DHFR Trp₃₀Arg expression increased the IC50 of trimethoprim. Surprisingly, neither DHFR Tyr₁₅₁Asp, nor the Trp₃₀Arg-Tyr₁₅₁Asp double mutant conferred resistance to the antibiotic (Figure 8 E). Indeed, both mutants conferred lower IC50 values than wild type DHFR. This effect was due to loss of hydrophobicity at the Tyr₁₅₁ position, as DHFR-Tyr₁₅₁Leu and DHFR-Tyr₁₅₁Phe conferred similar IC50 values as wild type DHFR (Figure 8 E). Since the Trp₃₀Arg-Tyr₁₅₁Asp mutation combination had evolved in a lon-deficient background, we next expressed these mutants in E. coli Δlon (Figure 8 D). Here too, DHFR Trp₃₀Arg conferred trimethoprim resistance, while the Tyr₁₅₁Asp mutant did not (Figure 8 E). Strikingly, the combination of Trp₃₀Arg and Tyr₁₅₁Asp conferred very high-level trimethoprim resistance (Figure 8 E) demonstrating synergistic epistasis that was contingent on the absence of Lon protease. The levels of expression of DHFR mutants in wild type and lon-deficient bacteria mirrored trimethoprim resistance levels indicating that clearance by Lon altered their genotype-phenotype relationship (Figure 8 F). Based on these data, we concluded that the action of proteostatic machinery generated additional demand for expression by clearing mutant-DHFR proteins. This additional demand drove the maintenance of folA duplications in bacteria evolving at very high drug pressures.

Discussion

Ohno’s hypothesis, proposed in 1970, provided a conceptual framework to think about the evolution and eventual fate of duplicated genes [55]. Several theoretical and experimental studies have built upon the original idea to delineate the evolutionary outcomes of gene duplication [18, 56, 57]. More recently, experimental evolution of microbes has been used to test some of the suggestions made by Ohno and others. This approach has the inherent advantage of greater control over the selection pressure applied and availability of genetic tools to interrogate molecular mechanisms. In most studies, an artificial system was constructed and then evolved under defined selection pressures [33, 58–61]. Our study began with serendipitous isolation of spontaneous trimethoprim-resistant mutants of E. coli that overproduced DHFR because of a folA gene duplication. We exploited this system to understand the relationship between selection pressure and gene duplication in the context of drug-resistance evolution. The selection pressure for higher gene expression, i.e. “expression demand”, can be understood as the amount of a gene product that is required to maintain fitness in a particular environment. Antibiotics directly increase expression demand by inhibiting the activity of their targets. Our study shows that increasing gene dosage is a rapid, but transient evolutionary response to ensure that expression demand is met. On the other hand, point mutations are relatively stable routes to mitigating expression demand. Over short time scales, changes in gene dosage can counter the effects of antibiotic. However, their inherent instability results in almost inevitable replacement with phenotypically equivalent point mutations. Thus, continued antibiotic exposure delayed but could not prevent reversion to low gene copy number in most of our experiments. Further, we discover a novel association between proteostasis and gene dosage evolution. In E. coli evolving at high trimethoprim pressure, bacterial protein quality control exacerbated expression demand by clearing mutant DHFR-s. As a result, drug-resistant mutants of DHFR with low proteolytic stability drove the maintenance of folA gene duplication (Figure 9). Based on these findings we propose that a combination of factors, which include environmental and physiological, alter the relationship between gene copy number and fitness, and drive gene dosage evolution.

Figure 9.

A key finding from our study was the two-fold influence that proteostatic machinery had on gene copy number. Firstly, Lon protease-deficiency promoted GDA mutations, possibly by allowing accumulation of IS-tranposase proteins in E. coli. This observation resonates with earlier work by Nicoloff and co-workers who found that large IS-element flanked genomic duplications in tetracycline-resistant E. coli were commonly preceded by spontaneous mutations that inactivated lon [62, 63]. Further, the role of Lon protease in regulating IS-element mobility by a similar mechanism has been previously demonstrated for IS1 and IS903 [64, 65]. We note however that other mechanisms by which Lon may modulate the frequency of IS-mediated GDAs in E. coli cannot be ruled out. For instance, Lon protease homologs from several organisms directly bind DNA [66–70], though the role of this binding in DNA transactions is unclear. Similarly, Lon is a regulator of the recombination protein, RecA, by its action on LexA [71] and in Salmonella RecA stimulated the rate of gene duplications mediated by IS3 [23]. The second impact of protein quality control on GDA was the action of proteases on trimethoprim-resistant mutants of DHFR, which promoted gene duplication by generating additional demand for expression (Figure 9). Importantly, which of these two effects dominated was determined by a combination of environment, i.e. drug concentrations, and other mutations acquired at the folA locus during evolution. In particular, mutations in the folA promoter could alleviate expression demand and hence reduce selection for high gene copy number.

The interplay between GDA and point mutations during evolution that we report here agrees in many respects with a recent study by Tomanek and Guet [61]. They use an elegantly designed artificial system [33, 61] in which expression of the galK gene in E. coli was driven by a very low activity promoter. This strain was then allowed to evolve in different concentrations of galactose. Like the present study, they found that the activity of IS-elements biased evolutionary trajectories in favor of duplications. They also showed that co-occurrence of duplications and point mutations was selected only at high expression demand, in agreement with our findings. Finally, they argued that gene duplication events may constrain evolution by point mutations. Our study corroborates these results at the short term in that the lon-knockout which showed a higher propensity for GDA mutations did not evolve point mutations in folA at early time points of evolution. However, given the relative low stability of large GDAs, these were inevitably lost over longer duration of evolution and replaced with point mutations. Thus, it is likely that the different mutational profiles of GDA-biased populations may be relevant only at short durations of evolution. It is important to note that though we have focused on the effects of mutations at the folA locus, which are directly under trimethoprim pressure, many of the isolates from our evolution experiments also harboured additional mutations at other loci. For instance, mutations in rpoS were seen throughout our evolution lines and at least in drug-free media were significant contributors to bacterial fitness. It is not yet clear what the role of other mutations are likely to be under antibiotic selection, though it may be envisaged that these additional mutations play the role of media adaptations. We did not take their effects into in this study account since we expected them to be minor contributors to fitness especially during evolution under trimethoprim selection. However, their role in modulating the stability and selection of GDA mutations in drug-free/low drug environments may be more significant and remains to be investigated in the future.

In summary, the present study uncovers the complex relationship between gene copy number, point mutations and environment, and identifies proteostasis as a regulator of gene dosage evolution in bacteria. These findings open up the possibility of alternative explanations for how gene duplications may be maintained over long periods of evolution and provides insight into how gene expression demand is mitigated by multiple mechanisms in bacteria. Moreover, our findings relating proteostasis in general and Lon in particular, to adaptive mutations add to a growing body of literature that define a causal link between protein quality control and mutational paths of evolution. For DHFR itself, decreasing the stringency of protein quality control changes the fitness landscape of single mutations and their combinations, altering the mutational repertoires accessible during adaptive evolution [49, 50, 72, 73]. Similar effects of protein homeostatic machinery on mutant fitness are reported from a diverse range of proteins that influence phenomena like viral evolution [74] and the fate of horizontally acquired genes [75–77]. Indeed, a central role for protein quality control machinery in molecular and organismal evolution is reflected in conserved signatures of the proteostatic network across the different kingdoms of life [78]. These links may be important to consider for a new wave of potential antibacterials that target proteostasis [79–81]. It appears critical to ensure that the impact of changes in evolvability are considered when designing such strategies for therapeutic intervention.

Materials and methods

Strains, plasmids and culture conditions

Escherichia coli K-12 MG1655 was used as wild type for all experiments. E. coli DH10B was used for plasmid construction and manipulation. E. coli strains were cultured in Luria-Bertani (LB) broth with shaking at 180-200 rpm or on LB agar plates at 37 °C. Ampicillin was used at a concentration of 100 μg/mL for selection of plasmid-transformed strains. Trimethoprim was added to growth media at required concentrations as needed. The strains and plasmids used in this study are listed in Table 2.

Table 2.

Laboratory evolution of antibiotic resistance

Short-term evolution of antibiotic resistance

For short term laboratory evolution of antibiotic resistance, wild type and E. coli Δlon populations were cultured in 96-well plates in LB medium incorporating sub-inhibitory concentrations of different antibiotics. Following concentrations of antibiotics were used: trimethoprim (300 ng/mL; MIC/3); spectinomycin (5 μg/mL; MIC/10); nalidixic acid (4.25 μg/mL; MIC/3); erythromycin (10 μg/mL; MIC/10); rifampicin (4.25 μg/mL; MIC/3). Control populations were also set up in antibiotic-free LB. Five replicates were included for each strain in each antibiotic. Total volume of culture in each well was 150 μL. 10 % of the population was passaged to fresh growth media every 24 hours, for 7 days, resulting in ∼25 generations of evolution. At the end of 7 days, the evolved populations were frozen as glycerol stocks.

To check for the evolution of resistance and isolation of resistant clones for further analyses, frozen populations were revived in drug-free LB broth and grown overnight (14-16 hours). Saturated cultures were serially diluted and spotted on LB agar plates supplemented with the appropriate antibiotic at MIC. Colonies were randomly picked, cultured overnight in drug-free LB and frozen in 50% glycerol at-80 °C until further use.

Long-term evolution at different trimethoprim pressures

For long-term evolution, trimethoprim-resistant isolates LTMPR1, WTMPR4 and WTMPR5 were first revived from frozen stocks in drug-free LB overnight. Triplicate evolving lineages were then established by sub-culturing 1% of the saturated culture of each ancestor into 2 mL of appropriate medium. For relaxed selection drug-free LB was used. For constant antibiotic pressure, trimethoprim was added at a concentration of 300 ng/mL. For increasing antibiotic pressure, trimethoprim was increased from 1 μg/mL to 10 μg/mL over course of the experiment as shown in Figure 7 B. Cultures were grown for 10-14 hours (∼7 generations) and then passaged (1%) into fresh media. This was continued for 36 passages, i.e. ∼252 generations. Aliquots of cultures were frozen periodically for further analyses.

For isolating individual clones, populations were revived in drug-free LB from frozen stocks at 252 generations of evolution in drug-free LB, serially diluted and then spotted onto LB agar plates supplemented with trimethoprim (1, 5 or 10 µg/mL). Plates were incubated for 18-20 hours at 37 °C and resistant colonies were picked randomly for further analyses.

Measuring trimethoprim resistance

Trimethoprim IC50

Trimethoprim Inhibitory Concentration - 50 (IC50), i.e. the antibiotic concentration needed for 50 % growth inhibition, was calculated using a broth dilution assay. Frozen glycerol stocks of appropriate bacterial strains were revived in 1-3 mL LB for 6-8 hours. In a 96-well plate, 150 μL of LB was added to individual wells. Trimethoprim stocks were added, and the antibiotic was serially diluted. No antibiotic was added to a few wells, which served as controls. Next, 1.5 μL of the appropriate bacterial culture was added to each well. The peripheral wells of the plate were filled with water to prevent dehydration. The 96-well plate was then incubated at 37LJ for 15-18 hrs. Optical density at 600 nm was measured for each well using a plate reader (Ensight, Perkin Elmer). Obtained data were first normalised to growth in drug-free wells (set to 1) and then fitted to a sigmoidal log(inhibitor)-response curve using Graphpad Prism software. The equation for this model is as follows:

Y= Min + (Max-Min)/(1+10 ^{(X-Log(IC50))})

where, Min: Minimum value of normalised OD and Max: maximum value of normalised OD Values of IC50 were estimated from the above curve fit and averaged across 3 independent measurements for every strain.

Colony forming efficiency

To calculate the colony forming efficiency of evolving populations at different trimethoprim concentrations, frozen glycerol stocks were first revived in 1 mL LB for 14-15 hours. Saturated cultures were serially diluted (10-fold). Diluted cultures were then spotted (10 mL) on LB agar plates containing 0, 1, 5 and 10 μg/mL of trimethoprim. Spots were allowed to dry completely and the plates were incubated at 37 °C for ∼18 hours. Colonies were counted from the least dense spot and multiplied by the appropriate dilution factor to estimate colony forming units per mL (CFUs/mL) of culture.

Growth curves

Appropriate frozen stocks were revived overnight in drug-free LB (1-3 mL) and 1μl of the saturated culture was added to 200 μL of media supplemented with trimethoprim (300 ng/ml) in the wells of a sterile uncoated 96-well plate. Optical density at 600 nm was measured during growth at 37 °C at 15-minute intervals over a 24-hour period on a Perkin-Elmer Ensight Multimode Plate Reader. Doubling time (DT) (in mins) was estimated from log phase, using the following formula:

Where OD1 and OD2 represent optical densities at the onset and end of log phase respectively, and T1 and T2 represent time elapsed (in mins), corresponding to OD1 and OD2.

Genome sequencing, identification of mutations and detection of GDA mutations

Genomic DNA was extracted from late log-phase cultures of appropriate strain/resistant isolates/evolving populations using phenol:chloroform:isoamyl alcohol (25:24:1) extraction as previously described [82]. Extracted genomic DNA was then cleaned-up using HiGenoMB DNA isolation and purification kit (Himedia, India) as per manufacturer’s instructions and quantitated spectrophotometrically. The integrity was checked by running purified DNA on an agarose gel before further analyses.

For whole-genome resequencing, paired-end whole-genome next-generation sequencing (NGS) was performed on a MiSeq system (Illumina, USA) with read lengths of 150-250 bps. Library preparation and sequencing services were provided by Eurofins (Bangalore, India) and Genotypic (Bangalore, India). Processed reads were aligned to the reference genome E. coli K-12 MG1655 (U00096) using bowtie2. All sequenced genomes had an average coverage depth of at least 100x. Variant calling and prediction of new junctions in the genome to identify large structural mutations was done using Breseq [83] using default settings in population or clonal modes as needed. All variants and new junctions that were already present in the ancestor were excluded. Variants unique to the test sample were classified as “mutations”.

For detection of GDA mutations, ‘coverage’ vs ‘coordinate in reference genome’ plot generated by the bam.cov command in Breseq was used [83]. Contiguous stretches greater than 10 kb with coverage depth of 1.5x or more, compared to the rest of the genome, were classified as GDA mutations. The presence of IS-elements or repeat sequences was manually checked using the coordinates of the GDA.

Calculating relative fitness of E. coli strains

Relative fitness of gene knockouts or evolved clones was estimated using competition assays. Saturated cultures of test strains and reference strains (E. coli ΔlacZ or E. coli ΔlacZΔlon) were mixed 1:1 to a total of 1% in 3 mL LB. Strains were allowed to compete for ∼7 generations (24 hours incubation). Mixed cultures were diluted and plated on LB agar supplemented with IPTG (1 mM) and X-gal (50 μg/mL) at the start and end of the competition experiment. Plates were incubated for 18-24 hrs, and blue and white colonies were counted to determine the CFUs/mL at initial and final time points.

Relative fitness (w) was calculated using the following formula:

where T_f and T_i are the final and initial CFUs/mL of the test strain and R_f and R_i are the final and initial CFUs/mL of the reference strain, respectively.

Estimation of DHFR expression by immunoblotting

Immunoblotting was used to determine the expression level of DHFR in various strains/mutants of E. coli as described in Patel and Matange (2021) [46]. Briefly, equal lysate protein (5 μg) was subjected to SDS-PAGE on a 15% gel and electroblotted onto PVDF membranes. The membrane was blocked with 5% BSA solution and cut into 2 halves at the 25 kDa molecular weight band. The lower part of the membrane was probed with anti-DHFR polyclonal IgG (100 ng/mL) and the upper part was probed with anti-FtsZ serum (1:50,000). HRP-linked anti-rabbit IgG (1:10,000) was used as the secondary antibody. Hybridised antibodies were detected using chemiluminescence. Band intensities were quantitated in ImageJ and normalised to appropriate control (set to 1). Fold change in DHFR was calculated as follows:

Gene expression analysis using RNA-sequencing

Transcriptome analyses were performed using RNA sequencing. For RNA-seq, saturated overnight culture of appropriate strains/isolates was inoculated (1 %) into 3 mL LB in triplicate and incubated at 37 °C for 3 hours with shaking at 180 rpm. The cultures were centrifuged and bacterial pellets were resuspended in 1 mL of RNAlater (Invitrogen) for storage. RNA extraction, sequencing and preliminary data analyses were performed by Eurofins (Bangalore, India). Total RNA was extracted using RNAeasy spin columns and quantitated using Qubit and Bioanalyzer. Bacterial rRNA was depleted using Ribozero kit. Bacterial mRNA was then reverse transcribed, library prepared and quality control performed using Tapestation platform. Paired end sequencing was performed on Illumina platform with read lengths of 150 bp read length. Processed reads were aligned to the reference genome (U00096.3) using HISTAT2(version 2.1.0). Abundance estimation was done using featureCounts(version1.34.0). Differential gene expression (DGE) analysis was done by comparing the expression of individual genes in test and reference samples. The output of the DGE analyses were log₂(fold change) values for each gene and P-values for statistical significance.

Cloning and site-directed mutagenesis

Cloning of InsB and InsL transposases

InsB was amplified from E. coli K-12 MG1655 genomic DNA by PCR using primers insB1_fwd_NcoI (5’-GGCGCCATGGATATGCCGGGCAACAGCCCG-3’) insB1_rev_EcoRI (5’- CCTTTGAATTCTTATTGATAGTGTTTTATG-3’). For amplification of insL, primers insL_EcoRI_fwd (5’- ATGCCGAATTCATGAATTACTCTCACGATAAC- 3’) and insL_HindIII_rev (5’- ATGCCAAGCTTTTAGTTCTTCTTTTCGGATCC-3’) were used. Amplified insB gene was digested with NcoI and EcoRI restriction endonuclease and ligated with similarly digested pPROEx-HtC plasmid. Amplified insL was digested with EcoRI and HindIII restriction endonucleases and ligated with similarly digested pPROEX- HTa plasmid. Positive clones were screening using restriction digestion. Clones were transformed into E. coli wild type or Δlon for further analyses.

Site-directed mutagenesis for DHFR

Site-directed mutagenesis was performed to introduce mutations into the coding sequence of folA in pPRO-folA as described in [54]. Briefly, PCR was performed using pPRO-folA or pPRO-folA Trp₃₀Arg as template and the appropriate forward and reverse mutagenic primers given below:

folA_Tyr₁₅₁Asp_fwd_BglII:

5’-CTCTCACAGCGATTGCTTTGAGATCTTGGAGCGGCGGGGCC-3’

folA_Tyr₁₅₁Asp_rev_BglII:

5’-CGCCGCTCCAAGATCTCAAAGCAATCGCTGTGAGAGTTCTGCG-3’

folA_Tyr₁₅₁Leu_fwd_BglII:

5’-TCTCACAGCTTATGCTTTGAGATCTTGGAGCGGCGGGGCC-3’

folA_Tyr₁₅₁Leu_rev_BglII:

5’- CGCCGCTCCAAGATCTCAAAGCATAAGCTGTGAGAGTTCTGCG-3’

folA_Tyr₁₅₁Phe_fwd_BglII:

5’- CTCACAGCTTTTGCTTTGAGATCTTGGAGCGGCGGGGCC-3’

folA_Tyr₁₅₁Phe_rev_BglII:

5’- GCCGCTCCAAGATCTCAAAGCAAAAGCTGTGAGAGTTCTGCG-3’

PCR products were digested with DpnI overnight, transformed into chemically competent E. coli DH10B and selected on LB agar plates containing ampicillin. Plasmids habouring the desired mutations were screened by restriction digestion and confirmed by sequencing. Sequenced clones were transformed into E. coli wild type or Δlon for further analyses.

In-cell stability of InsB and InsL

A chloramphenicol-chase assay was used to assess the in-cell stability of heterologously expressed InsB and InsL transposases as described in Matange (2020) [49] with a few modifications. Overnight cultures of E. coli wild type or Δlon harbouring plasmids pPRO- insB or pPRO-insL were diluted 1% in 5 mL LB supplemented with ampicillin and IPTG (0.1mM). After 4 hours of growth at 37 °C, 1 mL of culture was pelleted. The cell pellet was lysed in 100 μL of 1x Laemmli sample loading buffer by heating at 95°C for 10-20 mins. Chloramphenicol was added to the remaining culture at a final concentration of 25 μg/mL to stop protein synthesis and the culture was placed back in the incubator. After 15, 45, and 90 minutes 1 mL aliquots of the culture were pelleted and lysed. Lysates were subjected to SDS- PAGE followed by immunoblotting to detect the levels of InsB and InsL. Anti-Hexahis monoclonal antibody (Invitrogen) was used as the primary antibody for InsL. Anti-Hexahis polyclonal antibody (Invitrogen) was used as the primary antibody for InsB. Band intensities were calculated in ImageJ and normalised to control (set to 1) to calculate fraction of protein remaining after chloramphenicol treatment.

Supplementary Figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Acknowledgements

We acknowledge Dr. Manjula Reddy for providing us with the anti-FtsZ polyclonal antibody. We acknowledge Bruhad Dave for his valuable inputs during discussions at early stages of this project and Prof. Satyajit Rath, Prof. Sandhya Visweswariah and Dr. Girish Deshpande for their comments on the manuscript. Funding for this work was provided by DBT/Wellcome Trust India Alliance. NM is a recipient of the India Alliance Intermediate Fellowship. CJ is a recipient of a Senior Research Fellowship from the Council of Scientific and Industrial Research (CSIR), Govt. of India.

Additional files

Supplementary File 1. Gene expression changes in LTMPR1 and WTMPR4 isolates.

Supplementary File 2. Genomic changes associated with evolution of LTMPR1 in drug-free media.

Supplementary File 3. Gene expression changes in isolates derived from LTMPR1 evolution in drug-free media containing mutation in rpoS or GDA-2.

Supplementary File 4. Genomic changes associated with evolution of LTMPR1 in trimethoprim (300 ng/mL).

Supplementary File 5. Genomic changes associated with evolution of LTMPR1 in increasing trimethoprim pressure.

Supplementary File 6. Genome sequencing of isolates derived from WTMPR4 evolution at different trimethoprim pressures.

Supplementary File 7. Genome sequencing of isolates derived from WTMPR5 evolution at different trimethoprim pressures.