For the adaptation evolution experiment and the generation of samples for collection in this study, the bacterial strain used was MG1655, which is an E. coli K-12 MG1655 wild type (Blattner et al., 1997). Further characterization experiments were performed on successfully recombined strains to reintroduce given mutations (Table 1).

Microtiter plates saved from Munck et al. (2014) were defrosted for 30 min at 4°C. To prevent aerosol contamination, the plates were centrifuged at 300 ×g for 5 s before seal removal. To revive our samples, 150 μL were transferred to new tubes with 500 μL LB media and incubated at 37°C for 3 h. Genomic DNA was extracted using Genomic Mini Kits in accordance with the manufacturer’s protocol (A&A Biotechnologies, Gdynia, Poland).

The loci that were selected for sequencing were chosen based on knowledge of their involvement in antibiotic resistance and their mutation profiles in the end-point sequenced original evolved lineages, and the population level loci sequencing was performed on days 2, 4, 6, 8, 10 and 12 (Figure 1A). For this study, we monitored the following mutations of AMK: cpxA, fusA, and sbmA; for CHL, we monitored the following mutations: acrR, marR, rob, and soxR; for CIP, we monitored the following mutations: acrR, gyrA, and soxR; for AMK + CHL, we monitored the following mutations: acrR, cpxA, fusA, marR, rob, and soxR; and for CHL + CIP, we monitored the following mutations: acrR, gyrA, marR, rob, and soxR. We believe that these loci cover important genomic locations that could begin to unravel important interactions in intrinsic antibiotic resistance. The primers were designed with overhangs to allow for barcoding with Nextera XT indices in accordance with the Illumina 16S Metagenomic Sequencing Protocol (Illumina Protocol Online) to amplify 12 different genomic regions (Supplementary Table S1) in areas with known resistance-causing SNPs. All primers were tested on E. coli-K12 wild-type genomic DNA, and the PCR products were verified by 1% agarose gel electrophoresis with GelRed (Biotium, Freemont, CA, United States).

All samples were amplified using the primers listed in Supplementary Table S1 according to the descriptions in the above section for each drug condition. Each PCR reaction was performed in 0.2-mL sterile PCR tubes that contained 10-μL Phusion Flash High-Fidelity PCR Master Mix (Thermo Fisher Scientific, United States), 0.5 μM forward primer, 0.5 μM reverse primer, 1 μL DNA template and H₂0 to a total reaction volume of 20 μL. The PCR amplification consisted of an initial denaturation at 98°C for 30 s, 30 cycles consisting of 98°C for 10 s denaturation, 65°C for 10 s annealing, 72°C for 15 s elongation, 72°C for 60 s for final extension, and subsequent holding at 4°C. To validate correct amplification, gel electrophoresis was performed with 1% agarose gel stained with ethidium bromide. After verification, all the different locus amplicons that belonged to the same drug condition, lineage and time point were pooled together. The PCR products were cleaned with an Agencourt AMPure XP PCR purification system according to the manufacturer’s online protocol (Beckman Coulter 2013). The barcoding of the amplicons was achieved via a second PCR amplification that was performed in 0.2-mL sterile PCR tubes that contained 25 μL Phusion Flash High-Fidelity PCR Master Mix (Thermo Fisher Scientific), 5 μL Nextera XT index 1 (N7XX; Illumina), 5 μM Nextera XT index 2 (S5XX; Illumina), 5 μL of the pooled amplicons and H₂0 to a total reaction volume of 50 μL. The PCR amplification used an initial denaturation of 98°C for 30 s, 8 cycles consisting of 98°C for 10 s denaturation, 65°C for 10 s annealing, 72°C for 15 s elongation, 72°C for 60 s for final extension and subsequent holding at 4°C. The barcoded amplicons were cleaned according to the instructions of Agencourt AMPure XP PCR purification systems (Beckman Coulter 2013). The barcoded amplicons were measured for their DNA concentrations with Qubit dsDNA HS Assay kits on a Qubit Fluorometer (Thermo Fisher Scientific Inc. 2015) and average fragment size with Agilent DNA 1000 kits on an Agilent 2100 bioanalyzer instrument (Agilent Technologies 2013). Next, the barcoded amplicons were pooled based on fragment size and DNA concentration value to form a sequencing library. The sequencing library was sequenced on the MiSeq sequencing platform (Illumina, United States).

