Experimental evolution in morbidostat reveals converging genomic trajectories on the path to triclosan resistance.

Understanding the dynamics and mechanisms of acquired drug resistance across major classes of antibiotics and bacterial pathogens is of critical importance for the optimization of current anti-infective therapies and the development of novel ones. To systematically address this challenge, we developed a workflow combining experimental evolution in a morbidostat continuous culturing device with deep genomic sequencing of population samples collected in time series. This approach was applied to the experimental evolution of six populations of Escherichia coli BW25113 towards acquiring resistance to triclosan (TCS), an antibacterial agent in various consumer products. This study revealed the rapid emergence and expansion (up to 100% in each culture within 4 days) of missense mutations in the fabI gene, encoding enoyl-acyl carrier protein reductase, the known TCS molecular target. A follow-up analysis of isolated clones showed that distinct amino acid substitutions increased the drug IC90 in a 3-16-fold range, reflecting their proximity to the TCS-binding site. In contrast to other antibiotics, efflux-upregulating mutations occurred only rarely and with low abundance. Mutations in several other genes were detected at an earlier stage of evolution. Most notably, three distinct amino acid substitutions were mapped in the C-terminal periplasmic domain of CadC protein, an acid stress-responsive transcriptional regulator. While these mutations do not confer robust TCS resistance, they appear to play a certain, yet unknown, role in adaptation to relatively low drug pressure. Overall, the observed evolutionary trajectories suggest that the FabI enzyme is the sole target of TCS (at least up to the ~50 µm level), and amino acid substitutions in the TCS-binding site represent the main mechanism of robust TCS resistance in E. coli. This model study illustrates the potential utility of the established morbidostat-based approach for uncovering resistance mechanisms and target identification for novel drug candidates with yet unknown mechanisms of action.

Data Summary

