Resistance variation and bacterial interactions shape adaptation of a genetically diverse pathogen population to antibiotic therapy
Aditi Batra, Leif Tueffers, Kira Haas, Tabea Loeblein, João Botelho, Michael Habig, Daniel Schuetz, Gabija Sakalyte, Florian Buchholz, Ernesto Berríos-Caro, Hildegard Uecker, Daniel Unterweger, Hinrich Schulenburg
TL;DR
This study explores how antibiotic resistance evolves in mixed bacterial populations, showing that genetic diversity and interactions between strains influence resistance development.
Contribution
The study reveals that AMR evolution in polymicrobial infections depends on strain-specific resistance, ecological interactions, and mutation rates.
Findings
AMR evolution in Pseudomonas aeruginosa multistrain communities is influenced by strain-specific resistance profiles and ecological interactions.
Mutation rates for resistance affect the likelihood of new resistance emerging during antibiotic treatment.
Strain–strain interactions and genetic diversity play a central role in adaptation to antibiotics.
Abstract
Antimicrobial resistance (AMR) poses a major threat to global human health. The emergence and spread of AMR is usually studied for single pathogen lineages. Therefore, we currently have only limited knowledge on the causes and dynamics of resistance evolution in polymicrobial or multistrain infections that involve different pathogen species or strains, respectively, even though these kinds of infections are widespread. To address these current knowledge gaps, we here used the opportunistic human pathogen Pseudomonas aeruginosa as a model to investigate how AMR evolves in populations with different genetically distinct strains (multistrain communities). By using controlled evolution experiments, extensive phenotyping and genome sequence analysis, we demonstrate that the response to antibiotic selection is shaped by a combination of strain-specific resistance profiles, ecological…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6- —Deutsche Forschungsgemeinschaft10.13039/501100001659
- —Excellence Cluster Precision Medicine in Chronic Inflammation
- —Walter Benjamin
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsEvolution and Genetic Dynamics · Antibiotic Use and Resistance · Antibiotic Resistance in Bacteria
Introduction
Antimicrobial resistance (AMR) is an eco-evolutionary problem at its core [1]. Bacterial adaptation in response to antibiotic usage has resulted in widespread AMR, rendering many drugs ineffective [2]. As a consequence, we find ourselves in the middle of an AMR crisis, with approximately 4.95 million deaths worldwide associated with AMR in 2019 [3, 4].
The key to countering AMR evolution is to understand how it happens. AMR evolution is usually a consequence of the selection caused by antibiotic therapy [5] and can emerge within a few days of treatment [6–8]. AMR may be shaped by the number of infecting strains and species and thus be driven by a combination of ecological and evolutionary processes. Single strains dominate some infections, such as those related to sepsis [9] or uncomplicated urinary tract infections [10]. In these infections, available resistance and pathogen evolvability shape AMR. Many other diseases, however, can be polymicrobial. Infections in the abdominal cavity such as acute cholangitis [11] or acute pancreatitis [12] are frequently caused by subsets of the gut microbiome. The lungs of cystic fibrosis (CF) patients are often infected with multiple strains of Pseudomonas aeruginosa [13–16]. Based on this background, we here ask how AMR evolves in such multistrain infections?
To date, most work on AMR in multistrain infections is based on observational data collected from patients [13, 17, 18], for which the inference of cause–effect relationships is often a challenge. The few experimental studies on the topic focused on interactions between different bacterial species [19–22], thus overlooking the variation commonly encountered within single species [13, 14–16]. AMR in such intraspecific infections could be favoured by pre-existing AMR [23] or horizontal gene transfer (HGT) [15]. Intraspecific variation may further associate with beneficial, neutral, or competitive interactions between strains [24], which in turn can affect AMR spread, contingent on the exact nature of the interaction [25, 26, 27]. Metapopulation structure is an additional factor, particularly relevant for lung infections [28], that can lead to independently evolving subpopulations [29], facilitating AMR emergence and spread [30]. Overall, standing genetic variation, microbial interactions, HGT, and metapopulation dynamics are likely critical factors in shaping AMR in infecting microbe populations, yet, to date, their exact role in AMR evolution is largely unexplored.
To address knowledge gaps, we performed a proof-of-concept study that utilizes an experimental approach to assess the extent to which pathogen evolution is influenced by (i) standing genetic variation in AMR, (ii) microbial interactions, (iii) HGT, and (iv) metapopulation structure. Our unit of selection was a genetically diverse population (gdPop) that contained a mixture of 12 strains of the human opportunistic pathogen P. aeruginosa as a model. We first characterized variation in AMR and microbial interactions within the gdPop. Thereafter, the gdPop was subjected to an evolution experiment with different antibiotic treatments imposing different selective pressures and two different transfer protocols, simulating two levels of metapopulation structuring. We analysed growth dynamics and strain diversity during experimental evolution and characterized evolved AMR, genome sequence changes, and HGT at the end of the experiment. The insights from these measurements were tested and validated with an additional, independently performed evolution experiment and, furthermore, fluctuation assays to infer rates of AMR emergence.
Materials and Methods
Strains
The gdPop used for experiments consisted of 12 strains of P. aeruginosa. These were taken from the mPact strain panel representing the entire genomic diversity of the species [31, 32]. Each of the 12 strains was grown separately in 10 ml Luria Bertani (LB) medium for 20 h at 37°C and 150 rpm following which they were mixed in a 1:1 ratio (by volume - v/v) to obtain the final mixed population. No statistically significant difference in the CFU/ml was observed at the end of 20 h between the different strains (Kruskal–Wallis test, χ2 = 19, P = .456 [32]) This mixture was then used for experiments.
Media and antibiotics
All experiments with the gdPop were conducted in either LB (Carl Roth Art. No. X964.2) or minimal M9 medium as broth or combined with agar (1.5%, Carl Roth Art. No. 5210.5). M9 medium was supplemented with glucose (2 g/l, Carl Roth Art. No. 6780.1), citrate (0.58 g/l, Carl Roth Art. No. 4088.3), and casamino acids (1 g/l, Carl Roth Art. No. AE41.1). Antibiotics were added when required. Antibiotics used included the aminoglycoside gentamicin (GEN, Carl Roth, Order No. HN09.1) and a combination of the beta-lactam piperacillin (PIP, Sigma Aldrich, Code: P8396-1G) and the beta-lactamase inhibitor tazobactam (TAZ, Sigma Aldrich, Code: T2820-10MG). The inhibitor was included as some of the strains possess beta-lactamases [32] most of whose activity would be blocked by tazobactam. Piperacillin and tazobactam (PTZ) is the recommended first-line treatment for severe P. aeruginosa infections. GEN is used for severe infections such as nosocomial ventilator-associated pneumonia, malignant otitis externa, and in topical applications. PTZ were combined in an 8:1 ratio. Bacteria were incubated at 37°C.
Experiments with conditioned media
To determine contact-independent pairwise interactions, we cultured the 12 strains separately in LB and M9 media overnight at 37°C. After incubation, cultures grown in the M9 medium were filtered through 0.22 μm filters to generate bacteria-free supernatant. Strains grown in LB were subcultured until the mid-exponential phase and washed from the residual medium. We resuspended washed bacteria in the filtered M9 medium using the matrix-like format in a 96-well plate and normalized the initial OD to OD_600_ = 0.001. Cultures were incubated at 37°C for 62 h with shaking (162 rpm). Measurements of OD_600_ were taken every 15 min with a microplate reader (Tecan, Spark). We performed three independent biological repetitions for each combination of strain and conditioned media. Area under the curve (AUC) measures were calculated using the R package growthcurver. Strain interactions were determined using the formula:
Effect of Strain B on Strain A = AUC (Strain A cells under Strain B supernatant) – AUC (Strain A cells under Strain A supernatant)
AUC values for the effect of Strain B supernatant on other strains’ biomass (Focal strain - > Others) were summed up to obtain the cumulative effect of Strain B’s supernatant. AUC values for the effect of other strains’ supernatant on Strain B’s biomass (Others - > Focal strain) were summed to obtain the cumulative effect on Strain B’s biomass.
Minimum inhibitory concentration measurement using MIC test strips
Minimum inhibitory concentration (MIC) test strips (Liofilchem) [33] were used to assess variation in antibiotic resistance among individual strains. The individual strains were grown in 10 ml M9 for 18 h. Each strain was then diluted to an OD_600_ of 0.08 and spread onto an M9 plate using a cotton swab. MIC test strips (Liofilchem) were placed onto the plate, and the plate was incubated for 24 h, after which the MIC was read at the intersection of the zone of inhibition and the test strip. For testing the MIC of the gdPop, the 12 strains were grown individually in 10 ml LB for 20 h and then mixed in equal proportions. This mixture was then diluted to an OD600 of 0.08 and spread onto a M9 plate using a cotton swab, followed by MIC inference as above using MIC test strips (Liofilchem). For PTZ, we observed bacterial growth in the shape of wings above the intersection of the zone of inhibition and the test strip (Fig. S1B). This is likely due to the paradoxical effect [34]. It is important to note that this effect was not observed on PTZ when individual strains were tested. The MIC was read where the top edge of the wings intersected the test strip. The measured MIC reflects the antibiotic concentration needed to inhibit visible growth of bacteria and accounts for any pre-existing resistances due to mutations/enzymes. To measure MIC of the evolved populations from day 14, we thawed the frozen material from this timepoint, transferred 50 μl of the sample into 10 ml LB, followed by cultivation for 20 h at 37°C and 150 rpm. After incubation, the cultures were diluted to an OD_600_ of 0.08 and then prepared for strip-based MIC inference as above. Change in MIC of the evolved population was calculated as MIC_evolved population_ - MIC_ancestor_. Not all frozen populations could be recovered for MIC testing.
OD-based dose–response curves for calculating resistance changes of the evolved populations (broth microdilution)
The OD-based broth microdilution approach was used to assess changes in antibiotic resistance among ancestral and evolved populations. For the ancestral population, the 12 strains were grown individually in 10 ml LB for 20 h after which they were mixed in equal proportions. A total of 2 μl (2%) of this mixture was inoculated into the wells of a 96-well plate (100 μl total volume, Greiner Ref. 655 161) containing a range of concentrations of antibiotics covering the MICs of the individual strains. The plate was incubated at 37°C and 750 rpm for 24 h after which the OD_600_ was read in a BioTek (Epoch 2) plate reader. The area under the resulting dose–response curve was calculated using the drc package in R and was used as a measure of the antimicrobial susceptibility of the ancestral population. For measuring the antimicrobial susceptibility of the evolved populations, the populations frozen on day 14 were thawed and 50 μl was inoculated into 10 ml LB and then subjected to the same procedure described above. Change in antimicrobial susceptibility of the evolved population was calculated as antimicrobial susceptibility evolved population - antimicrobial susceptibility ancestor.
CFU-based dose–response curves for standardizing selection strength of antibiotic treatments (broth microdilution)
The CFU-based broth microdilution approach was used to infer the antibiotic concentrations that reduces cell number by three orders of magnitude compared to the corresponding no-drug control (Inhibitory concentration of 99.9% or IC99.9) and, thus, to determine the concentration to be used for imposing comparable selection in both evolution experiments, irrespective of the treatment conditions, antibiotic applied (both evolution experiments) or whether single strains or the gdPop were used (only second evolution experiment). Individual strains were grown overnight for 20 h in LB and then mixed in a 1:1 (v/v) ratio. A total of 40 μl (2%) of this mixture was inoculated into the wells of 12 well plates (Greiner Ref. 665180) containing a range of antibiotic concentrations chosen based on the MIC data of the individual strains. No-drug controls and blanks were included. All treatments were randomized on the plate. Total volume in each well was 2000 μl. Plates were incubated at 37°C and 900 rpm for 24 h. At the end of incubation, the plates were sonicated (Bandelin Sonorex) for 1 min to break any clumps of bacteria, and these were spotted onto plates to determine CFU/ml. This was done for GEN, PTZ, and their combination (Fig. S1A). For the combination DRC, each concentration tested contained an equal mixture of the two antibiotics. CFU-based dose–response curves were repeated for the second evolution experiment in the new experimental settings required for the experiment and again including the mixed gdPop as well as the individual strains H05 and H18 (Fig. S4).
Strain-specific polymerase chain reaction and estimation of strain diversity over time
To detect which strains survived postantibiotic treatment, strain-specific polymerase chain reactions (PCRs) were carried out. Each strain was uniquely identified by a primer pair (for a list of primer pairs, see Table S11). Samples were taken at concentrations that reduced the cell number by three orders of magnitude and subjected to PCRs. The reaction mix included Dream Taq buffer (1X, Thermo Fischer B71), dNTPs (0.25 mM, Thermo Fischer R0181), forward and reverse primers (0.5 μM each, Eurofins), Dream Taq (1.5 U, Thermo Fischer EP0702), and template DNA (10–100 ng). The mix was subjected to an initial denaturation of 45 seconds at 95°C, followed by 30 cycles of annealing for 30 seconds at 55°C and extension for 30 seconds at 72°C, and then a final extension for 10 min at 72°C (SensoQuest labcycler). To estimate strain diversity over time in the evolution experiment, strain-specific PCRs for each of the 12 strains were performed on 3 independent replicate populations from transfers 4, 7, 11, and 13. For the populations from transfer 13, DNA for the PCRs was isolated directly from the evolution experiment. For the rest, the frozen populations had to be regrown in LB for DNA isolation for the PCRs. The replicate populations were chosen randomly and from all 6 treatments. Once chosen, the same populations were tested at each transfer. Shannon index values were calculated from the strain occurrence at these transfers for every combination of treatment and meta-treatment.
Evolution experiment 1
The evolution experiment consisted of six treatments nested within two meta-treatments. The six treatments included three evolution-informed treatments (one combination and two sequential, each starting with a different antibiotic) and three controls (No drug, monotherapy with each antibiotic). Each treatment was initiated with eight independent biological replicate populations. Treatments were fully randomized on each plate. Each of the 12 strains was grown for 20 h in LB and then mixed in equal proportions (v/v). From independent mixtures, bacteria were inoculated into 12 well (Greiner Ref. 665180) antibiotic plates containing 2000 μl of medium. They were then grown for 24 h at 37°C and 900 rpm (Heidolph titramax 101), following which the plates were sonicated to break any bacterial clumps. 40 μl (2% v/v) was then transferred to another plate containing fresh medium. For the No-mixing meta-treatment, we ensured a 1-to-1 correspondence between wells during transfer. The “mixing” meta-treatment was introduced to limit strain loss due to genetic drift and simulate metapopulation structuring. For this, all eight replicates within a treatment were mixed in one 50 ml tube, and this mixture was used to inoculate all wells of this treatment for the next growth season (thereby the eight replicate populations of the mixing treatment became linked). Each of the meta-treatments had 6 treatments, resulting in a total of 12 treatments and 96 bacterial populations. Growth of the bacteria was followed by measuring the endpoint OD_600_ over time. The experiment lasted a total of 15 days with 14 days of antibiotic exposure followed by 1 day of growth in drug free medium to assess population extinction. Populations were frozen in DMSO on days 4, 7, 11, and 14. Populations that grew after the removal of antibiotics were frozen on day 15.
Calculation of population extinction and adaptation parameters
Both population extinction and adaptation of the surviving bacterial populations to the treatment were quantified at the end of the evolution experiment. A population was identified to be extinct if in the last growth season without drug its OD_600_ was less than 0.05. Two parameters were used for quantifying the pattern of adaptation. Time to adaptation is the time taken (in number of growth seasons) for the bacterial population under treatment to attain the same level of growth as those from the simultaneously evolving no-drug control. Robustness of adaptation is the number of growth seasons the bacterial population under treatment is able to maintain the growth level of the simultaneously evolving no-drug control (i.e. its OD_600_ was identified to be within the mean ± 2× standard deviation of the no-drug control in that growth season), after this level was first reached. This is a measure of how well adaptation was maintained once acquired.
Determination of strain abundance, horizontal gene transfer, and sequence variants
DNA was isolated from the surviving evolved populations at transfer 14 using the NucleoSpin Tissue kit from Macherey-Nagel and then subjected to short read sequencing with a HiSeq 3000 System (Illumina). DNA quality control for sequencing was done by the sequencing centre. Only DNA that passed the quality check was sequenced. Quality control of the Illumina reads was performed using MultiQC v1.12 [35] with default parameters. To perform strain-level abundance estimation in our metagenomic samples (the term is used here and below to indicate the genome sequencing data obtained for mixed-strain populations rather than single isolated clones), we generated variation plots, using StrainFLAIR v0.0.1 with default parameters [36]. Analysis of horizontal gene transfer was done using a subset of genomes sequenced with both short and long reads, generated with a HiSeq 3000 and a Sequel II (PacBio), respectively. Illumina short reads were assembled using SPAdes v3.15.5 in meta-assembly mode [37]. PacBio long reads were assembled using Flye v2.9.3-b1797 in metagenome assembly mode [38]. Short and long read assemblies were annotated using Bakta v1.9.2 with default parameters against the complete database [39]. We used the wild-type genomes from the gdPop collection and the short-read assemblies to build a pangenome using Panaroo v1.4.2 [40]. We ran the analysis in strict mode, with the flag to remove invalid genes and a protein family sequence identity threshold of 90%. To investigate the genetic rearrangements of the phage plasmid acquired from H13, we compared the relevant genomic regions in the wild-type genomes of H13 and H18 and the long-read assemblies of the metagenomic samples using clinker v0.0.20 [41] and a minimum alignment sequence identity of 90%. For samples with only short reads, the phage replicase from the phage plasmid was compared using DIAMOND (v 2.1.9) [42] with the blastx mode against all short-read assemblies.
The analysis of sequence variants was based on the short reads data. Reads were trimmed using Trim Galore (v 0.6.10). To assign mutations accurately, we leveraged the fact that each replicate at the end of the experiment contained a limited number of strains—specifically, one dominant and, in a few cases, one nondominant strain. Sequencing reads from each replicate were partitioned using BBsplit (v 38.18), which evaluated each read against all 12 reference genomes (the 12 strains, the evolution experiment was started with) simultaneously. Reads were assigned to the genome yielding the highest alignment score based on k-mer matches; only in instances where multiple genomes produced equally high scores were the reads assigned to several matching references. This approach ensured that only the highest-scoring reads were retained for the dominant or nondominant strains, effectively removing any reads with a superior mapping score to any of the other 11 ancestral genomes.
Following partitioning, reads assigned to the dominant or nondominant strains were re-aligned to their respective references using Bowtie2 (v 2.4.1) [43], and SNPs were called using FreeBayes (v v1.3.2-dirty). SNPs and INDELs were filtered for a quality score > 50 using BCFtools (v 1.14) [44] and GATK (v 4.0.5.1). SNPs/INDELs detected in the ancestral strain, as well as those found when mapping reads from the nondominant strain onto the reference strain, were removed using a combination of BCFtools and BEDTools (v v2.30.0). This pipeline allowed for the high-resolution assignment of mutations to specific strains within each sequenced replicate based on a population sequencing approach. To account for the high proportion of nondominant reads in samples NMDF1 and NMDF2 (both No Drug control), the quality threshold was increased to 150 upon visual inspection in these samples. The orthologues of known resistance genes (mexY, fusA, mexN, csrA, nalC, ftsI, mexR, ampC, ampR, mpl, dacC, nalD, dacB, ompQ, galU, amrR, mexX, parR, parS, PA14_45870, phoQ, ampDh3, mraY, ampD, pmrB, cysQ, amgS, ampDh2, orfH) involved in resistance against GEN or PTZ were identified using OrthoFinder [45]. SNPs and INDELs overlapping these genes were determined using BEDTools.
Evolution experiment 2
To test the effect of genetic diversity and population size (and thereby mutational supply) on adaptation to PTZ and GEN, we performed a second evolution experiment. The effect of diversity was tested by comparing the adaptation of the mixed population to the single strains H05 and H18. An increase in population size was achieved by conducting the experiment in 2 ml and 4 ml volumes while maintaining the same selective pressure (IC99.9). The IC99.9 was the concentration that inhibited the population by 99.9% relative to the no drug control of the same population. In total, three treatments were tested: no drug, GEN monotherapy, and PTZ monotherapy. Each of these was initiated with 20 independent biological replicates per treatment. Treatments were fully randomized on 48-deep-well plates (Sigma Aldrich AXYP9605). Each of the 12 strains was grown for 20 h in LB and then mixed in equal proportions (v/v). From independent mixtures, bacteria were inoculated into 48-deep-well plates containing either 2 ml or 4 ml of medium (40 μl inoculum for 2 ml and 80 μl inoculum for 4 ml). They were then grown for 24 h at 37°C and 750 rpm, following which the plates were sonicated to break any bacterial clumps. 2% (v/v) inoculum was then transferred to another plate containing fresh medium. Growth of the bacteria was followed by measuring endpoint OD_600_ and CFU/ml over time. The experiment lasted 5 days.
Fluctuation assay
We performed a classical fluctuation assay to measure the rates of resistance mutations towards GEN and PTZ [46, 47]. In short, one colony of the tested bacterium was inoculated into 10 ml of M9 and incubated for 20 h at 37°C and 150 rpm. Similar to the other experiments, the gdPop was prepared by mixing the overnight culture in equal ratios. Parallel cultures were initiated in 96-deep well plates, each well containing 1 ml with a 10^3^–10^4^ CFU/ml founding population. To control for strong bottlenecks affecting the gdPop rates, we used a larger founding population. Each strain and antibiotic combination included two independent experiments. The deep-well plates were incubated at 37°C and 300 rpm for 20 h. Thereafter, we plated a 1:10 dilution of the cells on M9 agar plates containing 4× MIC of either GEN or PTZ at a density of 10^5^–10^6^ cells/cm^2^. The mutants were counted after 48 h of growth at 37°C. Plates containing more than 250 mutants were counted as 250 mutants. Resistance rates were determined using the Lea-Coulson model with partial plating in webSalvador 0.1 [48].
Statistical analysis
Statistical analysis was carried out in R [49] using packages car, lawstat, lme4, vegan, MASS, and lmPerm. For evolution experiment 1, we used full factorial randomization analysis of variance (ANOVA) to evaluate the influence of “Meta-Treatment,” “Treatment,” and their interaction on either of the response variables “Time to adaptation,” “Robustness of adaptation,” “MIC on GEN,” and “MIC on PTZ” (aov(Response ~ Meta.Treatment * Treatment, data = datadf)). To assess statistical significance, we generated a null distribution of F ratios using 10 000 randomizations of the data. Post hoc testing was carried out using pairwise Wilcoxon rank sum tests, or for population survival using a logistic regression model: glm(Survival ~ Meta.treatment + Treatment, data = survival, family = binomial).
The effect of the meta-treatment on Shannon index was analysed using general linear models with the formula lmer(Shannon ~ meta-treatment + (1|season), data = diversity). To assess the influence of the meta-treatment and the treatment on the survival frequency of the gdPop strains, a PERMANOVA (adonis2 function of the R package vegan) with 10 000 permutations was conducted. To avoid modelling strain correspondence, the Bray–Curtis dissimilarity metric of the postexperimental survival (as a measure of beta diversity) was calculated instead of the original strain proportions, and these distances were used as the response in a full-factorial model. To assess whether pre-existing resistance and bacterial interaction significantly contribute to explaining the variance in strain survival, linear regressions were performed on the mixing and no mixing meta-treatments separately. A full factorial model was initially used, and stepwise model selection was conducted using the Akaike Information Criterion (AIC; stepAIC function from the R package MASS). Due to nonparametric and heteroscedastic residuals, the significance of the identified regression model was verified using a permutation test (lmp function from R’s lmPerm package).
Resistance rates were compared with the likelihood ratio test, using webSalvador 0.1, and a variety of initial values for the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\hat{m}_{c}}\end{document} , which denotes the combined estimate of the expected number of mutants from both strains under the null hypothesis of equal mutation rates (e.g. 0.1, 0.4, and 2.1), thus following analysis recommendations [50]. These variations all converged to the same likelihood ratio test statistic. All *P-*values were adjusted for multiple testing with the false discovery rate.
In the second evolution experiment, the effect of volume, antibiotic treatment, and strain on survival was assessed using logistic regression with the model glm(Survival ~ Volume + Antibiotic + Strain, data = datadf, family = binomial) whereas adaptation parameters were analysed with a randomized ANOVA using the model aov(Response ~ Volume * Strain, data = datadf) followed by post hoc analysis using pairwise Wilcoxon rank sum tests.
Results
Genetically diverse strains vary in their antibiotic resistance, pairwise interactions, and horizontal gene transfer potential
The gdPop consisted of 12 strains from the previously characterized mPact strain panel [31, 32], covering the known genomic diversity of P. aeruginosa (Fig. 1A). Before performing the evolution experiment with the gdPop, we experimentally characterized the strains’ antibiotic resistance profiles and their pairwise interactions with one another. We focused on two antibiotics, GEN and PTZ, which are used to treat P. aeruginosa infections and, more importantly, which are known to exert strong selective pressure on this pathogen, because—as we demonstrated previously—they interact synergistically with one another [51], they display reciprocal collateral sensitivity upon resistance evolution to each of the drugs [52], and they further express directional negative hysteresis (i.e. a short exposure to PTZ induces a transient physiological change that enhances the killing efficacy of subsequently applied GEN [53]). For these reasons, these two antibiotics provide well-characterized drug models that are ideally suited for the current experimental evolution analysis of AMR. We used the MIC based on antibiotic test strips as a quantitative measure for resistance (rather than a qualitative definition of resistance based on reaching clinically relevant breakpoints). We found that the strains differed in their resistances to the two antibiotics used, GEN and PTZ, with larger differences observed towards PTZ than GEN (Figs 1B, S1).
Strain characteristics of the P. aeruginosa gdPop. (A) The gdPop is composed of 12 strains from the major P. aeruginosa clone type (mPact) strain panel. These strains belong to the species’ two major phylogroups and are representative of its genomic diversity. (B) Resistance of the gdPop strains to gentamicin (GEN) and piperacillin/tazobactam (PTZ), shown as MIC. Bars show the mean, and error bars the standard error of the mean (SEM, n = 3). (C) Pairwise interaction matrix of the 12 gdPop strains. A negative interaction was defined as reduced biomass of recipient Strain A in donor Strain B’s conditioned medium, whereas increased biomass indicated a positive interaction. Relative biomass was calculated by dividing growth in nonself-conditioned medium by biomass in self-conditioned medium (n = 3). (D) The interaction matrix was used to calculate the effect of each focal strain’s supernatant (donor strain) on the biomass of other strains (recipient strains) and vice versa. The summation of the effect of focal strain supernatant on all other strains (one row in the matrix in C) gave the cumulative effect of focal strain supernatant (focal - > others). The summation of the effect of other supernatants on focal strain biomass (one column in the matrix) gave the cumulative effect on focal strain biomass (others - > focal). Bars show the mean, and error bars the SEM (n = 3). Data provided in Supplementary Data Tables 1 and 2.
We characterized bacteria–bacteria interactions by assessing the effect of any secreted compound of a donor strain (e.g. released metabolites) on the growth and survival of a recipient strain, thereby evaluating processes that act across space, independent of direct contact, and thus represent an important component of how bacteria interact with each other, following the rationale of similar previous analyses of bacteria–bacteria interactions [54–56]. These interaction effects between the strains were quantified in the absence of antibiotics as relative bacterial biomass in conditioned media (i.e. the supernatant of a culture of one of the bacteria). Bacterial biomass of a strain in its own conditioned medium (i.e. same–same conditions) was used as a baseline. The relative reduction of biomass of recipient Strain A in donor Strain B’s conditioned medium compared to the corresponding value of Strain A in its own conditioned medium was used as an indication for a negative interaction, whereas a relative increase in biomass was used as an indication for a positive one. It is possible that the same-same conditions may stimulate growth and could then introduce a bias when used as a baseline for the comparisons (e.g. a neutral interaction with another strain may then become negative). Even if this may happen, the relative effects remain unchanged and provide insights into the effects on biomass during an interaction. Taking this possible limitation into account, our analysis identified 69% of the 144 one-way interactions to reduce biomass of the recipient strain, 8% had no effect, and 22% increased recipient biomass. We further calculated two types of cumulative interaction effects, focusing either on the cumulative effect of a focal donor strain on all others (focal - > others) or vice versa (others - > focal). The cumulative effect of the focal donor strain’s supernatant (focal - > others) was calculated by summing the effect of its supernatant on all other strains. The cumulative effect on a focal recipient strain’s biomass (others - > focal) was calculated by summing the effect of other supernatants on focal strain biomass. These cumulative measures revealed substantial variation across the gdPop (Fig. 1D). As a last point, we note a high potential for HGT for the gdPop strains, demonstrated by our previous work that identified numerous integrative conjugative elements (ICEs), plasmids, as well as other mobile genetic elements for the strains (Table S1) [32].
Population mixing enhances adaptation under antibiotic selection
We evolved the gdPop under various antibiotic selection conditions to test the impact of standing genetic variation, bacterial interactions, HGT, and metapopulation structure on resistance evolution. Six antibiotic treatments were included—two monotherapies, two switching treatments, a combination treatment, and a no-drug control (Fig. 2A). GEN and PTZ were chosen because their sequential and combination treatments impose strong selective constraints on bacteria [51–53]. The selection strength of each treatment was standardized to reduce the cell number of the mixed population of 12 strains by three orders of magnitude compared to the no drug control (i.e. inhibitory concentration of 99.9% or IC99.9). This level of inhibition was chosen to exert a strong selection (near the MIC) on the whole population. Each treatment was replicated eight times. To assess the importance of metapopulation dynamics, we used two separate transfer protocols, resulting in two meta-treatments (Fig. 2B). The no-mixing meta-treatment involved transferring bacteria from one well of the growth plate to the same well in the new growth plate containing fresh medium, thereby simulating separated populations. In the mixing meta-treatment, all bacteria from replicate wells of a particular treatment were mixed in equal proportions at the end of each growth season, followed by transfer of this mixture to the fresh plates, thereby simulating populations connected through migration. The experiment included 15 transfers with the last day in drug-free media to assess population survival (Fig. 2A). The dynamics of evolutionary adaptation were quantified during the evolution experiment using two parameters: (i) Time to adaptation was the time taken until a population under antibiotics had reached the same biomass (measured as optical density, OD) as the simultaneously evolving no drug control, whereas (ii) the robustness of adaptation was calculated as the number of growth seasons the population then remained at biomass levels of the simultaneously evolving no-drug control after having already adapted. Time to adaptation indicates how fast populations are able to adapt to a specific selective pressure, as imposed by the antibiotic treatment. In treatments with antibiotic changes, however, there are changes in selective pressure across time. The parameter robustness of adaptation quantifies how robust the adaptation is to these changing conditions.
Experimental evolution of the gdPop under antibiotic selection. (A, B) To test the impact of standing genetic variation, bacterial interactions, HGT, and metapopulation structure, the gdPop was evolved under six different specific treatments (A), including monotherapies, switching treatments (i.e. alternation of GEN and PTZ, starting with either of the drugs in separate treatments), a combination treatment, and a no-drug control. All treatments were standardized to inhibit the gdPop by 99.9% and were initiated with eight independent biological replicates. An additional meta-treatment (B) was imposed on all six specific treatments to simulate metapopulation structure. The no-mixing meta-treatment involved a well-to-well correspondence during transfer after each growth season. The mixing meta-treatment involved mixing of all replicate populations within a particular specific treatment (shown in B) and inoculation of the next plate from this mixture. (C) Growth dynamics during the evolution experiment. Biomass (OD600) relative to the no-drug control is plotted. Thin lines represent replicate populations, the thicker line the mean, and the dotted line the no-drug control. n = 8 per treatment. (D) Adaptation response of the gdPop across specific treatments and meta-treatments. Adaptation was quantified using two parameters, the time to adaptation (i.e. the number of transfers until the growth level of the simultaneously evolving no-drug control was reached), and robustness of adaptation (i.e. the number of transfers the population reached the growth level of the no-drug control after having adapted once). Each data point is one replicate. Bar is the mean of all replicates per specific treatment. (E) Number of surviving populations per specific treatment out of a total of eight at the end of the evolution experiment. (F, G) Log2MIC fold change of evolved gdPop relative to the ancestral gdPop for PTZ and GEN. Two to eight populations were recovered from frozen material per treatment (No population was recovered for the PTZ - > GEN treatment in the mixing meta-treatment). Dots are the individual replicate populations, the bar is the mean, and the dotted line is the ancestral log2MIC. The adaptive responses in D, F, and G were analysed with a permutation ANOVA, followed by pairwise Wilcoxon rank sum post hoc tests. Statistically significant differences between specific treatments are indicated by different letters (see line denoted “stats”). Detailed results on the statistical analysis of the data are given in Tables S2-S10. All data for the figure panels are provided in Supplementary Data Tables 3–7.
The inferred evolutionary dynamics varied significantly between the two meta-treatments and the individual treatments (Fig. 2D and E, randomized ANOVA, P < .0001 for both parameters). Under no-mixing conditions, adaptation to GEN monotherapy happened in 2 days, whereas PTZ monotherapy exerted stronger selective constraints, with bacteria unable to adapt until the end (Fig. 2D). The switching treatments increased the time to adaptation compared to the GEN monotherapy but not to the PTZ monotherapy. The combination treatment also delayed adaptation. Populations that took longer to adapt also appeared to have a less robust adaptation. Adaptation was faster and more robust under the mixing conditions. In particular, the gdPop adapted significantly faster to PTZ monotherapy and the combination in the mixing rather than the no-mixing meta-treatment (Fig. 2D). The mixing meta-treatment was also associated with a significant increase in population survival (logistic regression, P = .001; Fig. 2F).
We characterized AMR of the evolved populations using MIC test strips (Fig. 2G and H) and broth microdilution (Fig. S1B). Relative resistance was defined as the increase in MIC of the evolved population over the ancestral gdPop. PTZ resistance under no-mixing conditions was similar and not higher than that of the ancestor (Fig. 2H). In the mixing meta-treatment, the no-drug treatment and GEN monotherapy remained at ancestral PTZ MIC level, whereas the others were significantly higher. Resistance to GEN was observed in all antibiotic treatments in both mixing and no-mixing meta-treatments except for PTZ monotherapy under mixing. Overall, mixing meta-treatment increased survival, adaptation, and AMR. Evolutionary adaptation was observed to be most constrained when PTZ was included.
Metapopulation structure, standing genetic variation, and bacterial interactions determine strain abundance under antibiotic selection
We assessed to what extent the treatments and meta-treatments (metapopulation structure) affected the strain composition of the evolving gdPop. We used strain-specific PCRs on three randomly chosen replicate populations from transfers 4, 7, 11, and 13 (Fig. S2) and calculated a Shannon index for each treatment combination. We identified temporal variation in the Shannon index within treatments (Fig. 3A) and a significantly higher Shannon index under mixing than no-mixing conditions (Fig. 3C, GLM, P < .001), observed until day 11 of the evolution (Fig. 3A).
Strain diversity and HGT in the gdPop. (A) Strain diversity during experimental evolution. Strain-specific PCRs (Table S11) were performed on three replicate populations from transfers 4, 7, 11, and 13 to obtain presence/absence data for the 12 strains, followed by calculation of Shannon index indices for each treatment and time point. Open circles are the individual replicates, and filled circles represent the mean. Error bars show standard deviation. (B) Whole genome sequences for populations from transfer 14 were used to determine strain abundance. Each bar corresponds to a replicate population. Strains are indicated by the different colours (see left) and also indicated by their strain names once in the graph. Empty columns refer to populations that went extinct (cf. Figure 2F) or could not be used for sequencing due to low DNA quality. HGT was only identified for a phage plasmid and indicated by green boxes right below the individual bars. (C) Results of the general linear mixed model used to assess the impact of meta-treatment on diversity. (D) Results of the PERMANOVA used to test the effect of experimental design on strain abundance at transfer 14. (E) Results of the permutation-based linear regression used to test the effect of pre-existing resistance and bacterial interaction parameters on strain abundance, assessed for the two meta-treatments separately. Only significant coefficients are shown. Detailed results of statistical analyses are given in Tables S12–S15, and the data in Supplementary Data Tables 8 and 9.
We assessed strain abundance and HGT using whole genome sequence data for populations from transfer 14 (Fig. 3B, Fig. S3). Only 4 out of initially 12 strains were still present at the end, and usually one strain dominated a particular population. Strain abundance significantly varied with the meta-treatment, treatment, and their interaction (Fig. 3D, PERMANOVA, P < .05 for all). We next tested if strain abundance was associated with starting strain MIC (i.e. standing genetic variation in MIC) and bacterial interaction parameters, as determined above prior to the evolution experiment (Fig. 1B–D). We found that the cumulative effect, which results from the conditioned medium of the focal strain on others (focal - > others), significantly decreased strain abundance (Fig. 3E, PERMLM, P < .0001) in both mixing and no-mixing meta-treatments. This is consistent with the observation that most strains negatively affected others. The MIC on PTZ increased strain abundance under no-mixing conditions (PERMLM, P < .0001). The interaction between PTZ MIC and the cumulative effect on focal strain growth (others - > focal) also significantly affected abundance in both meta-treatments (Fig. 3E, PERMLM, P < .0001). Finally, the three-way interaction between the PTZ MIC, focal - > others, and others - > focal was a significant predictor of strain abundance in the no-mixing meta-treatment (Fig. 3E, PERMLM, P < .0001).
If we consider the four predominant strains, then these observed patterns can be explained most easily for H18, which receives the strongest benefit from the interaction with other strains (Fig. 1D, bottom panel), and simultaneously outcompetes the others (Fig. 1D, top panel). Another strain, H08, mainly dominates in the GEN- > PTZ switching treatment under mixing conditions. This is possibly because H08 benefits from many of the pairwise interactions and is one of the stronger competitors (Fig. 1C and D). Strain H11 dominates in single replicates in treatments with PTZ under no-mixing conditions (Fig. 3B). This strain benefits from some of the interactions (Fig. 1C and D, bottom panel) and, importantly, shows one of the highest resistance levels towards PTZ (Fig. 1B, bottom panel). Finally, the strain H05 dominates the no-drug control and the GEN monotherapy treatments under mixing conditions and also a few replicate populations of different multidrug treatments under no-mixing conditions (Fig. 3B). This pattern cannot be easily explained with the available data, as this strain does not perform very well in the pairwise interactions and it also only shows somewhat higher resistance levels towards PTZ, but not GEN, when compared to the other strains (Fig. 1B–D). Overall, we conclude that the increased abundance of some albeit not all strains during antibiotic selection can be explained by the type of treatment, the pre-existing resistance to the antibiotics as well as the interactions between bacteria.
We identified HGT, in all cases of a phage plasmid from strain H13, which spread to new hosts in 59% of the populations (Figs 3B, S3). The genes located on this phage plasmid have not been reported to affect antibiotic resistance (Fig. S3). Therefore, the spread of the phage plasmid appears to be a consequence of selfish behaviour and seems unlikely to have contributed to AMR evolution.
Increases in antimicrobial resistance are associated with mutations in known antimicrobial resistance genes and variation in rates of de novo resistance emergence
We observed increases in GEN and PTZ resistance in populations dominated by the same P. aeruginosa strain. For example, the strain H05 dominated in the control treatment and the GEN monotherapy under mixing conditions, whereby we only observed an increased GEN resistance in the latter but not the former (Figs. 2G and 3B). Similarly, H18 was most abundant in PTZ monotherapy under mixing conditions, but also prevailed in several replicates of the control treatment under no-mixing conditions; an increased resistance was only recorded for the former but not the latter (Figs. 2H and 3B). These observations suggest that the resistance increase is not only due to changes in gdPop strain composition but additionally mediated by newly emerged resistance mutations. To assess this point, we characterized the distribution of mutations in known P. aeruginosa AMR genes across treatments and meta-treatments. We found that AMR gene mutations only occurred in treatments with antibiotics and, in this case, in almost all of the replicate populations, whereas no AMR gene mutation was identified for the no-drug controls (Fig. 4A). Moreover, the mixing meta-treatment led to the uniform spread of the same variants across replicate populations of a particular treatment, including variants in the gene ampR whenever a treatment contained PTZ and also variants in pmrB or parS for treatments with GEN (Fig. 4A). The no-mixing meta-treatment produced more variation in the favoured AMR gene mutations, including variation in the detected mutations in pmrB (Fig. 4A). Based on these results, we conclude that the mixing meta-treatment favoured a uniform selection of resistance variants; multiple beneficial mutants likely arose in the different populations but only the one with the highest competitive fitness dominated and distributed across the replicates. Under the no mixing conditions, similar mutants may have been selected for as under the mixing conditions. However, the separation of the populations prevented direct competition and individual populations fixed different mutant alleles.
The distribution of mutations in known resistance genes and variation in resistance rates towards GEN and PTZ. (A) Mutations in known AMR genes identified for populations from different treatments at the end of the evolution experiment. The known AMR genes are given along the vertical axis, whereas the different treatments and meta-treatments are shown horizontally. The size of the circles corresponds to the number of replicate populations with a variant in a particular gene. Within each circle, different background genotypes are indicated by colour, following the same colour scheme used in Fig. 3B (H05, yellow; H08, light red; H11, dark red; H18, purple). Different mutations within a specific genotype are shown as separate pieces. (B) Resistance rates towards GEN or PTZ. Resistance rates were inferred using the fluctuation assay with either the single strains H05 or H18 or the whole gdPop. The results are shown as the mean of two independent runs. Error bars represent confidence intervals (CI95). The statistical comparison between antibiotics was based on likelihood ratio tests. The four stars indicate a significantly lower resistance rate towards PTZ than GEN (P < .0001; detailed statistical results in Table S16 and the data for this figure in Supplementary Data Tables 10–12).
Although our findings can explain the observed resistance increases, they cannot explain the lack of resistance evolution against PTZ (Fig. 2G) or the increased extinction rates in treatments with PTZ (Fig. 2F). We postulated that these results may be due to a small mutational target size for PTZ resistance and thus a small rate of spontaneous PTZ resistance mutations. We tested this idea with standard fluctuation assays [46, 47] for H05, H18 (the two dominant strains at the end of evolution), and gdPop on both GEN and PTZ. We found that the resistance rates towards PTZ were indeed significantly lower than those towards GEN—consistently in both tested strains and the gdPop (Fig. 4B, Likelihood Ratio tests, P < .0001 for all). Therefore, the comparatively lower PTZ resistance rates appear to underlie the observed lower rates of resistance evolution and the higher extinction rates during experimental evolution in PTZ-containing treatments. Furthermore, the genome sequence analysis of populations from monotherapy treatments suggests that mutations in a larger number of genes appear to contribute to resistance to PTZ (up to four genes) rather than GEN (only one gene; Fig. 4A). Hence, the lower resistance rates may be due to either a smaller number of possible resistance mutations in these genes or a higher cost of the suitable mutations.
Strain diversity increases population survival under antibiotic selection in a second independent evolution experiment
We observed that metapopulation structure, standing genetic variation, and bacterial interactions affected evolution through their effect on population composition (Fig. 3D and E). Metapopulation structure, however, also affected the tempo of adaptation; well-mixed populations adapted much faster than the non-mixed populations (Fig. 2D). We postulated that the observed higher rate of adaptation was related to strain diversity as the mixing meta-treatment maintained a higher strain diversity across time (Fig. 3A and C). The higher strain diversity could result in a higher rate of adaptation due to (i) a reduced random loss of strains during the transfer bottlenecks, including those that already express resistance or show a high potential to evolve resistance, and/or (ii) the maintenance of beneficial bacterial interactions that somehow facilitate resistance due to maintenance of the contributing strains. The mixing conditions could, however, have also further enhanced adaptation by producing an overall larger population size (1 well in no-mixing vs 8 connected wells in mixing) and thus a higher likelihood for the occurrence of favourable new mutations. Any new mutation would then spread easily to the other connected populations. To explore the relevance of these alternative hypotheses, we performed a second independent evolution experiment that was specifically focused on testing the importance of strain diversity and population size. To test the effect of diversity, the evolution experiment was set up with either only single strains (H05 or H18) or the gdPop. The two considered single strains were most frequently observed at the end of the first evolution experiment, indicating their genetic potential to be maintained during experimental selection. If these were now unable to adapt, but the gdPop did adapt, it would suggest that (i) diversity played a role in adaptation in the first evolution experiment, and (ii) that the effect of diversity was not due to the retention of successful strains, but rather resulted from selective benefits arising from bacterial interactions. To explore a possible effect of population size, the experiment was carried out in two volumes, 2 ml and 4 ml, whereby we ensured that antibiotic selective pressure (i.e. antibiotic concentrations) and also cell density were identical across the two volumes (Fig. S4). If there were faster adaptation in the 4 ml treatment that would indicate the relevance of population size. The two single strains and the gdPop were taken in two different volumes and then subjected to three selective treatments across five days (Fig. 5A; 20 replicates per treatment combination). During the evolution experiment, we monitored bacterial growth by measuring biomass (OD_600_) (Fig. 5B and E) and colony-forming units (CFU/ml) (Fig. S5).
Validation evolution experiment with gdPop and two single strains in different volumes under antibiotic selection. The faster adaptation occurring in mixed-strain populations could result from either (i) higher diversity in the early part of the evolution and/or (ii) a larger combined population size (see text for more details). (A) Design of the repeat evolution experiment to explore a possible influence of population size. Population size was varied by evolving bacteria in two different volumes, 2 ml and 4 ml, allowing us to double the population size in the larger volume and simultaneously ensuring identical cell densities. Two single gdPop strains, H05 and H18 (the two most frequent survivors of the main evolution experiment, and therefore most likely to adapt on their own), and the whole gdPop were subjected to three treatments (two antibiotic monotherapies and a no-drug control) over five days (all with 20 replicates). Faster adaptation in the 4 ml volume would support the relevance of population size. Adaptation of the gdPop but not the single strains would support the relevance of diversity for adaptation, likely through its effect on maintenance of beneficial bacterial interactions. (B, E) Growth dynamics during the evolution experiments for either small (B) or large (E) volumes. Biomass (OD600) over time is plotted. Thin lines are the replicate populations, and the thicker lines are the means of the surviving (blue) and extinct (red) replicates (n = 17–20 per treatment combination). (C, F) Surviving populations per strain and treatment for small (C) or large (F) volumes. The total height of the bar represents the total replicates tested (with a possible maximum of 20). Black and white represent surviving and extinct replicate populations, respectively. (D, G) Adaptation parameters for populations growing in the GEN treatment for either small (D) or large (G) volumes. Adaptation was quantified as the time to adaptation and robustness of adaptation (see legend to Fig. 2). Each data point is one replicate. Bar is the mean of all replicates. Detailed statistical results are provided in Tables S17–S20, and data for the figure panels in Supplementary Data Tables 13–21.
We found that strain diversity but not strain identity or volume influenced the dynamics of adaptation of bacteria, contingent on the antibiotic used. Under GEN monotherapy, the proportion of extinct replicate populations varied between strains (Fig. 5C and F) and further depended on the initial CFU/ml that usually showed some variation under the used experimental conditions. In this case, we now observed that replicates with low initial CFU/ml were more likely to go extinct (Fig. S5). On PTZ, however, strain diversity positively impacted survival (logistic regression, P < .05): all single-strain replicates went extinct, whereas at least some gdPop replicates survived (Fig. 5C and F). Therefore, the presence of strains with the potential to spread during experimental selection (e.g. because of expressing antibiotic resistance) was not sufficient. By contrast, strain diversity was advantageous, likely due to selective benefits arising from bacterial interactions. We next assessed the pattern of adaptation across the two tested volumes, for the antibiotic treatment with little extinction, the GEN treatment, using the parameters of time to adaptation and the robustness of adaptation (Fig. 5D and G). The time to adaptation was not significantly different between volumes, but it showed a statistical trend of a difference among the three strains and the interaction of strain and volume (randomization ANOVA, P = .59 for volume, P = .08 for strain, and P = .056 for their interaction). The robustness of adaptation was significantly different between volumes and the strains (randomized ANOVA, P < .05 for both), and it was higher across strains in the 2 ml volume.
Our combined results strongly suggest that the higher survival and faster adaptation observed in the mixing meta-treatment was a consequence of higher strain diversity. The increased diversity enhanced adaptation of the bacterial populations not by limiting the random loss of strains with the genetic potential to survive, but most likely by retaining bacterial interactions that contribute to AMR.
Discussion
Many infections are polymicrobial and often include strain variation within species [11–13, 15]. In our proof of principle study, we here used a controlled experimental approach to test the effect of strain variation in resistance, strain–strain interactions, HGT, and metapopulation dynamics on AMR evolution in a genetically diverse pathogen population. We demonstrate that high strain diversity, as a consequence of the high migration treatment, increases population survival under antibiotic selection (Figs 5C, F, and 6), most likely as a consequence of strain interactions within the population. Indeed, we also found that the selectively favoured strains at the end of evolution are those that benefit the most from interacting with other bacteria and have high pre-existing resistance (Figs 3B and 6). In the following, we will discuss the importance of metapopulation dynamics, bacterial interactions, and standing genetic variation on antibiotic resistance evolution, and also address the unusual observation of low resistance evolution towards PTZ.
Strain variation in resistance, bacterial interactions, and metapopulation dynamics determine selection outcome. Summary of the main findings of our study. We assessed to what extent strain variation in antibiotic resistance (determined as the MIC), bacterial interactions, horizontal gene transfer (HGT), and metapopulation dynamics shape adaptation of a gdPop of P. aeruginosa to antibiotic selection. We found that pairwise bacterial interactions, MIC, and metapopulation dynamics are key predictors of strain composition. Bacterial interactions maintained through a lack of spatial separation (mixing) were further found to contribute to population survival through their beneficial effect on the focal strain H18. Created in BioRender [57].
Metapopulation dynamics—here understood as the subdivision of a species into local populations connected by migration—occur in various infections and depend on the exact infection characteristics. In CF lungs, the spatially structured environment allows migration of the infecting pathogen like P. aeruginosa, although usually at low frequency [58]. By contrast, lung-infecting Mycobacterium tuberculosis is frequently disseminated between lung regions and to extrapulmonary sites [59]. Nonpulmonary infections, such as those by Helicobacter pylori, migrates between distinct compartments of the stomach [60]. Such dynamics play a key role in the adaptation of microbial communities, affecting biodiversity and competition [61]. The influence of such migration-mediated structure on adaptation has been addressed comprehensively using theoretical approaches [62–64]. When migration is limited, populations are separated, and beneficial mutations spread independently within subpopulations, potentially extending the duration of selective sweeps [65] and facilitating the coexistence of multiple genotypes [66, 67]. Spatially separated populations connected by low levels of migration have indeed been shown to maintain high genetic diversity under migration and serial dilution [29], and, when beneficial mutations are rare, migration can also boost their spread [68]. In our experiments, the no-mixing and mixing meta-treatments represent two extreme cases of migration: in the no-mixing meta-treatment, replicate wells evolved as separated populations without migration between wells, whereas in the mixing meta-treatment all replicate wells of a given treatment were repeatedly mixed, thereby simulating populations connected through migration (Fig. 2B). Under these conditions, populations evolving without migration showed a decrease in the speed of adaptation (Fig. 2D and E) and also produced a larger diversity of sequence variants between populations (Fig. 4A left). At the same time, populations connected through migration showed a more uniform spread of sequence variants (Fig. 4A right) and an at least initial increase in strain diversity (Fig. 3A), the latter possibly because of a reduced likelihood of random strain loss due to bottlenecks at transfer. Taken together, these findings offer a proof of concept for the general importance of meta-population dynamics in AMR evolution and highlight that such structuring can play a critical role in shaping AMR dissemination dynamics. Higher adaptation rates are associated with higher strain diversity under mixing conditions (Figs. 2D, E, and 3A, and especially Fig. 5B–G), likely due to maintaining beneficial interactions between the bacteria. Bacterial interactions were previously implicated in modulating antibiotic resistance [25, 26]. For example, cross-feeding can decelerate de novo resistance evolution when one partner always has to “wait” for its cross-feeding partner to evolve higher resistance [19, 69]. As another example, the presence of carbapenemase-producing Stenotrophomonas maltophilia allowed coexisting P. aeruginosa to survive and subsequently evolve imipenem resistance [21]. Studies assessing the impact of a more complex community on resistance evolution are fewer and usually follow the dynamics of a focal strain within a community subjected to antibiotics. Interspecific competition increased the costs of resistance in focal Escherichia coli resulting in higher minimal selective concentrations [20]. However, a subsequent increase in antibiotic concentration led to the competitive release of the focal strain. Another study found that resistance evolution in focal E. coli was suppressed in human faecal communities [70]. In our study, the ultimately dominant strain H18 (Fig. 3B) benefitted from the others the most and simultaneously outcompeted them (Fig. 1C and D). Thus, the overall distribution of interactions impacts how the whole population adapts to novel selective pressure, like that imposed by antimicrobials.
Strain abundance in our gdPop was significantly influenced by pairwise bacterial interactions as well as standing genetic variation in AMR (Fig. 3B and E). Our finding of the importance of microbial interactions is in line with previous work that showed that competition could predict survival in multispecies [71] or multistrain communities [72]. Moreover, standing genetic variation can also contribute to adaptation by providing the necessary diversity to adapt [73]. For example, in patients colonized by multiple P. aeruginosa strains, the pre-existence of resistant strains facilitated the spread of AMR [23]. Altogether, strain interactions as well as standing genetic diversity are important predictors of population dynamics.
The inability of bacteria to adapt to PTZ as compared to GEN was likely caused by a low rate of PTZ resistance mutations, as indicated for both individual strains and the mixed-strain gdPop by the results of the fluctuation assay (Fig. 4B). Since we standardized concentrations of both antibiotics to inhibit the growth of individual strains or populations to the same extent during our experiments, the observed difference in adaptation rates are unlikely due to variable expression of antibiotic-degrading enzymes during growth and, thus, are best explained by differences in the rate of emergence of resistance mutations. As a consequence, the inclusion of PTZ in switching and combination treatments increased their efficacy (Figs. 2D–F and 5B–G). These findings are consistent with our previous demonstration that the inclusion of an antibiotic with low resistance rates in switching treatments with only beta-lactams increased bacterial eradication [74]. In the current study, we confirmed this effect, using different antibiotics than before, and observed in both switching and combination treatments.
Our focus on a highly diverse 12-strain community is generally consistent with multistrain infections in patients. Such infections span a wide range of strain diversities, from apparent single-strain infections to heterogeneous populations with multiple clone types [14, 75]. Even if infections are initiated by a single clone, they can rapidly diversify in vivo through mutation, spatial segregation, and adaptation [15, 76], generating multiple coexisting lineages with distinct resistance profiles and interaction phenotypes [58]. In consideration of these previous observations, we expect the general principles identified here—namely the influence of strain-specific resistance, meta-population structure, and bacteria–bacteria interactions on infection characteristics and resistance evolution—to apply to in vivo conditions. These principles should also shape communities with fewer strains, although the resulting outcome may differ. In less diverse communities, strain composition and resistance should depend more strongly on the presence or absence of highly resistant strains, the specific structure of interaction networks, and also the antibiotics used (see divergent results for GEN and PTZ). As a consequence, we anticipate that the evolutionary dynamics in low-diversity communities will be more stochastic, with reduced buffering against bottlenecks and antibiotic stress. In contrast, higher strain diversity increases the likelihood of the presence of resistant strains and beneficial interactions, thereby stabilizing population-level responses under strong selection.
Predicting how populations or communities change over time is of particular interest in microbial ecology and additionally of high relevance in medical infectiology. Our work demonstrates that the presence of strain diversity enhances the ability of bacteria to counter the high selection imposed by antibiotics (Fig. 6). Further, we show that inter-strain interactions are as important as inter-species interactions in their effect on evolution. The consideration of these ecological processes is key to understanding antibiotic resistance spread in the widely occurring polymicrobial infections. The role of intraspecific variation during bacterial AMR evolution clearly warrants further investigation covering aspects that were not addressed in our current work. Future work could explore longer time periods for experimental evolution, additional clinically relevant antibiotics, and targeted manipulation of strain interaction networks, for example, by including or excluding central strains. More explicit tests on the relevance of HGT, using strains with additional types of mobilizable genetic elements, as well as an assessment of AMR evolution within a host, would further strengthen links to the clinical context and help clarify how host-associated environments shape bacterial interactions and resistance evolution.
Supplementary Material
Supplementary_data_file_wrag039
Supplement_Batra_etal_revised_Clean_wrag039
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Merker M, Tueffers L, Vallier M. et al. Evolutionary approaches to combat antibiotic resistance: opportunities and challenges for precision medicine. Front Immunol 2020;11:1938. 10.3389/fimmu.2020.0193832983122 PMC 7481325 · doi ↗ · pubmed ↗
- 2Davies J, Davies D. Origins and evolution of antibiotic resistance. Microbiol Mol Biol Rev 2010;74:417–33. 10.1128/MMBR.00016-1020805405 PMC 2937522 · doi ↗ · pubmed ↗
- 3Murray CJ, Ikuta KS, Sharara F. et al. Global burden of bacterial antimicrobial resistance in 2019: a systematic analysis. Lancet 2022;399:629–55. 10.1016/S 0140-6736(21)02724-035065702 PMC 8841637 · doi ↗ · pubmed ↗
- 4Naghavi M, Vollset SE, Ikuta KS. et al. Global burden of bacterial antimicrobial resistance 1990–2021: a systematic analysis with forecasts to 2050. Lancet 2024;404:1199–226. 10.1016/S 0140-6736(24)01867-139299261 PMC 11718157 · doi ↗ · pubmed ↗
- 5Andersson DI, Balaban NQ, Baquero F. et al. Antibiotic resistance: turning evolutionary principles into clinical reality. FEMS Microbiol Rev 2020;44:171–88. 10.1093/femsre/fuaa 00131981358 · doi ↗ · pubmed ↗
- 6Hjort K, Jurén P, Toro JC. et al. Dynamics of extensive drug resistance evolution of Mycobacterium tuberculosis in a single patient during 9 years of disease and treatment. J Infect Dis 2022;225:1011–20. 10.1093/infdis/jiaa 62533045067 PMC 8921999 · doi ↗ · pubmed ↗
- 7Tueffers L, Barbosa C, Bobis I. et al. Pseudomonas aeruginosa populations in the cystic fibrosis lung lose susceptibility to newly applied β-lactams within 3 days. J Antimicrob Chemother 2019;74:2916–25. 10.1093/jac/dkz 29731355848 · doi ↗ · pubmed ↗
- 8Woods RJ, Barbosa C, Koepping L. et al. The evolution of antibiotic resistance in an incurable and ultimately fatal infection: a retrospective case study. Evol Med Public Health 2023;11:163–73. 10.1093/emph/eoad 01237325804 PMC 10266578 · doi ↗ · pubmed ↗