All raw sequencing reads from our amplicon sequencing have been deposited in the Sequence Read Archive (SRA) under the BioProject accession number PRJNA328094. The sequencing reads were mapped to the reference genome (GenBank accession NC_000913) using the CLC genomics workbench (Qiagen). These files were exported and used for coverage analysis. The CLC genomic workbench (Qiagen) and the statistical program R (R Core Team, 2014) were used to verify the sequencing coverage mapping read files. Coverage plots were generated, and any amplicons that were below our set threshold of 1000 read coverage were re-done until our threshold was achieved to ensure good quality data (Supplementary Figures S1–S5). To detect variants, all mappings were analyzed with the basic variant detection calling function of the CLC Genomic Workbench (Qiagen). The files were further analyzed using the program R (R Core Team, 2014). Our detected variants are listed in Supplementary Table S2.

For the single mutant variants, the ancestral wild-type E. coli K-12 MG1655 was transformed with the transient mutator plasmid PMA7sacB (GenBank accession KT285941) (Lennen et al., 2016). For each single variant mutation, two recombineering cycling rounds were performed with E. coli K-12 MG1655 harboring the PMA7sacB plasmid according to the MAGE cycling with the single-stranded oligonucleotides (ss-oligos) with verification primers (Supplementary Table S3) and plasmid curing described by Lennen et al. (2016). Direct mutagenesis verification was confirmed by Sanger sequencing performed with Eurofins Mix2Seq kits according to the manufacturer’s instructions (Eurofins Genomics, Esbersberg, Germany). After verification, the double-mutant variants were created by taking one of the single-mutant variants from the observed double-mutant pair and transforming it with the PMA7sacB plasmid. Recombineering with the other mutant variant ss-oligo was performed in two rounds as previously described. Plasmid curing was also performed.

The relative growth rate was used as an indicator of fitness and was assessed for all the recombineered strains and wild type. To measure relative growth fitness, overnight cultures in LB medium were diluted to 1 × 10⁶ CFU/mL and then added to 96-well microtiter plates with a medium negative control. All procedures were performed in eight biological replicates with two technical replicates to allow for the calculation of the relative growth rate. The growth of the samples was done at 37°C with shaking for 10 s and kinetic measurements on ELx808 absorbance microplate reader with Gen5 software (BioTek Instruments Inc., Winooski, VT, United States). Optical density measurements were taken at 630 nm (OD₆₃₀) every 10 min. Exponential growth was defined by OD₆₃₀ values between 0.02 and 0.1. To calculate relative growth, R statistical software was used to find the doubling time (Td) of each of the strains and the standard deviation (Td SD) and then to calculate relative growth rate, and the recombineered strain Td value was divided by the wild-type Td value. Selective advantage was measured by performing IC90 determination assays with the antibiotic drug condition in which the mutations were observed (R Core Team, 2014). IC90 determination was performed in 96-well microtiter plates with a two-fold drug gradient that consisted of 10 drug concentrations that were prepared in LB broth; each well had a volume of 100 μL and was inoculated with 1 μL bacterial culture to provide each well with an initial inoculum of 1 × 10 CFU/mL. Each plate had eight biological replicates. The plates were incubated for 18 h at 37°C, and the OD₆₀₀ values were read on a Synergy H1 plate reader (BioTek Instruments Inc., Winooski, VT, United States). The data were analyzed with Prism and IC90 values were calculated with point-to-point curve analysis Prism 6 (GraphPad Software Inc.).