Illumina sequencing data in FastQ format have been deposited at the National Center for Biotechnology Information (NCBI) Sequence Read Archive [1] and can be downloaded from the NCBI BioProject database using accession number PRJNA472810 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA472810). The supplementary data (text, figures and tables) are available on Figshare (https://doi.org/10.6084/m9.figshare.13076264). A detailed description of the morbidostat implementation, the scripts for NGS data processing and the tool for calling IS element reallocation (iJump) can be found in the following GitHub repositories: https://github.com/sleyn/morbidostat_construction, https://github.com/sleyn/Triclosan_EE_paper and https://github.com/sleyn/ijump, respectively.

The current article describes the application of a morbidostat-based experimental evolution workflow for studying the dynamics and mechanisms of acquired antibiotic resistance, which is critical for the development and optimization of antimicrobial therapies. The established workflow enables accelerated studies of antimicrobial resistance via continuous culturing of target bacteria under automatically adjustable drug pressure combined with monitoring of the dynamics of mutation accumulation via deep sequencing of evolving bacterial populations in time series. We report the validation of this workflow in a model system of evolving resistance to a common biocide triclosan (TCS), which revealed distinct evolutionary trajectories in six parallel reactors converging to a series of amino acid substitutions in the active site of enoyl-acyl carrier protein reductase FabI, the known TCS target. The practical significance of this work lies in the detailed description of the hardware and software components of a ready-to-implement workflow enabling antibiotic resistance studies for novel drug candidates.

Introduction

Antibiotic resistance is a major threat to public health. According to US and European Centers for Disease Control and Prevention reports, infections by drug-resistant bacteria cause more than 30 000 deaths per year [2, 3]. At the same time, the current pace of development and regulatory approval of new antibiotics is too slow to curb the spread of acquired multidrug resistance. Thus, of 51 antibiotics and 11 biologics counted as being in development by the World Health Organization in 2017, only 10 are expected to be approved within the following 5 years [4]. Knowledge of the dynamics and mechanisms of acquired drug resistance across a range of target bacterial pathogens is one of the bottlenecks in the development and approval of new therapies. Historically, these features of new drug candidates have been addressed at a late stage of development, mostly because the traditional methodology of studying acquired resistance (via serial dilutions) is slow and labour-intensive [5]. The recent introduction of continuous culturing methodologies combined with highly efficient genome sequencing is expected to change the traditional drug discovery paradigm. It will enable the includision of ‘resistibility’ profiling at an earlier stage of testing and triaging drug candidates, on a par with conventional efficacy and safety testing. More efficient methods or experimental evolution-based drug resistance studies are also required to accelerate the identification of molecular targets of novel phenotypically selected antibacterial agents with unknown mechanisms of action.

Realization of the advantages of automated (computer-controlled) continuous culturing over manual serial dilution methods of experimental evolution [6] was enabled by the recently enhanced accessibility of 3D printing technology and programmable microcontrollers (such as Arduino and Raspberry Pi). One of the rapidly emerging custom-built continuous culturing devices [7-10], termed ‘morbidostat’, is a modification of the chemostat in which, instead of nutrient supply, the culture density is controlled by selective drug pressure [11, 12]. This approach enables gradual genetic adaptation of evolving bacterial populations to higher concentrations of an antibiotic and generally recapitulates the knowledge of resistance mechanisms from previous studies [11, 13, 14].

Among recently reported continuous culturing devices, eVOLVER [9] and omnistat [7] provide technologically sophisticated multi-purpose platforms supporting a broad range of microbial evolutionary experiments in chemostat, turbidostat and morbidostat modes. The key feature of the eVOLVER platform is the capability for highly parallel processing of multiple samples over a variety of individually controlled experimental conditions. By the addition of specialized circuit boards and microfluidics, it can be customized to implement a morbidostat evolutionary mode. The omnistat platform features a simpler approach to automated liquid handling by bundling together a number of commercially available laboratory instruments. Our implementation of the morbidostat-based workflow (Fig. 1) included the engineering of a simple and robust device that uses inexpensive commercially available components, does not require sophisticated robotics, and allows Arduino-based programmable process control and data acquisition (Fig. 2). Our current morbidostat design supports parallel evolution in six 20 ml reactors (glass tubes). It has a minimal number of parts requiring low-resolution 3D printing, and it may be seamlessly replicated in any research laboratory based on the detailed hardware and software description provided on GitHub (https://github.com/sleyn/morbidostat_construction).

Fig. 1.

The experimental evolution of antibiotic resistance workflow includes four main stages. (1) Bacterial populations are evolving under increasing selective drug pressure in six parallel reactors of the morbidostat, a continuous culturing device. Time series of samples are collected from each reactor. (2) Total genomic DNA from each sample is sequenced with high coverage to obtain quantitative representation of low- and high-frequency sequence variants in bacterial populations at every time point. (3) Possible resistance mechanisms are deduced from the dynamics of observed mutations. (4) Follow-up experiments aimed to test bioinformatic predictions and elucidate genotype-to-phenotype associations include verification of identified mutations and measurement of the extent of acquired resistance (MIC fold change) in selected representative clones.

Fig. 2.

Morbidostat approach implementation. (a) Morbidostat control logic: cultures with OD LT, CT is decreased to 15 min. At OD > drug threshold (DT) the ‘medium with drug’ (pump 2) is added if the OD at the end of the current cycle is greater than at the end of the previous cycle. Otherwise the culture is diluted with medium without drug (pump 1). (b) General scheme of morbidostat. Two bottles: one containing medium without drug and the other containing medium with drug connected to two peristaltic pumps (pumps 1 and 2, respectively). Each pump is connected with six reactors. A valve on each tube controls media flow to individual reactors. Each reactor has three ports: (i) input for media and air, (ii) sampling port and (iii) waste line. The OD in each reactor is monitored by laser and photodiode. (c) Morbidostat device: side view with open box.

At the end of the first stage of the workflow, a selective outgrowth in the morbidostat (Fig. 1), increased antibiotic resistance is tested experimentally in collected microbial population samples (one–two samples per day until maximum resistance is reached, typically within 3 to 8 days). At the next stage, total genomic DNA isolated from these samples is used for the preparation of non-amplified genomic libraries and deep sequencing (up to 1000× genomic coverage). Primary data processing using an assembled computational pipeline (see the Methods section and https://github.com/sleyn/Triclosan_EE_paper) yields a list of identified mutations (SNVs, small indels, genomic rearrangements including ‘jumping’ of IS elements), along with their fractional representation (%) in each sample. Further filtration with chosen confidence criteria (thresholds) and downstream analysis of gene-by-gene distribution and dynamics of accumulation of mapped variants across all reactors and time points lead to tentative identification of genetic features associated with antibiotic resistance. The final stage of the workflow is aimed at experimental testing of inferred genotype-to-phenotype associations via mapping of mutations and determination of MIC values for selected individual clones.

Here we report the application of this workflow for the experimental evolution of triclosan (TCS) resistance in a model system of BW25113. TCS is an antimicrobial agent (biocide) commonly used in a variety of consumer products for dental care, cosmetics and soaps [15]. Extensive, indiscriminate use of TCS has led to its accumulation in the environment, in some soil samples reaching levels comparable to its MIC for and other susceptible bacterial pathogens (e.g. and ) [16]. An observed spontaneous acquisition of TCS resistance by pathogenic bacteria not only jeopardizes the utility of this biocide, but it may also pave the way toward their cross-resistance to some clinical antibiotics [17, 18].

The mechanism of action of TCS was established via the isolation of resistant clones with missense mutations in the gene fabI [19]. This essential gene encodes NADH-dependent enoyl-[acyl carrier protein] reductase, which is involved in one of the four elongation steps in fatty acid synthesis in most bacteria. TCS acts as a FabI inhibitor [19, 20], forming a stable and catalytically inactive TCS–NAD+–FabI complex [21-23]. In all eukaryotes the function of enoyl-[ACP] reductase is performed by a non-homologous domain of a multi-domain fatty acid synthase, which makes FabI attractive as a highly selective antimicrobial drug target. In some Gram-positive bacteria (such as ) the gene fabI is replaced by a non-homologous isofunctional gene fabK, which explains their intrinsic TCS resistance [19, 24]. In addition to some amino acid substitutions in the FabI target (most notably of Phe203 in ecFabI), an acquired TCS resistance was reported to result from overexpression of the fabI gene [25] and deregulation of efflux transporters [26]. However, according to some reports, TCS may also act on additional yet unknown targets at higher concentrations of the drug [27, 28].

Therefore, one goal of our study was to assess all possible mechanisms of acquired TCS resistance via experimental evolution across a broad range of TCS concentrations, from sub-MIC (<1 µm) to the highest concentration achievable in aqueous media (50 µm). The observed convergent evolution trajectories unambiguously pointed to the FabI enzyme being the sole protein target of TCS. We have identified a broad range of resistance-conferring amino acid substitutions in FabI (including all of the previously reported mutations) as the only tractable mechanism of robust TCS resistance in . All identified substitutions are localized in the vicinity of the FabI–TCS binding site and confer different levels of acquired resistance, ranging from ~3 to 16× MIC. This model study also points to the potential utility of the established morbidostat-based approach for uncovering resistance mechanisms and identifying targets for novel drug candidates with as yet unknown mechanisms of action.

Methods

Morbidostat design (Fig. 2)

We designed and assembled the first version of the morbidostat device used in this study based on the principles described by Toprak et al. [12]. A detailed description of a more advanced current version of hardware and accompanying software is available from GitHub (https://github.com/sleyn/morbidostat_construction). Briefly, the main components of the device (Fig. 2c) include: (i) 6×20 ml bioreactors – glass tubes with magnetic stir bars and cap assemblies containing air-tight needle ports for feeding dilution media, constant air flow, and displacement of excess liquid after dilution; (ii) glass tubes fixed by silicone rings in 3D-printed plastic housings containing a laser and a sensor diode for measuring turbidity; (iii) a mini air-pump, which provides constant aeration for growing cultures (through sterile filters) and enables liquid displacement from reactors (over a fixed level corresponding to a total volume of 20 ml); (iv) a thermoregulated heater and mini-fan, which control the temperature inside the morbidostat enclosure box; (v) two 2 l feed bottles controlled by 2 peristaltic pumps connected by sterile tubing to an assembly of 12 check valves (2 valves per reactor, each connected to 1 of the 2 pumps) providing the flow of media (with and without antibiotic); (vi) a 6-position magnetic stir plate, which enables mixing (upon dilution) and agitation. An Arduino microcontroller is programmed to control the main parameters of the run: (i) temperature; (ii) time between dilutions (cycle time, CT); (iii) selection of the volume delivered by pump 1 (medium without drug) or pump 2 (medium with drug) for each culture dilution, depending on the culture turbidity and growth rate in each tube. A user interface for parameter manipulation and real-time status display (including growth curves) is run on a PC using MegunoLink software (https://www.megunolink.com/).

The major differences between our implementation of morbidostat hardware and the originally published design [12] include (Fig. 2b):

using two pumps and a system of check valves to control the media flow instead of an array of multiple individual pumps, which allowed us to significantly simplify and increase robustness of the fluidics component;

a permanent sterile-filtered air flow (disrupted only for short periods of time during media addition and mixing), which allows us to maintain constant volume (20 ml) in the reactor and create an air gap between media and culture; the latter allows us to avoid cross-contamination between reactors and input lines.

The encoded morbidostat logic (Fig. 2a) includes the following major steps:

an optical density at the end of the current cycle, before dilution (OD1), is compared with two parameters: predefined drug threshold (DT) and optical density reached during the previous cycle (OD0);

if OD1 ≥DT and (OD1−OD0) ≥0, the dilution is made by drug-containing medium (pump 2);

otherwise drug-free medium (pump 1) is added.

An additional step was introduced on top of the original morbidostat logic [12]. It operates with an additional predefined parameter, the lower threshold (LT) (see Fig. 2a), instructing the system to skip periodic dilutions (typically every 10–20 min for rapidly growing cultures such as ) if OD1

Parental strain

We used the BW25113 ΔuxaC::kan strain (JW3063) from the KEIO collection [29] as a starting point for the evolution of TCS resistance. The choice of this strain (instead of the wild-type BW25113) with a transposon-replaced nonessential gene (hexuronate transporter UxaC) was dictated by two considerations: (i) to take advantage of transposon-encoded kanamycin (Kan) resistance to guarantee sterility, and (ii) to leverage the KEIO collection of gene deletion mutants with the same genetic background for follow-up validation experiments. A total genomic DNA extract from the unevolved culture used to inoculate all six reactors of the morbidostat on day 1 was sequenced with mean coverage of 530×.

Evolutionary runs (Fig. 3)

On the first of the four consecutive evolutionary runs (day 1), growth in all six reactors of the morbidostat was started using the same inoculum of the chosen parental strain [ BW25113 ΔuxaC::kan strain (JW3063)]. This inoculum was obtained by propagation of the primary glycerol stock from the KEIO collection [29] in Luria-Bertani (LB) medium with 50 µg ml−1 Kan and 2 % DMSO. The overnight culture was diluted 100-fold, grown to OD600 ~0.2 and diluted 10-fold with the same medium, which was also used as drug-free (pump 1) medium for culture dilutions in the morbidostat over the course of the entire study. A drug-feeding (pump 2) medium had the same composition with 10 µm TCS (during day 1) or 50 µm TCS (days 2–4). In the following consecutive runs (days 2–4), glycerol stocks made from samples collected at the end of the previous run were used to inoculate respective reactors. A total of four one-day (continuous 24 h) runs were performed. After each run, six culture samples were collected; reactors and all feedlines were sterilized prior to starting the next run. The main parameters used in each run were: dilution volume V=2 ml (1/10 of total volume); LT=0.15; DT=0.4; CT=15 min. In addition to preparing glycerol stocks that were used for inoculation, MIC analysis, and later isolation of individual clones, a portion of each collected sample was used to prepare and freeze cell pellets for genomic DNA extraction.

Whole-genome sequencing (WGS)

Genomic DNA from the starting culture and 24 samples of evolving bacterial populations in time series (6 reactors ×4 days) was extracted from cell pellets with the Sigma Aldrich GenElute Bacterial Genomic DNA kit (NA2110). Sequencing libraries were prepared with the Illumina TruSeq DNA PCR-Free High Throughput Library Prep kit (#20015963) according to the manufacturer’s instructions. DNA was sequenced with the Illumina NextSeq 500 using 2×150 reads to reach up to 500–1000-fold genomic coverage in each sample (see Table S1, available in the online version of this article, for data quality and statistics). Reads were cleaned from adapter sequences, poly-G sequences and were trimmed by Phred base quality score with the Trimmomatic [30] (with ‘SLIDINGWINDOW : 4 : 15 HEADCROP : 10 CROP : 145 MINLEN : 65′ parameters). Cleaned reads were mapped onto an BW25113 genome downloaded from the PATRIC database (genome ID 679895.18) [31] with the BWA aligner [32] (using the ‘-M’ parameter). Alignments were corrected with lofreq viterbi [33] (using ‘--keepflags’ parameter). Quality scores were recalibrated with GATK BaseRecalibrator [34]. Sequencing data in fastq format have been deposited at the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA), BioProject PRJNA472810.

Variant calling and downstream bioinformatics analysis

Variant calling was performed using lofreq [33] (with the ‘--call-indels’ parameter). The scripts are available in the https://github.com/sleyn/Triclosan_EE_paper GitHub folder. Variants in repeat regions were filtered out. The repeat regions were identified via self-alignment of the BW25113 genome using NUCmer [35] according to the ‘Identifying repeats’ section of the MUMmer manual (http://mummer.sourceforge.net/manual/) and filtered out. To eliminate sequencing biases manifest in unusual enrichment of low-frequency A→C and T→G variants, we applied additional filters to the lofreq output (see Methods S1). Structural variants were further analysed using the Breseq pipeline [36] in ‘population’ mode. To assess IS element rearrangements in populations, we have developed a new iJump tool (a description and the code are available on https://github.com/sleyn/ijump) [37]. The list of IS elements for tracking their rearrangements was created by using blastn of the reference genome against the ISfinder database [38]. The IS element blast hits were filtered to keep ones with e-values <10−30 and the highest bitscore among overlapping hits (see Table S2). We kept iJump insertion predictions that reached at least 1 % frequency in any reactor after summarizing frequencies estimated for all copies of IS elements. iJump predictions were inspected manually in BAM alignments with samtools [39] and Tablet [40]. To assess evolutionary dynamics, Muller plots were constructed for a subset of predominant mutations (selected from the complete list of observed confident structural variants provided in Table S3) and their combinations were tentatively deduced from synchronous patterns of variant frequencies. The ‘MullerPlot’ R package [41] was used for plot visualization. Initial plots were further adjusted based on the analysis of mutations in isolated individual clones.

Selection and characterization of individual clones

We used LB agar plates with a range of TCS concentrations (from 0 to 1.25 µm) to perform colony-forming unit (c.f.u.) analysis of 24 collected samples (6 reactors × 4 time points) of evolving populations in comparison with the unevolved starting culture of the BW25113 ΔuxaC::kan strain (see Fig. S2). From this initial screening, a total of 288 clones were selected for the initial assessment of MIC values by a microplate growth curve method using endpoint growth/no-growth measurements at 12 and 24 h (see below) (Table S6). Based on the diversity of MIC values and the representation of distinct sequence variants across all reactors and time points, 49 clones were ultimately selected for mapping of candidate mutations. This was achieved using PCR amplification and Sanger sequencing of selected regions in the genome (covering genes fabI, cadC, fimE, yncD, motB, rpoA) that contained the most prominent mutations. The obtained PCR products were sequenced by Eton Bioscience, Inc. (https://www.etonbio.com/). Insertions of IS elements in the yifE gene and the flhD–uspC intergenic region were tested by comparing the size of the PCR products from additionally selected clones with the wild-type PCR product. The 18 clones for IS element analysis were selected separately from previously described samples; these clones were sourced from reactor 3 at 24 h (n=8) and reactor 5 at 48 h (n=10).

MIC determination (IC50 and IC90)

The acquired TCS resistance of selected clones was measured using a serial dilution method. Clones were grown on LB/Kan (50 μg ml−1) agar plates. Colonies were resuspended in 1 ml of MHB/Kan solution and 5–10 µl of suspension was then incubated in 1 ml of MHB/Kan and 2 % DMSO in 96-well plates with a range of increasing TCS concentrations. For the lower range (closer to the IC90 of the unevolved parental strain=0.65 µm) the final concentrations of TCS were: 0, 0.3125, 0.625, 1.25, 2.5, 5, 10 and 20 µm. For 22 h, the OD of each sample at each concentration was measured at 600 nm wavelength by a BioTek ELx808 plate reader at 37 °C. For each well, the area under the growth curve (AUC) was calculated using the ‘growthcurver’ R package [42]. The IC90 values were calculated from the AUC–TCS concentration data (in independently grown duplicates) using the ‘drc’ R package and a curve-fitting procedure [43]. We used the LL4 model for curve-fitting. The same method was used to estimate changes in the apparent values of MIC in evolving populations (see Fig. S1), except that diluted glycerol stocks of collected samples, instead of individual colonies, were used to inoculate the microtitre plates.

Results

Populational dynamics of experimental evolution in the morbidostat

Pre-existing mutations in the unevolved population

Of the 10 pre-existing variants (listed in Table S4) revealed by comparison with the published genomic sequence of BW25113 (RefSeq assembly GCF_000750555.1), only 2 mutations affected the whole population (~100 % frequency): the G470S missense variant in the nnr gene encoding bifunctional NAD(P)H-hydrate repair enzyme and the disruption of lrhA gene (encoding a transcriptional regulator) mediated by IS2 mobile element.

Among other variants represented in the initial population (within the 1–20% frequency range), two variants displayed interesting dynamics during the course of experimental evolution (Fig. 4). First, a L65Q missense mutation in gene fimE encoding regulatory recombinase of the type I fimbriae formation expanded from 20 up to ~50 % of populations during day 1 in all six reactors. This expansion was followed by elimination of the FimE:L65Q variant at a later stage in all reactors except R4, where it remained at the level of ~10 %, and R5, where the main resistance-driving variant FabI:A21T emerged on FimE:L65Q background and extended to cover the entire bacterial population. Notably, the population representation of the FimE:L65Q variant showed good correlation with the frequency of the ON vs OFF position of the fim switch (Pearson correlation coefficient r=0.92, see Fig. S3). This observation suggests that FimE:L65Q likely has a lower recombinase activity, which controls the regulation of type I fimbriae expression via on-to-off recombination of the fim switch region [44]. This conjecture is consistent with the predicted localization of the L65 residue near the active site in a FimE 3D model (Fig. S4). Second, a low-frequency (~1 %) disruption of the yifE (maoP) gene by the IS2 mobile element also expanded up to 31 and 48 % of populations during day 1 in two out of six reactors. The remaining low-medium frequency pre-existing variants quickly disappeared from populations by day 2–3 in all or most of the reactors, along with the emergence and expansion of TCS-resistant FabI mutants.

Fig. 4.

Muller plots of evolution under triclosan stress in six reactors. Colours indicate various phenotypes. The labels on plot indicate emerged mutations. If a label begins with ‘+’, then a mutation is emerging, with the other mutation in the background.

Dynamics of TCS resistance and mutations in evolving bacterial populations

The dynamics and extent of TCS resistance acquired over the course of experimental evolution (days 1–4) were estimated using two approaches. First, the apparent ‘population-wide’ IC90 was assessed using a standard growth curve method without isolation of individual clones. The estimated IC90 showed a general trend, increasing from 2- to 4-fold in days 1–2 up to 10–20-fold by day 4 (Figs 3 and S1). Second, we performed c.f.u. analysis on agar plates over a range of TCS concentrations. This analysis revealed a gradual increase in the frequency of clones growing at TCS concentrations above the MIC of the unevolved strain (Fig. S2). Thus, in most reactors, between days 2 and 4, the extent of populational resistance increased from 12–15 % of clones growing at the highest TCS concentration, corresponding to 1× MIC, up to 40–70 % growing at 4–8× MIC. Up to 300 clones representing different time points and extents of evolved TCS resistance across all 6 reactors were picked for further analysis.

Fig. 3.

Experimental evolution of TCS resistance. (a) Six parallel reactors were inoculated from a single stock of BW25113 ΔuxaC::kan strain from he KEIO collection [29] and cultured in the morbidostat at 37 °C in four 24 h daily cycles (a–d). By the end of each daily cycle, the cultures were collected and split between cell pellets (for DNA extraction) and glycerol stocks used for inoculation of the next cycle. IC90 was measured in collected samples of evolving bacterial populations to confirm the evolution of TCS resistance. (b) Total genomic DNA was purified from 25 samples [unevolved inoculum + (6 reactors × 4 daily cycles)] and used for Illumina sequencing. After primary data processing and variant calling, a total of 60 variants were identified and analysed for evolutionary dynamics and mechanistic implications. The mutations in enoyl‐acyl carrier protein reductase fabI are shown as codons with mutated nucleotide and amino acid mutations. (c) For a set of 67 (49+18) isolated clones, 20 distinct variants (single, double and triple mutants) were mapped by PCR amplification/Sanger sequencing of selected loci. BW25113 (KEIO collection parent strain), BW25113 ΔuxaC::kan (parent strain for evolution runs) and 6 KEIO collection isolates (ΔcadC::kan,ΔcadA::kan,ΔmotB::kan,ΔyncD::kan,ΔlysP::kan,ΔcadB::kan) were added and IC90 values for TCS were determined using the growth curve method.

Following verification of gradual acquisition of drug resistance, total genomic DNA from all 24 collected populational samples was sequenced with a mean coverage of ~1000× (sd ~300×) (see Table S1). Upon primary data analysis, variant calling and filtration, we identified 60 distinct sequence variants that appeared in at least 1 reactor (see Table S3). Among them, were 36 missense and 7 synonymous mutations, 2 frameshifts, and 3 intergenic region mutations and small indels. In addition, 12 IS element insertion events were mapped using a newly developed iJump programme.

A tentative reconstruction of the dynamics of the appearance and disappearance of the main sequence variants created from deconvolution of populational sequencing data (Table S3), illustrated here by Muller plots (Fig. 4), revealed some notable trends. At an early stage of experimental evolution (days 1–2, up to 48 h), we only see a limited accumulation of new mutations. As already mentioned, the predominant adaptation processes included expansion within day 1 of the two pre-existing variants, FimE:L65Q (in all reactors) and yifE:IS2 (in reactors 3 and 5). While in most cases, these variants were outcompeted and eliminated from populations by day 3, in one reactor (R5), the FimE:L65Q variant expanded to 100 % and was sustained over the entire course of evolution due to the emergence of additional FabI mutations on the FimE:L65Q background.

Several low-frequency missense mutations (1–7 %) in genes cadC, lysP and motB that appeared within day 1 in at least one reactor (up to four reactors in the case of the CadC:M322T variant) either disappear by day 2 or combine with one of the emerging FabI mutations and expand at a later stage. In most reactors at an early stage (days 1–2), FabI mutations do not emerge at all (as in R1 and R5) or accumulate at a relatively low level (≤30 %, as in R2, R3 and R4). However, in R6, a FabI:F203L variant reaches ~80 % of population at day 2 as a triple mutant combined with missense mutations in CadC and YncD (see Fig. 4, Table S3). The observed faster accumulation of FabI in R6 is consistent with more rapid acquisition of robust TCS resistance at an earlier stage (day 2) of evolution (see Figs S1 and S2). The initial tentative assignment of double (and some triple) mutants that become predominant at a later stage (starting day 3) in most reactors (except R1), was based on parsimonious analysis of parallel trends in the frequency of mutations (from Table S3). However, most of the assignments shown in Fig. 4 were later verified by sequencing of selected individual clones (see Table 1).

Table 1.

IC90 and 95% confidence Interval (CI) values for clones isolated in this study (bold) and selected KO strains from the KEIO collection (gene ID in parentheses). The variants are highlighted by distinct positions of FabI mutations and sorted in ascending order by IC90 values. The colour in the second column shows the IC90 gradient from low (green) to high (red)

Evolved variants and KO strains	IC90 (μm)	CI	Fold change
I. Evolved FabI variants
FabI:T194S+CadC:M322T	1.99	1.79–2.18	3
FabI:T194S	2.23	1.74–2.72	3
FabI:I192F	3.52	3.25–3.78	5
FabI:A21T+CadC:M322T	3.56	3.36–3.77	5
FabI:A21T+FimE:L65Q	4.1	3.70–4.51	6
FabI:A21T	4.12	3.76–4.48	6
FabI:I192F+CadC:A220T	4.52	4.31–4.74	7
FabI:F203L	4.85	4.48–5.22	7
FabI:G93S	5.63	5.23–6.03	8
FabI:F203V+CadC:D506E+MotB:L192F	6.17	5.63–6.70	9
FabI:F203L+CadC:M322T+YncD:N319Y	6.72	6.23–7.21	10
FabI:F203C	7.37	6.84–7.90	11
FabI:M159T+FimE:L65Q	8.65	7.83–9.48	13
FabI:F203V	9.48	8.58–10.39	14
FabI:F203S	10.65	9.96–11.35	16
II. Other evolved variants
yifE:IS2*	0.63	0.58–0.69	0.9
us_flhD<>us_uspC:IS5†	0.70	0.62–0.77	1.1
III. Unevolved parental strain
uxaC-KO (JW3063)‡	0.65	0.62–0.67
uxaC-KO (JW3063) +FimE:L65Q‡	0.7	0.63–0.78
WT	0.85	0.75–0.96
IV. KO strains from KEIO collection
cadC-KO (JW4094)	0.6	0.58–0.61	0.9
cadA-KO (JW4092)	0.73	0.70–0.75	1.1
motB-KO (JW1878)	0.8	0.78–0.83	1.2
yncD-KO (JW1446)	0.91	0.87–0.94	1.4
lysP-KO (JW2143)	0.98	0.95–1.01	1.5
cadB-KO (JW4093)	1.01	0.95–1.06	1.5

*Evolved clone with IS2 insertion in yifE gene of unknown function.

†Evolved clone with IS5 insertion in the intergenic region between genes flhD and uspC.

‡Two clones with SNP in gene fimE were isolated from the inoculum of the unevolved parental strain.

§Wild-type BW25113 strain (ATCC).

Importantly, the dynamics of the appearance and expansion of double mutants allow us to confidently reconstruct some of the evolutionary trajectories. These dynamics also permit us to distinguish mutations that could only have emerged in the course of evolution, as opposed to those resulting from selective propagation of pre-existing micro-SNPs. Thus, in addition to the example of R6 outlined above, another mutation, FabI:A21T, emerged in R4 at day 2 in the background of a CadC:M322T variant (initially detected in day 1 without any FabI mutations). Starting at day 2, both mutations expanded in parallel from ~10 % ion day 2 up to >50 % on day 4 (see Fig. 4, Table S3). The inferred co-occurrence of both mutations in one clone was later verified directly (Table 1).]

Nevertheless, the most important and obvious trend is the evolution of robust TCS resistance via a range of alternative amino acid substitutions in the known TCS target, the enoyl-[ACP]-reductase FabI target, as a single general mechanism. Despite different observed dynamics and a diverse array of specific mutations and combinations thereof, 97–99 % of bacterial populations in all reactors at day 4 contained at least one FabI variant (most commonly F203L or F203C). The emergence of a FabI:Q40L variant in the background of the dominant FabI:F203L mutant on day 4 in R3 was inferred from population sequencing data (Fig. 4), but the respective clone was not isolated due to low frequency (~6 %).

Mechanistically relevant variants

Overall, we identified 12 missense variants in 8 positions in the gene fabI that are distinct at a nucleotide level, resulting in 10 amino acid substitutions. All of the respective amino acid residues are located in the vicinity of the TCS-binding site in the known 3D structure of the FabI-NADH-TCS complex (Fig. 5) [23]. Five out of 12 variants are substitutions at FabI:F203 (2 versions of F203L; F203C; F203V; and F203S). Notably, these variants are predominant at day 4 in nearly all reactors except R5, which is dominated by FabI:A21T. However, the A21T mutation of FabI is the second most frequent among all FabI variants and only reaches 99 % in the reactor lacking any FabI:F203 variants (R5). This observation is consistent with a superior intrinsic resistance/fitness of FabI:F203 variants. Of note, a double mutant containing a RpoA:L289F variant in the FabI:A21T background emerges and expands up to ~1/3 of the population in R5 at day 4 (Fig. 4). This variant was not characterized at the clonal level, but since no other mutations in RpoA (alpha subunit of DNA-directed RNA polymerase) or any other genes involved in DNA/RNA metabolism were detected, its impact on TCS resistance/fitness appears unlikely.

Fig. 5.

Three-dimensional structure of FabI (PDB ID: 1QSG [23]) with NAD (yellow) and TCS (orange). Amino acid residues affected by mutations identified in this study are coloured red. For mutations where the IC90 is known it is shown by extension lines in µm.

Beyond FabI target-based adaptation, only two low-frequency events disrupting the acrR gene, the known transcriptional repressor of the multidrug efflux transporter AcrAB [45], have a clear link to a possible mechanism of TCS resistance. An acrR frameshift mutant emerged at day 4 and reached 7 % in the population of R3 on the background of the predominant FabI:F203L variant. The second variant with IS5 mobile element insertion in acrR gene reached maximum frequency (6 %) at day 3 in R1, but was outcompeted by the FabI:F203L single mutant by day 4. None of these variants were characterized at the clonal level due to frequency, and their actual contribution to TCS resistance remains unknown.

Other evolved variants

As noted above, mutations in other genes reached appreciable levels either at an early stage and/or when combined with one of the FabI variants at a later stage of experimental evolution. Known functions of these genes do not point to any obvious links with TCS resistance. However, the observed multiple missense mutations in gene cadC point to a possible relevance of this integral membrane protein, a pH-sensing transcriptional activator of cadaverine pathway CadAB, for adaptation to a relatively low TCS pressure or another type of stress at the early stage of evolution. Notably, all five missense variants of CadC (A220T, M322T, R333S, N466Y and D506E) are found in the C-terminal signal-sensing domain of CadC 3D structure (PDB ID: 3LY7) [46]. Another observed early stage variant, LysP:G80S, observed in three reactors on day 1 (at ~3.5 % frequency), affects lysine permease LysP, which is known to interact and negatively regulate CadC. Due to low frequencies of respective single mutants, none of these variants were characterized at the clonal level.

Another interesting variant, insertion of IS5 element upstream of the flhDC operon, emerged and peaked in two reactors on day 3 (30 % in R4 and ~70 % in R5). This operon encoding transcriptional activators of the flagella biosynthesis pathway was disrupted by IS5 insertion at two positions: (i) –95/–99 interfering with LrhA (–89/–129) and OmpR (–101/–163) binding sites, and (ii) –315/–318 upstream of them (relative to the transcription start site according to Baker et al. [47]). These positions are known hotspots for IS element insertions leading to upregulation of flagellar and mobility genes in [47]. Notably, previously discussed IS2 insertions in the gene yifE also have a tentative association with the flagella pathway, as this gene is located upstream of hdfR, the regulator of flagellar master operon flhDC [48]. Analysis of isolated clones representing both variants (see below) revealed no appreciable impact on TCS IC90 (see Table 1), suggesting that a yet unknown adaptive role of these variants (if any) is not directly related to drug resistance.

Verification and characterization of individual evolved clones

A total of 49 clones representing the maximum diversity of TCS resistance (determined from preliminary screening of 288 clones, Fig. S2, Table S6) and anticipated mutations (deduced from evolutionary dynamics, Table S3) were selected for detailed characterization. This included accurate determination of TCS IC90 and mapping of potential driver mutations by PCR amplification of the respective segments of the cadC, fabI, yncD, motB, fimE and rpoA genes. These results obtained for 15 evolved clones with mapped distinct mutations and 2 distinct sub-clones isolated from the unevolved population of BW25113 ΔuxaC::kan strain (with and without FimE:L65Q variant) are provided in Table 1. In addition, two clones representing IS element insertions in the yifE gene and the flhD–uspC intergenic region were isolated from the respective samples, verified by PCR/Sanger sequencing and tested for TCS resistance. A panel of characterized clones covered nearly all of the FabI single mutant variants and some of the observed double and triple mutants. The observed IC90 values vary substantially in a range from ~2 up to >10 µm over the entire panel of 15 FabI variants (section I in Table 1). This corresponds to a fold change from 3 to 16× as compared to the averaged IC90 of the parental strain sub-clones of BW25113 ΔuxaC::kan with and without the FimE:L65Q variant (section III in Table 1). These variations depend on the position and specific amino acid substitution (as in the case of multiple FabI:F203 variants) but, with a few exceptions, show almost no dependence on the presence or absence of mutations in other genes. The observed variations in TCS resistance are generally consistent with the outcome of convergent evolution, which included the competition of distinct FabI variants within the same reactors. As already emphasized, most day 4 populations are dominated by the FabI:F203-driven variants that show the highest level of TCS resistance. Not surprisingly, at day 4 in R6, the FabI:F203S single mutant, which has the highest IC90(10.7 µm), ultimately outcompeted a triple mutant (FabI:F203L+CadC:M322T+YncD:N319Y) with a somewhat lower IC90 (6.7 µm).

Although CadC single mutant variants were not isolated due to their low frequency, we used a series of strains from the KEIO collection [29] to assess the effect of individual gene knockouts in the cad pathway (genes cadC, cadA, cadB, lysP) on TCS resistance (Table 1). CadC is a transcriptional activator of lysine/cadaverine antiporter cadB and lysine decarboxylase cadA [49], and lysine transporter lysP is a posttranslational negative regulator of CadC [50]. Therefore, deletion of cadC might be expected to have a negative effect, similar to that of the deletion of cadA or cadB, and opposite to the effect of lysP deletion. However, the observed IC90 values were either nearly the same (for CadC-KO and CadA-KO) or only ~1.5-fold higher (for CadB-KO, LysP-KO) than for the isogenic control strain (UxaC-KO).

Cross-resistance between biocides and antibiotics is a concern for the overuse of TCS in consumer products, and exposure to biocides was found to decrease the susceptibility of to multiple antibiotics, including broad-spectrum quinolones [51]. Since their presence did not increase TCS resistance directly, a series of CadC variants were examined for possible cross-resistance effects of such mutations. Five single FabI-mutant clones and seven double CadC-FabI mutant clones (with shared FabI mutations) were tested for resistance to various antibiotics (ampicillin, chloramphenicol, ciprofloxacin, erythromycin and tetracycline). All clones had the same MIC for each tested drug, exhibiting no cross-resistance trends with the presence or absence of CadC mutations.

Discussion

In this study we have evaluated the morbidostat-based continuous culture method to study the evolution of antibiotic resistance. For this purpose, we have established an experimental evolution workflow (Fig. 1), which included engineering a morbidostat device based on a published conceptual design (Fig. 2) [12], and applied the workflow and device to study the evolution of TCS resistance in (Fig. 3). This analysis revealed a notable diversity of emerging mutant variants expanding and contracting in six parallel bacterial cultures and showing distinct evolutionary trajectories. All these trajectories converged within 4 days into a limited number of solutions (Fig. 4) sharing one general mechanism of robust TCS resistance (up to 16-fold increase in IC90) driven by a range of adaptive amino acid substitutions in the vicinity of NAD/TCS-binding site of the known target enzyme enoyl-[ACP]-reductase FabI (Fig. 5). The 8 detected amino acid residues comprising a total of 12 distinct substitution variants included all 3 FabI positions (G93, M159 and F203) previously reported to confer robust TCS resistance [19, 24]. These three positions, most notably F203, featuring four distinct substitutions, are associated with FabI variants dominating bacterial populations at the latest stage of evolution and showing the highest extent of TCS resistance measured in isolated clones (Table 1). Other FabI variants identified only in this study (affecting A21, Q40, I192, T194 and A197 residues) were generally less competitive (except for a special case of A21T) and showed lower IC90 for TCS. This observation illustrates an advantage of the dynamic genomic profiling-based approach (over a more conventional clone isolation) for the detection of transient variants with lower drug resistance and/or fitness. Remarkably, most of these ‘low-ranking’ residues are located at the NAD side of the FabI–NAD–TCS complex (Fig. 5), and some of them (I192 and T194) were shown to form hydrogen bonds with NAD [23]. Therefore, it is tempting to hypothesize that the respective mutations are more likely to confer TCS resistance via shifting NAD position in the FabI active site, which may be a less efficient strategy than a direct interference with TCS binding. Indeed, the three ‘high-ranking’ residues (G93, M159 and F203) are located on the TCS side of the FabI–NAD–TCS complex and are known to form hydrophobic interactions with the TCS molecule [23].

Beyond FabI variants, we observed only two events potentially contributing to TCS resistance via upregulation of the acrAB efflux pump operon due to inactivation of the AcrR repressor. Both reached only a low frequency (6–7 %) at the later stage (days 3–4) of experimental evolution. Although none of the respective clones were characterized, in previous studies, acrAB overexpression was shown to provide only a twofold increase in TCS-resistance [26]. These observations generally suggest that AcrAB-driven efflux upregulation provides a degree of resistance or fitness that is less competitive under selective conditions in the morbidostat.

Of the early stage events (Fig. 4, Table S3), several low-frequency missense mutations in the gene cadC (and one in lysP) appeared to be the most promising candidates for an alternative mechanism of TCS resistance. This conjecture was based on the observed variety of independent, distinct and exclusively missense mutations that are all localized in the C-terminal periplasmic domain of the pH-sensing CadC transcriptional activator of the cadBA operon. Although the exact function and mechanism of action of CadC have not been fully elucidated, it has been implicated in acid stress response, which could provide a mechanistic basis for the adaptation to TCS-induced stress. However, a comparative characterization of the five pairs of isolated FabI single mutant variants with respective (FabI +CadC) double mutants failed to provide conclusive support for this hypothesis (see Table 1, Fig. S5). Therefore, at this point a mechanistic role of the observed cadC mutations in the evolution of TCS resistance (if any) remains unclear. Interestingly, a point mutation, CadC-Y504H, was implicated in the increased resistance of avian pathogenic to some plant-derived antibacterial agents, such as carvacrol [52].

The observed early-stage expansion of pre-existing variants (FimE:L65Q+YfiE:IS2), as well as the emergence of the (flhD–uspC:IS5) variant, point to a possible contribution of fimbriae and flagella biogenesis to TCS adaptation. This conjecture is consistent with previous reports on (i) activation of flagellar assembly genes in response to TCS [53, 54] and (ii) the contribution of fimbriae biosynthesis (genes fimABCDEI) to survival in the presence of TCS revealed by transposon mutagenesis [55]. However, none of the isolated clones representing these variants showed any improvement of TCS IC90 (Table 1). Likewise, a double mutant (FabI:A21T+FimE:L65Q) variant showed the same level of TCS resistance as a respective single mutant FabI:A21T. Overall, the obtained clonal data did not provide support for the involvement of flagella or fimbriae pathways in the evolution of TCS resistance. The observed transient expansion of respective variants may be explained by the general adaptation to growth in the morbidostat or by compensation for the pre-existing IS2-disruption of the gene lrhA (Table S4), a master regulator of motility, chemotaxis, fimbriae and flagella biosynthesis genes [56].

In summary, in this study we have successfully used the established morbidostat-based workflow to assess the dynamics of experimental evolution and identify and verify a range of mutations driving acquired TCS resistance in . Despite the diversity of evolutionary trajectories and mutations across six parallel reactors, we have observed parallel evolution toward a single mechanism of robust TCS resistance (up to 16× increase of IC90). This level of resistance is achieved by the substitution of at least one of several amino acid residues in the active site of FabI, the essential enzyme in fatty acid synthesis and the known TCS target. The established approach allowed us to detect all previously reported FabI mutations conferring TCS resistance and identify a broader set of previously undetected ones.

Overall, the obtained results provide a proof of concept and illustrate the potential advantages of the developed methodology for target identification (or verification) and the elucidation of yet unknown mechanisms of action and resistance for established antibiotics and novel drug candidates. The prevalence of FabI mutations is in line with recent studies of bacterial isolates (human- and soil-associated) and metagenomes, where it was shown that mutations in FabI or substitutions to resistant analogues are the most abundant mechanisms for TCS resistance (see Table S7) [24, 57].

Click here for additional data file.

Click here for additional data file.