To calculate the doubling time of each of the bacterial strains, overnight cultures in LB medium were diluted to 1 × 10⁶ CFU/mL and then added to 96-well microtiter plates to which a two-fold serial dilution of antibiotics was added with a medium negative control; all were performed with four biological replicates. Growth of the samples was performed at 37°C with shaking for 10 s, and kinetic measurements performed on ELx808 absorbance microplate reader with Gen5 software (BioTek Instruments Inc., Winooski, VT, United States). Optical density measurements were taken at 630 nm (OD₆₃₀) every 10 min. Exponential growth was defined as OD₆₃₀ values of between 0.02 and 0.1. To calculate relative growth, R statistical software was used to find the doubling time (Td) of each of the strains, as well as the standard deviation (Td SD) (R Core Team, 2014). Plots were generated with prism 6 (GraphPad Inc.).

We specifically targeted putative resistance loci that were identified by Munck et al. (2014) through the whole-genome sequencing of evolved isolates. As a template, we used DNA isolated from the evolving populations that were sampled every 48 h during the experiment. The following genomic loci used for each drug condition based on previously end-point sequenced mutant alleles were: cpxA, fusA, and sbmA (AMK); marR and rob (CHL); acrR, gyrA, and soxR (CIP); cpxA and marR (AMK + CHL); and acrR, marR, gyrA, and soxR (CHL + CIP). These genes are implicated in either multi-drug resistance or specific-drug resistance. The genes implicated in multi-drug resistance include global regulators that induce the up-regulation of acrAB-TolC multi-drug efflux systems, i.e., acrR, marR, rob, and soxR (Chubiz et al., 2012). Mutations in these global regulators often have generalized cross-resistance effects across multiple antibiotic classes. By contrast, the genes implicated in specific-drug resistance are gyrA mutations (Fu et al., 2013) that confer CIP resistance and mutations of cpxA (Girgis et al., 2009), fusA (Johanson and Hughes, 1994), and sbmA (Chubiz et al., 2012; Puckett et al., 2012) that confer AMK resistance.

Our sample collection consisted of 90 samples obtained from 15 bacterial populations (triplicate lineages of five drug conditions). Three of these drug conditions involved single drugs, i.e., AMK, CHL, and CIP, and the remaining two drug conditions involved the drug combinations of AMK–CIP and CHL–CIP. The triplicate parallel lineages were assigned letters from A–C and were sampled at six time points (days 2, 4, 6, 8, 10, and 12; Figure 1A). We extracted genomic DNA from these samples, amplified the loci of interest, verified the products by gel electrophoresis, pooled the amplicons comprising targeted alleles based on drug condition and time point, barcoded the amplicons with indices, and pooled these together to form a sequencing library. The library was sequenced using the MiSeq sequencing platform to establish the frequency of mutants in each of the selected loci (see Materials and Methods and Figure 1B). A minimum coverage of 1000 sequencing reads per amplicon was required before further data processing (Supplementary Figures S1–S5).

From this data, we compared the allele dynamics of the single-drug evolved populations with those of the combination drug-evolved populations. It is assumed that antibiotic combinations can suppress antibiotic resistance evolution compared with drug monotherapy. However, we have previously shown that this is not always the case; rather, only combinations with collateral sensitivity interactions seem to suppress drug evolution (Munck et al., 2014). We performed a Holm–Sidak multiple comparison two-way ANOVA statistical test to compare the cumulative frequency of mutations for the sequenced loci between the different drug conditions (Figure 2). We found spurious differences between the single- and combination-drug conditions in cumulative frequency of mutations for specific loci. However, two loci exhibited consecutive significant differences over two or more time points. These were the fusA gene and the sbmA gene in which mutations were significantly less frequent in the AMK–CHL drug-evolved populations compared with the AMK drug-evolved populations (P-value ≤ 0.001, Holm–Sidak multiple comparison two-way ANOVA test; Figure 2B). Similar results were observed in the work of Munck et al. (2014) and likely results from collateral sensitivity to CHL caused by these mutations (Macvanin et al., 2005; de Cristobal et al., 2008).

To focus on important alleles within our populations, we selected all alleles that appeared in at least one time point with a population frequency ≥ 30% (Supplementary Figure S6). We observed five mutations that occurred in parallel evolution events in different adapted populations (AMK: cpxA Trp184Arg, fusA Pro610Gln; AMK–CHL: acrR Lys55Glu; CIP and CHL–CIP: gyrA Ser83del; CHL–CIP: soxR Arg20Cys). Of these mutations, two have not previously been reported in the literature, namely, acrR Lys55Glu in AMK–CHL populations and cpxA Trp184Arg in AMK populations. We speculate that the acrR Lys55Glu mutation confers resistance to the CHL component in the AMK–CHL populations via the up-regulation of the acrAB-TolC multi-drug efflux pump, whereas cpxA Trp184Arg confers AMK resistance by modifying the membrane stress response.

The most prevalent mutation in the evolved bacterial populations was gyrA Ser83del, which was observed in three CHL-CIP and two CIP populations at several time points. The gyrA gene often mutates at amino acid position 83; however, it has only been reported once that the mutation was due to a deletion of the codon (Jaktaji and Mohiti, 2010), as seen within our data. We observed a high prevalence of this mutation, especially in the CHL-CIP lineage A; on day 12, this mutation was observed at a frequency of 100% in the population. However, this mutant allele was not detected by Munck et al. (2014) possibly due to the use of whole-genome sequencing of a representative isolate from the population that was performed on day 14. To support our findings, we have provided a picture of our population result (Supplementary Figure S7) and an additional Sanger sequencing result of the population to re-confirm our findings (Supplementary Figure S8).

Within our data, we observed three distinct allele dynamics (Figure 3A), which we termed as follows: interfering allele dynamics in which the mutant alleles have opposing dynamic trajectories (Figure 3A); super-imposing allele dynamics in which the mutant alleles have the same frequency trajectories (Figure 3B); and fixation allele dynamics in which one allele fixes in the population and provides the genetic background in which new alleles appear (Figure 3C).

We observed interfering allele dynamics in the following three populations: CIP lineage B, AMK lineage B, and CHL Lineage C. We assume that clonal interference is the most likely cause of our observations, but in some cases, it could be due to mutational incompatibility. For instance, we were able to demonstrate direct antagonistic interactions between the two mutations in the CIP lineage B because both mutations were captured within a single read. We observed that the individual reads from this lineage either had an amino acid change in gyrA at position Ser83del or Asp87Gly for the time points from 4 to 12. However, further experimentation is required to deduce which effect is responsible for our other observations of interfering allele dynamics.

We observed super-imposing allelic dynamics in the following three populations: AMK lineage C, CHL–CIP lineage A, and CIP lineage A. For all three populations, the two mutant alleles exhibited similar population frequencies over multiple time points (Figure 3C). Such fixation patterns likely result from the quick succession of two separate mutations within a cell. Because our population sequencing was performed only every 48 h, we likely were unable to resolve the successive occurrence of these mutations. These observations could be the result of two possible effects. The first is the hitchhiking effect in which a neutral mutant allele co-occurs with a beneficial mutant allele and, due to the selective advantage provided by the beneficial allele, the neutral mutation is therefore irreversible (Taddei et al., 1997; Tenaillon et al., 1999). The second is that both mutations contribute to the selected phenotype (Blount et al., 2008; Fogle et al., 2008; Wong and Seguin, 2015). We are not able to make a distinction between these two scenarios based on our allele frequency data alone.

For the fixation allele dynamics, we observed one allele that fixes at 100% in the population followed by the establishment of a subsequent allele. This pattern was observed in the following three populations: AMK lineage B, CHL lineage B, and CIP Lineage A.

In summary, we observed interfering allele dynamics in five populations, superimposing allele dynamics in seven populations interactions and fixation allele dynamics in eight populations (Table 2). This highlights that diverse population dynamic patterns can be observed during laboratory antibiotic resistance evolution.

To follow up on our observations, we used recombineering to engineer strains with our observed mutant alleles. This process was performed with a subset of detected mutations that exhibited superimposing establishment or interfering allele dynamics. By measuring the relative growth rates and antibiotic tolerance of the isogenic mutants compared with those of the ancestral wild type, we were able to deduce the basis of some allele dynamics (Figure 4).

In the CHL lineage C population, we observed interfering allele dynamics. In particular, we observed that the rob G467A allele became dominant in the population and that the marR C328T was lost from the population from day 6 onward (Figure 3A). Therefore, we constructed three bacterial strains: MG1655: rob G467A, MG1655: marR C328T, and MG1655: marR C328T+rob G467A, and tested the relative growth rate and CHL IC90 for each strain relative to the wild type (Figure 4A). We find that the rob G467A strain was selected due to its increased antibiotic resistance compared with the marR C328T strain (Figure 4A). The rob G467A strain exhibited a 6-fold increase in the IC90 relative to the ancestral wild type, whereas the marR C328T strain exhibited a 5.5-fold increase. Even though this difference in IC90 values is small it suggests that the observed fixation patterns in the CHL lineage C population resulted from clonal interference and not mutational incompatibility. To validate this assumption, we constructed the double mutant marR C328T+rob G467A; this strain was viable and exhibited CHL IC90 11-fold higher than the ancestral wild type. The double mutant did exhibit a reduced growth rate in absence of antibiotics compared to the single mutants and the ancestral wild type. However, when quantifying the growth rate of the strains in the presence of different CHL concentrations we observed that the double mutant was able to sustain growth up to concentrations of 8 mg/L CHL (doubling time of 54 min), whereas both single-mutant allele strains could only sustain growth in CHL concentrations up to 2 mg/L (the doubling time for marR C328T and rob G467A was 40 and 60 min, respectively). Finally, the ancestral wild type could sustain growth in CHL concentrations up to 2 mg/L with a doubling time of 57 min (Figure 4B). Accordingly, the double mutant could have been selected for in the experiment at CHL concentrations above 2 mg/L.

For our super-imposing allele dynamics, we tested our observations in the CHL–CIP lineage A population in which gyrA 247_249delTCG and marR T251A appeared at the same frequency from day 8 (Figure 3B). These observed mutations likely complement each other since gyrA 247_249delTCG would confers target mediated resistance against CIP (Fu et al., 2013), and marR T251A induces up-regulation of the acrAB-TolC multi-drug efflux pump that causes collateral resistance to both CHL and CIP (Chubiz et al., 2012). Therefore, we recombineered three bacterial strains, i.e., MG1655: gyrA 247_249delTCG, MG1655: marR T251A, and MG1655: gyrA 247_249delTCG + marR T251A, and quantified the doubling time and CHL–CIP tolerance relative to the ancestral wild type (Figure 4A). We found that the double mutant exhibited strong selective advantage compared to the single mutants and the ancestral wild type in the presence of CHL–CIP. The IC90 of the double mutant was 11-fold higher than that of the ancestral wild type, and the corresponding values were 4.3-fold and 6.6-fold for the single mutants gyrA 247_249delTCG and marR T251A, respectively. The doubling time of the double mutant was 59% higher than the ancestral wild type when grown without antibiotic. Yet, in the presence of CHL–CIP the double mutant had a substantially lower doubling time than the single mutants and the ancestral wild type (Figure 4). The double mutant was able to sustain growth at concentrations of CHL–CIP up to 8.0116 mg/L (doubling time of 98 min). Both of the single-mutants could sustain growth at concentrations of CHL–CIP up to 2.004 mg/L (doubling time for marR T251A and gyrA 247_249del TCG was 58 and 124 min, respectively). The ancestral wild type sustained growth at CHL–CIP concentrations up to 2.004 mg/L (doubling time of 89 min). Based on these data it is obvious that if a double mutant was present in a population containing both single-mutants and wild type cells it would rapidly outcompete these other strains as antibiotic exposure increased.