Interspecies interaction alters the trajectory of antibiotic resistance evolution by amplifying negative fitness epistasis
Suraya Muzafar, Ramith R Nair, Dan I Andersson, Omar M Warsi
TL;DR
Interspecies competition in microbial communities can shape how antibiotic resistance evolves in Salmonella by altering genetic pathways.
Contribution
The study shows interspecies competition amplifies negative fitness epistasis, leading to higher antibiotic resistance in Salmonella.
Findings
Interspecies competition selects S. enterica mutants with higher streptomycin resistance.
Competition increases mutations following negative fitness epistasis trajectories.
The aadA gene's expression is enhanced due to interspecies interaction.
Abstract
Interspecies interactions can influence the physiology of competing species, shaping their long-term evolutionary trajectories. Although interspecific competition’s role in community dynamics is well-documented, its impact on evolutionary outcomes and mechanisms is less explored. Here, we investigate how interspecies competition affects antibiotic resistance evolution in the gut pathogen Salmonella enterica within synthetic microbial communities. Specifically, we examine how the presence of an interspecific competitor, Escherichia coli, modulates resistance evolution at low streptomycin concentrations. Our findings reveal that interspecies competition results in the selection of S. enterica mutants with higher resistance levels by increasing the likelihood of accumulating resistance mutations that follow a trajectory of negative fitness epistasis. We show that this effect is driven by…
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- —Swedish Research Council10.13039/501100004359
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 · Evolutionary Game Theory and Cooperation · Mathematical and Theoretical Epidemiology and Ecology Models
Introduction
Interspecies interactions are fundamental to ecological and evolutionary dynamics, shaping species’ coexistence, adaptation, and community structures [1–4]. Although several studies have investigated the ecological consequences of these interactions, their evolutionary implications, particularly in shaping responses to abiotic environmental conditions, are emerging as a focus of interest [5, 6]. For example, studies have shown that interspecies interactions could interfere with adaptation to abiotic conditions due to fitness trade-offs, as has been demonstrated for adaptation to elevated CO_2_ [7] and cold temperatures [8]. Similarly, the presence of competitors has been shown to alter a species’ evolutionary trajectory by affecting the availability of resources [9, 10]. Furthermore, interspecific interactions can also affect the effective population sizes of different species, impacting the emergence and fixation rates of beneficial mutations [5, 11–13]. Despite these studies highlighting the role of interspecific interactions in shaping evolutionary responses, how competition-induced physiological changes influence these responses, particularly within a fitness-phenotype landscape, remains poorly understood.
Understanding how physiological changes induced by interspecies interactions alter the evolutionary responses of a given species is crucial as studies have shown that these interactions have a profound impact on physiology. These effects manifest as induction of stress responses, alteration of gene expression patterns, or metabolic shifts [14–16]. For example, prey species of snowshoe hares (Lepus americanus) exhibited increased stress hormone levels in the presence of avian and mammalian predators, resulting in a decline in reproduction rates [17]. Similarly, perennial plant species of Achillea millefolium and Veronica chamaedrys showed altered physiologies in response to ozone exposure when grown with interspecific competitor Poa pratensis compared to when grown alone [18]. These examples underscore how interspecies interactions can reshape physiological processes across a wide range of taxa. These physiological changes and their effects on evolutionary trajectories are intricately linked to the properties of the corresponding fitness-phenotype landscapes. These landscapes serve as a multidimensional representation of how different phenotypes interact with one another to affect fitness [19, 20]. These interactions can manifest in distinct ways, either through additive effects, where phenotypes independently contribute to fitness [21], or through epistatic effects, where interactions affect fitness non-linearly [22–24]. They can thus result in either facilitating or constraining the paths available to organisms for achieving higher fitness by altering the fitness effects of subsequent mutations, thereby reshaping the set of mutational paths available to evolving populations [25, 26]. Given these studies, it is crucial to explore the properties of these fitness landscapes to understand the mechanisms that govern how interspecies interactions affect a species’ evolutionary response.
Synthetic microbial communities, owing to ease of manipulation and the ability to use -omics methodologies, provide an ideal model for studying the mechanisms underlying the evolutionary consequences of interspecific interactions. Unlike many macroorganisms, microbes can be manipulated at the genetic, physiological, and population levels, allowing detailed mapping of fitness landscapes under controlled competitive conditions [27–29]. To this end, we use a synthetic microbial community comprising Salmonella enterica serovar Typhimurium LT2 (hereafter referred to as S. enterica) and Escherichia coli (hereafter referred to as E. coli) to investigate how interspecific competition affects antibiotic resistance evolution under low concentrations of the aminoglycoside antibiotic streptomycin (STR). S. enterica and E. coli have the potential to interact in the mammalian gut, and under certain microbiota or environmental contexts, E. coli (and other commensals) can suppress S. enterica colonization by competing for limiting nutrients such as sugars or iron, or by altering gut conditions [30, 31]. This pair, therefore, provides a biologically relevant model to investigate how interspecific competition influences the evolutionary dynamics of antibiotic resistance. Several studies have demonstrated the selection of de novo arising high-level resistance mutations at very low antibiotic concentrations [32–34]; however, how this process is influenced by interspecies interactions remains an open question.
Our findings reveal that interspecific interactions significantly alter the fitness landscape experienced by S. enterica when exposed to low STR levels. In particular, competition with E. coli increases the likelihood of resistance mutations exhibiting negative fitness epistasis (the fitness effect of combined mutations is lower than expected from additivity), resulting in larger phenotypic shifts in resistance levels. We further demonstrate that these effects are mediated by stringent starvation response pathways activated under competitive conditions, which enhance the expression of a cryptic aminoglycoside resistance gene (aadA). By linking ecological competition-induced physiological changes to genetic adaptation, our study provides new insights into the evolutionary dynamics of antibiotic resistance.
Material and methods
Growth conditions and bacterial cultures
Two bacterial species were used in this study that included the common human gut commensal E. coli K-12 MG1655, and the gut pathogen S. enterica serovar Typhimurium LT2. The E. coli used contained a chromosomal copy of a gene encoding a yellow-fluorescent protein, whereas the S. enterica contained a chromosomal copy of a gene encoding a blue-fluorescent protein. The nutrient-rich medium used in this study was Lysogeny broth (LB; 1 g Yeast extract, 1.5 g Tryptone, and 1.5 g Sodium chloride in 150 ml deionized water). All strains were grown at 37°C.
Evolution experiments and determining the frequency of resistant mutants
Evolution experiments were performed as described in our previous work [35]. In this earlier work, we also showed that in the nutrient-rich media used here, E. coli and S. enterica interact only through resource competition and coexist via negative frequency dependence. Briefly, 24 independent lineages were evolved in parallel, with both strains grown either as monocultures or in co-cultures. These experiments were conducted in microtiter plates, with each population growing in a volume of 200 μL of LB containing 1 μg/ml STR at 37°C (below the minimum inhibitory concentration (MIC) for both species, Supplementary Fig. 1). After 24 h of growth, 1 μL was transferred into 200 μL of fresh STR containing LB medium, resulting in ~7.5–8 generations per transfer. After every 8–10 transfers, plates were stored at −80°C by adding 20 μL of DMSO to each well, mixing thoroughly, and sealing the plates with adhesive foils (VWR 60941–076). Two lineages in the co-culture conditions were found to be contaminated and were thus excluded from further analysis (Supplementary Fig. 2).
To determine the emergence and frequency of resistant mutants, 1 μL from each lineage was plated on STR plates containing 16 μg/ml or 64 μg/ml STR. Because the focus was on the general emergence of resistant mutants, four lineages were plated together on a single plate. These frequencies were also compared with control experiments in which S. enterica monocultures and S. enterica–E. coli co-cultures were evolved for 640 generations in the absence of streptomycin. At the end of the experiment, each lineage was also plated separately at 64 μg/ml STR to isolate resistant mutants individually, which were then further analyzed for levels of resistance and fitness measurements. In a small number of cases (three each in mono- and co-cultures), where 1 μL did not yield resistant mutants at 64 μg/ml, 5 μL was plated. Complete lawn growth on selective plates was considered indicative of resistance mutations reaching fixation. Throughout the evolution experiment, frequencies of E. coli and S. enterica were measured using flow cytometry (MACSQuant VYB, Miltenyi Biotec; details provided below), and the frequencies of S. enterica and E. coli from the appropriate time points were then used to determine the frequency of resistant mutants (Supplementary Tables 1, 2, 3, and Supplementary Fig. 3).
Measuring the minimum inhibitory concentrations, exponential growth rates, and stationary phase densities
Resistance levels in the mutants were determined by measuring the MICs for STR using Etests (BioMérieux), which are quantitative antimicrobial susceptibility tests that use a plastic strip impregnated with a gradient of antibiotic to determine the MIC of a microorganism (Fig. 1A). Two biological replicates were used for each mutant (Supplementary Tables 4 and 7). For S. enterica resistant mutants isolated from lineages evolving under co-culture conditions, MICs were measured in the absence of cognate E. coli (i.e. E. coli from the same lineage), because at streptomycin concentrations corresponding to the MIC of S. enterica resistant mutants, cognate E. coli does not grow and therefore cannot influence the measurement.
Effect of interspecies interaction on resistance levels and growth characteristics of S. enterica mutants selected at sub-MIC of STR. (A) Resistance levels of mutants were quantified using Etests. Two biological replicates were used in each case. The Y-axis represents the minimum inhibitory concentration (MIC) for STR, and the X-axis represents the absence (monoculture) or presence (co-culture) of interspecific competition. A Wilcoxon rank-sum test was performed to determine statistically significant differences. (B) The fitness costs of resistance mutations were determined by measuring relative fitness, relative exponential growth rate, and relative stationary phase density (both normalized to the values for the ancestral S. enterica) for 24 and 22 resistant clones isolated from monoculture and co-culture conditions, respectively. Relative fitness values represent the average of four biological replicates per clone, with the fitness of the co-culture isolated mutant measured in the presence of the cognate E. coli (isolated from the same lineage). Exponential growth rate and stationary phase measurements are based on the average of three biological replicates. The line in the center represents the median value. The Y-axis represents selection coefficients (where above 0 implies that the resistant mutant wins, and below 0 implies that the wild-type ancestral strain wins), relative exponential growth rate, or relative stationary phase density, respectively. The X-axis represents the absence (monoculture) or presence (co-culture) of interspecific competition. A permutation-based ANOVA was performed to determine statistically significant differences.
To determine the cost of resistance mutations, exponential growth rate and stationary phase density were determined by generating growth curves for the ancestral and mutant S. enterica strains using a BioscreenC analyzer. The experiment was initiated by making a 1:1000 dilution of an overnight-grown culture in fresh LB medium without STR. The cultures were grown at 37°C with shaking and optical density measurements were taken at 600 nm (OD_600_) every 4 min. The exponential growth rate was measured by fitting the OD_600_ values between 0.02 and 0.09 to an exponential growth equation N = N*e^rt^ using the KaleidaGraph software (v4.5.3), where r (min^−1^) represents the exponential growth rate. The stationary phase density was calculated using the R-package growthcurver [36]. Three biological replicates were used in each case. Costs were determined as the reduction in growth rates and stationary phase densities of the mutant strains relative to the values of ancestral S. enterica and are presented as such (Fig. 1B). To measure growth parameters for the ancestral strains in a co-culture setup, we grew them together and measured their population sizes using flow cytometry, as both these strains are marked with different fluorescent markers. Measurements were taken after 0, 1, 2, 2.75, 3.5, 4.25, 5, 5.75, 6.5, 7.25 and 24 h. The exponential growth rate was calculated as the maximum slope across the time-course data. The frequencies at 24-h time point were used as stationary phase density values in each case (Supplementary Fig. 4).
Measuring gene expression for genes iraP and aadA using qPCR
To determine if interspecies competition with E. coli results in upregulation of the stringent starvation response and of the aminoglycoside transferase gene, qPCR analyses were performed for the genes iraP and aadA. The gene iraP is induced under nutrient-limiting conditions, is ppGpp-dependent, and serves as an indicator of an activated stringent response [37]. Total RNA was isolated under three conditions. First, S. enterica cultures were grown either alone or co-cultured with E. coli in LB at 37°C until late stationary phase (OD600 ~ 2, Fig. 2A). Second, cultures were grown alone or in co-culture until the exponential phase (~3.5 h post-inoculation of 1 μL of overnight cultures of each species into 200 μL of medium; Fig. 2A and Supplementary Fig. 4). Third, S. enterica was grown in either its own spent media or spent media from overnight-grown E. coli, with the spent media reconstituted with fresh nutrients. Bacterial cells were obtained by centrifugation, resuspended in RNAprotect reagent, incubated at room temperature for 10 min, and then pelleted again. RNA was then isolated using Qiagen’s RNeasy kit, and the cDNA was then generated using High Capacity Reverse Transcription Kit (Applied Biosystems), following previously described protocols [38]. To ensure that the qPCR primers for aadA, iraP, and the housekeeping gene hcaT were S. enterica-specific and would not amplify E. coli sequences under co-culture conditions, colony PCR was conducted on the E. coli strain used in this study. No products were observed for aadA, iraP and hcaT. Two biological and two technical replicates were used in each case. The average of the two technical replicates was used as the final value for each biological replicate. The error bars represent standard deviations.
Interspecific competition results in the upregulation of the stringent starvation response and the aadA gene. (A) Relative expression levels of genes iraP and aadA, normalized to the expression of the housekeeping gene hcaT. Expression levels in monoculture are used as the baseline, with their value set as 1, and all values are plotted relative to the average expression level observed under monoculture conditions. Two biological replicates were used in each case, with the value of each biological replicate obtained by averaging two technical replicates. Error bars represent the standard deviation. A Student’s t-test was used to determine statistically significant differences. (B) Outcomes of competition between an ancestral S. enterica and the ΔaadA S. enterica mutant (plotted as S. enterica/ΔaadA S. enterica), and (C) between these strains and E. coli (plotted as S. enterica/E. coli and ΔaadA S. enterica/E. coli), in the presence and absence of STR. The Y-axis represents the selection coefficient. Four replicates were used in each case, and error bars represent the standard deviation. Positive selection coefficients indicate that the strain listed first outcompetes the second, whereas negative values indicate that the second strain outcompetes the first. A permutation-based ANOVA was performed to determine statistically significant differences.
Measuring the fitness of resistant mutants by performing competition experiments
Fitness measurements were performed by competing appropriate fluorescently-tagged strains and measuring the rate of change in frequencies using flow cytometry (MACSQuant VYB, Miltenyi Biotec), as has been described elsewhere [35]. Briefly, for each competition experiment, 1 μL of overnight-grown cultures of each strain was mixed in 200 μL of fresh media. The competition experiment was conducted over 2 days (~15 generations), with 1 μL of overnight-grown competition mixture transferred to 200 μL of fresh media daily. The frequencies of different strains were measured by inoculating 3 μL of cultures to 200 μL of phosphate-buffered saline and measuring the frequencies of the fluorescently tagged strains using flow cytometry. Four biological replicates were used unless stated otherwise. The selection coefficient (s) was calculated by plotting natural logarithms of ratios of population sizes against time [35]. All competitions involving co-culture isolated mutants were done in the presence of the cognate evolved E. coli, isolated from the same population. For competition experiments between E. coli and S. enterica, both species were able to coexist in LB. To assess changes in fitness of the evolved clones relative to each other, we ensured similar starting densities across all competitions and measured the rate of change in their population densities over ~15 generations, before these reached a steady state (Fig. 2C). This approach is similar to what has been described elsewhere [35, 39].
Determining epistatic interactions between the gene aadA and the enriched resistance mutations
To determine epistatic interactions between the aminoglycoside transferase gene aadA and the resistance mutations accumulated in our experiments, we deleted the aadA gene from a total of 20 resistant mutants (10 each from monoculture and co-culture conditions) isolated in our study. The gene deletion was performed by replacing the aadA gene with a mini-Tn10dTet carrying a tetracycline resistance marker, using P22 transduction. The P22 transduction protocol was followed as has been described elsewhere [38]. The lysate for the transduction was generated on a previously constructed S. enterica where the aadA gene was replaced by the mini-Tn10dTet [40], and the transductants were selected on plates containing 10 μg/ml tetracycline. The location of the aadA gene on the chromosome was verified in relation to other mutated genes to ensure that the chance of cotransduction of wild-type alleles at other locations was negligible. The aadA gene was also deleted from the ancestral S. enterica strain, which was subsequently used in competition experiments against all the other constructed strains. The fitness of the enriched resistant mutants, in the presence of aadA, was calculated by competing the resistant mutant against the ΔaadA S. enterica and is termed s_res + aadA_. The fitness advantage conferred by the aadA gene was determined by competing the ΔaadA S. enterica with the ancestral S. enterica and is termed s_aadA_, whereas the fitness advantage conferred by the enriched resistance mutations selected in our study was determined by competing the ΔaadA S. enterica with ΔaadA-resistant mutant and is termed s_res_. The sum, s_aadA_ + s_res_, gives the predicted value of fitness of the enriched resistant mutants, whereas the standard deviation for this predicted value was calculated using the standard deviation measurements of s_aadA_ and s_res_. Positive and negative fitness epistasis are defined when the predicted selection coefficient is either lower or higher than the observed selection coefficient, respectively (Fig. 3).
Fitness-epistasis effects mediated by resistance mutations and the aadA gene. (A) Outline of competition experiments performed between resistant mutants and the ΔaadA S. enterica mutant, ΔaadA-resistant mutants and the ΔaadA S. enterica mutant, and S. enterica and the ΔaadA S. enterica mutant. In each case, the selection coefficient is calculated for a given strain with respect to the ΔaadA S. enterica mutant. (B) These competition experiments were used to determine the prevalence of positive and negative fitness epistasis for resistant mutants isolated from monoculture and co-culture conditions. The Y-axis represents the difference between observed (Sres + aadA) and predicted (Sres + SaadA) selection coefficients. Values above 0 represent positive fitness epistasis, and values below 0 represent negative fitness epistasis. A Wilcoxon rank-sum test was performed to determine statistically significant differences. Comparison of selection coefficients determined for (C) resistant mutants and (D) ΔaadA-resistant mutants isolated from co-culture and monoculture conditions, when competed against the ΔaadA S. enterica mutant. A permutation-based ANOVA was performed to determine statistically significant differences for the data presented in C and D.
Whole genome sequencing of resistant mutants
To identify the resistance mutations enriched in our experiments, 18 S. enterica STR-resistant clones were whole-genome sequenced. 1 ml of overnight culture was used for DNA extraction using the MasterPure Complete DNA & RNA Purification Kit (Epicentre) according to the manufacturer’s instructions. Illumina’s Nextera XT kit was used to make libraries (2 × 300) that were sequenced on a MiSeq System (Illumina). Samples were dual-indexed and pooled together. The average whole-genome coverage per sample was ~30X. Single-nucleotide polymorphisms (SNPs) and short insertions/deletions were determined by mapping the fastq files obtained from MiSeq sequencing to the reference genome of the ancestral S. enterica using CLC genomics workbench version 8. To identify structural rearrangements, insertions, and deletions, Breseq (version 0.38.1) was used [41], with these results further being confirmed using the CLC genomics workbench. The fastq files used for whole genome analysis for the resistant clones can be found on NCBI’s sequence read archive (SRA, BioProjectID: PRJNA1219796).
Statistical testing
To assess statistical significance, Student’s t-tests, Wilcoxon rank-sum tests, and permutation-based ANOVA were applied as appropriate. Data were first evaluated for normality using the Shapiro–Wilk test and for homogeneity of variances using the F-test; when these assumptions were violated, non-parametric tests were used. Specifically, Student’s t-tests were used to evaluate differences in gene expression levels of iraP and aadA in the presence and absence of interspecies competition (Fig. 2A). Wilcoxon rank-sum tests were used to assess differences in MICs between mutants isolated from monoculture and co-culture conditions (Fig. 1A), and for determining whether the prevalence of fitness epistatic patterns differed between co-culture and monoculture isolated mutants.
A permutation-based ANOVA was performed to determine statistically significant differences between the costs to relative fitness, exponential growth rates and stationary phase density observed for the resistant mutants isolated from monoculture and co-culture conditions (Fig. 1B), between selection coefficients for competition experiments involving the ancestral S. enterica and ΔaadA mutant (Fig. 2B), ancestral S. enterica and E. coli (Fig. 2C), ΔaadA S. enterica mutant and E. coli (Fig. 2C), between the fitness values for resistant mutants isolated from monoculture and co-culture conditions (Fig. 3C), and between the fitness values of these resistant mutants with the gene aadA deleted (Fig. 3D). To identify predictors of antibiotic resistance, MIC values were log₂-transformed prior to analysis to account for their two-fold dilution structure and to improve model assumptions. A multiple linear regression was fitted with log₂(MIC) as the response variable and stationary-phase density, exponential growth rate, S. enterica population density at the first and final cycle, and the frequency of high-level resistant mutants from the respective populations as predictor variables. Model assumptions were assessed, with the residuals not deviating significantly from normality (Shapiro–Wilk test: W = 0.973, P = .36), and no evidence of heteroscedasticity was detected (Studentized Breusch–Pagan test: BP = 5.73, df = 5, P = .33). To assess predictors of epistatic deviation, an additional multiple linear regression was performed using the same explanatory variables, together with log₂(MIC) of each strain. Model assumptions were assessed, with the residuals not deviating significantly from normality (Shapiro–Wilk test: W = 0.97, P = .77), and no evidence of heteroscedasticity was detected (Studentized Breusch–Pagan test: BP = 6.8 df = 6, P = .33). All statistical analyses were performed in R (version 2024.04.2 + 764), and a significance threshold of α = 0.05 was applied. Correlation between growth parameters of the evolved clones in the presence and absence of the antibiotic was evaluated by determining Pearson’s correlation coefficient.
To determine if the whole genome sequencing results varied between evolved clones isolated from co-culture and monoculture conditions, mutational profiles were converted into a binary presence/absence matrix, with each gene coded as “1” if mutated in an isolate and “0” otherwise. The Principal Coordinate Analysis (PCoA) was performed on the binary mutation matrix using the Jaccard distance metric in R via the vegan::vegdist function. Sample-level coordinates along the first two PCoA axes were used to visualize separation between isolates, whereas correlations of individual genes with PCoA2 were calculated to identify mutations contributing most strongly to the observed variation. To test whether co-culture and monoculture-evolved isolates differed at the group level, a permutation-based ANOVA was performed on the Jaccard distance matrix of mutational profiles. Isolates were distinguished either by growth condition (monoculture vs. co-culture) or by type of epistasis (negative, additive or positive). Pairwise co-occurrence of mutated genes was tested using Fisher’s exact test applied to all 2 × 2 presence–absence tables. Resulting P-values were corrected for multiple testing using the Benjamini–Hochberg FDR procedure. Odds ratios were used to quantify the direction and magnitude of association. Gene pairs with FDR < 0.05 were considered significantly co-occurring.
Results
Resistant mutants with higher levels of resistance and reduced growth rates are enriched in lineages evolving in co-cultures compared to monocultures
To determine how interspecific interaction affects resistance evolution, we performed evolution experiments for ~640 generations in single-species (monocultures) and two-species (co-culture) systems at sub-MIC of STR (1 μg/ml). At this concentration, STR has a stronger inhibitory effect on E. coli compared to S. enterica, consistent with the higher MIC of S. enterica (16 μg/ml) relative to E. coli (10 μg/ml), and supported by their respective dose–response curves (Supplementary Fig. 1). Despite this, across all the replicates, we observed co-existence throughout the course of the experiment (Supplementary Fig. 2). Emergence and fixation of resistant mutants were determined by plating 10^6^ cells per evolving lineage from different time points (~0, 160, 440, and 640 generations) on plates with STR concentrations of 16 μg/ml and 64 μg/ml. Resistant mutants of S. enterica were observed after 150 generations of growth in both conditions, and no difference was observed in the emergence or fixation of resistant mutants at 440 and 640 generations, either when lineages were plated together (Supplementary Fig. 3) or individually (Supplementary Table 1, Gen_400_: t_44_ = 1.66, P = .10; Gen_640_: t_44_ = 1.05, P = .3). To confirm that the resistant mutants resulted from adaptive evolution in our experiments, we measured the frequency of resistant mutants in S. enterica monocultures and S. enterica–E. coli co-cultures evolved for 640 generations without streptomycin. The significantly higher frequency of resistant mutants in our experimental conditions supports their origin through adaptive evolution (Freq_monoculture_: t_10_ = 3.89, P = .003; Freq_co-culture_: t_10_ = 2.32, P = .04; Supplementary Table 2). In contrast, during the course of our experiment, antibiotic resistance evolution in E. coli was limited, both in co- and monocultures, with no high-level resistant mutant isolated from co-culture conditions (Supplementary Table 3). Because we were interested in studying the evolution of high-level resistant mutants, we focused on the mutants from S. enterica.
To investigate if the resistance levels were different between the S. enterica-resistant mutants isolated from lineages evolving in co-cultures compared to monocultures, we isolated resistant mutants from both these conditions (24 from monoculture and 22 from co-culture conditions) after ~640 generations and measured levels of resistance for each mutant by determining the MIC of STR. Resistant mutants isolated from lineages evolving under co-culture conditions had a higher level of resistance as compared to mutants isolated from monoculture conditions (Fig. 1A, Wilcoxon rank-sum test P = .001, Supplementary Table 4). We then investigated whether the fitness costs of resistance mutations, determined in the absence of STR, were also different between these groups of mutants. We observed lower relative fitness (Fig. 1B, permutation-based ANOVA P < .0001), lower exponential growth rates (Fig. 1B, permutation-based ANOVA P = .0004), and lower stationary phase densities (Fig. 1B, permutation-based ANOVA P = .008) for mutants isolated from co-culture conditions as compared to monoculture conditions. Thus, resistant mutants with higher levels of resistance and higher fitness costs were more prevalent in S. enterica lineages evolving in co-culture compared to monoculture. We next investigated if the exponential growth rates and stationary phase densities for the mutants isolated from monoculture and co-culture conditions differed in the presence of sub-MIC of streptomycin. Similar to the results in the absence of antibiotics, the relative exponential growth rate was lower for the evolved clones isolated from co-culture conditions as compared to those isolated from monoculture conditions (Supplementary Fig. 5, permutation-based ANOVA P = .037); however, the stationary phase densities were no longer significantly different between the two groups (Supplementary Fig. 5, permutation-based ANOVA P = .68). We found a strong positive correlation between exponential growth rates and stationary-phase density for evolved resistant clones isolated from co-culture conditions when comparing antibiotic-free and streptomycin-containing environments (Supplementary Fig. 6A; Pearson r_exp. growth rate_ = 0.80, P < .0001; Pearson r stat.phase den. = 0.80, P < .0001). In contrast, for resistant clones isolated from monoculture conditions, a positive correlation was detected only for exponential growth rates measured across the two environments (Supplementary Fig. 6B; Pearson r_exp. growth rate_ = 0.93, P < .0001), whereas stationary-phase density showed only a modest correlation (Pearson r stat.phase den. = 0.39, P = .06).
To identify predictors of levels of resistance in these mutants, we fitted a multiple linear regression model with MIC as the response variable and growth rate, stationary-phase density, population size at cycles 1 (~8 generations) and 80 (~640 generations), and the frequency of resistant mutants (from which the clones were isolated) as predictors (Supplementary Table 5). The overall model was significant (F_5,40_= 5.34, P = .00075). MIC was significantly negatively associated with growth rate (β = −3.05 ± 1.02 SE, P = .0048) and positively associated with the frequency of resistant mutants (β = 2.19 × 10^5^ ± 9.61 × 10^4^ SE, P = .028). In contrast, stationary-phase density and population size at generations ~8 and ~640 had no significant independent effects on MIC (all P > .19).
Competition with E. coli results in the induction of the stringent starvation response and the upregulation of the aminoglycoside transferase enzyme
We examined how the growth rates and the stationary phase density for S. enterica and E. coli changed when these were grown together under sub-MICs of streptomycin (Supplementary Fig. 4). S. enterica showed a small reduction in growth rate in co-culture compared to monoculture (relative exponential growth rate_co/mono,_ = 0.95 ± 0.02, permutation-based ANOVA P = .02), and E. coli showed the opposite trend (relative exponential growth rate_co/mono,_ = 1.04 ± 0.03, permutation-based ANOVA P = .04). In contrast, the stationary phase density for both E. coli and S. enterica was lower in the co-culture compared to the monoculture conditions. This effect was more pronounced for E. coli where the relative stationary phase density_co/mono,_ was 0.28 ± 0.003 (permutation-based ANOVA P = .03), whereas the relative stationary phase density_co/mono_ for S. enterica was 0.85 ± 0.01 (permutation-based ANOVA P = .016).
Previous work has shown that interspecific competition results in upregulation of the stringent starvation response in S. enterica [42], resulting in an increased expression of the aminoglycoside transferase gene (aadA), which is known to confer resistance to STR [40, 42]. To determine if this was also true for the growth conditions examined here, we performed gene expression analysis for iraP (upregulated under stringent starvation response [43] and a stabilizer of the master regulator sigma S of the stress-response pathway [37, 44]) and aadA, under conditions of monoculture (grown alone) and co-culture (grown with E. coli). Our analysis showed that in S. enterica, the expression of iraP and aadA was upregulated both during exponential growth phase (Fold-changeiraP = 1.9, P = 0.15; Fold-changeaadA = 2.32, P = .11), and during stationary phase (Fold-changeiraP = 3.5, P = .04; Fold-changeaadA = 2.5, P = .01) under co-culture conditions as compared to monoculture, though this effect was only statistically significant in the stationary phase (Fig. 2A). To determine whether these observations were due to resource competition or the secretion of metabolites by E. coli, we measured the expression of both genes in S. enterica grown in E. coli’s spent media that had been restocked with nutrients. Expression values were normalized to those measured when S. enterica was grown in its own restocked spent media. For both genes, changes in expression were not statistically significant (Supplementary Fig. 7A; fold-change iraP = 1.6, P = .2; fold-change aadA = 1.5, P = .31), indicating that resource competition is the primary driver of the increased gene expression.
To determine how the aadA gene impacted fitness at sub-MIC of STR, we performed three different competition experiments. First, we competed an ΔaadA S. enterica mutant with the ancestral S. enterica at sub-MIC of STR (1 μg/ ml). We observed, as expected, that the ancestral S. enterica outcompeted the ΔaadA mutant. Thus, the selection coefficient of the ancestral S. enterica with respect to the ΔaadA mutant was 0.009 ± 0.003 gen^−1^ in the absence of STR, whereas it was 0.12 ± 0.007 gen^−1^ in its presence (Fig. 2B, permutation-based ANOVA P < .0001). Next, we competed the ancestral S. enterica and the ΔaadA mutant individually against E. coli (Fig. 2C). We observed that the deletion of aadA resulted in S. enterica not having the fitness advantage over E. coli in the presence of STR. Thus, the selection coefficient of ancestral S. enterica over E. coli was significantly different between conditions of absence and presence of STR, changing from −0.09 ± 0.02 gen^−1^ to −0.01 ± 0.013 gen^−1^ (permutation-based ANOVA P = .02). However, the selection coefficients for the ΔaadA mutant over E. coli were similar when comparing growth in the absence (−0.112 ± 0.006 gen^−1^) and presence (−0.113 ± 0.006 gen^−1^) of STR (permutation-based ANOVA P = .659). Finally, we competed the ancestral and ΔaadA S. enterica mutant against each other in the presence of E. coli at a STR concentration of 1 μg/ ml (Supplementary Fig. 7B). Similar to earlier results, the ancestral S. enterica outcompeted the ΔaadA mutant in the presence of STR (0.068 ± 0.002 gen^−1^). These results show that interspecific competition upregulates the stringent starvation response and the expression of the aadA gene, with the latter being important for fitness at low STR concentrations.
Resistance mutations selected in lineages evolving under co-culture conditions interact with aadA to produce a negative fitness epistatic effect
Our earlier results showed that interspecies competition with E. coli results in an increased expression of the aadA gene in S. enterica (Fig. 2A). To determine if this increased expression of aadA influences the nature of resistance mutations selected in S. enterica evolving at sub-MICs of STR, we investigated the fitness contribution of the enriched resistance mutations by themselves and in combination with the aadA gene. To this end, we deleted the aadA gene in 20 different resistant mutants, 10 each from lineages that evolved in co-culture and monoculture. We then competed the original 20 mutants and their respective aadA knockouts individually against the ancestral ΔaadA S. enterica mutant at sub-MIC of STR (1 μg/ml, Fig. 3A), giving us the values of s_res + aadA_ and s_res_, respectively (Materials and methods). Similarly, s_aadA_ was calculated by competing the ΔaadA S. enterica mutant against the ancestral S. enterica. If the effect of the resistance mutations with aadA is additive, then the sum of selection coefficients s_aadA_ and s_res_ should equal s_res + aad_ (Fig. 3A). Any deviation from additivity would imply epistatic interaction at the level of fitness (Fig. 3A). We observed that out of the 10 resistant mutants isolated from lineages evolving under monoculture conditions, eight demonstrated positive fitness epistasis (s_aadA_ + s_res_ < s_res + aadA_), whereas for the remaining two mutants, the effects were additive (Supplementary Table 6). In contrast, six out of the ten resistant mutants isolated from co-culture conditions demonstrated negative fitness epistasis (s_aadA_ + s_res_ > s_res + aadA_), four demonstrated additivity, and none demonstrated positive epistasis (Supplementary Table 6). Across both the treatments, we observed that resistant mutants isolated from co-culture conditions were more likely to demonstrate a pattern of negative fitness epistasis compared to those isolated from monoculture conditions (Fig. 3B, Wilcoxon rank-sum test P = .002). We next determined if the fitness advantage conferred by the resistance mutations, in the presence (s_res + aadA_) and absence (s_res_) of aadA, was different between the mutants isolated from co-culture and monoculture conditions. We observed that although the fitness of resistant mutants isolated from lineages evolving under co-culture conditions was significantly lower than that of mutants isolated from lineages evolving under monoculture conditions (Fig. 3C, permutation-based ANOVA P = .003), this difference became less pronounced in the absence of the aadA gene (Fig. 3D, permutation-based ANOVA, P = .7). Finally, we investigated how the resistance levels of the mutants changed when the aadA gene was deleted. For all the mutants, resistance was reduced substantially in the absence of the aadA gene (Supplementary Table 7), with there no longer being any difference between the resistance levels of the mutants isolated from the two conditions (permutation-based ANOVA P = .12). To examine whether variability in epistatic patterns was associated with other clone- or population-level traits, we performed a multiple linear regression with MIC, stationary-phase density, exponential growth rate, S. enterica population size at ~generation 640, and resistant mutant frequency as predictors. None of these variables significantly explained variation in epistatic patterns (F_6,13_ = 1.29, P = .32, R^2^ = 0.08). Individually, none of the predictors were significantly associated with epistasis (P > .25 for all, Supplementary Table 8).
Taken together, our results show that the observed differences between the resistance levels and the fitness effects of resistance mutations selected under mono- and co-culture conditions are mediated by the interaction of these mutations with the aadA gene.
Whole-genome sequencing reveals distinct clustering of evolved clones under co-culture and monoculture conditions
To investigate if the type of resistance mutations enriched at sub-MIC of STR were affected by interspecies competition, resistant mutants of S. enterica (18 each from mono- or co-culture conditions) were whole-genome sequenced. We first investigated if the total number of mutations accumulated differed between these conditions. Clones isolated from co-culture conditions had a total of 67 mutations compared to 98 mutations in clones from monoculture conditions (Fig. 4A). To account for confounding effects due to differences in S. enterica population sizes between monoculture and co-culture conditions, we normalized the total number of mutations observed in each case with the relative population sizes. On average, S. enterica population sizes in co-culture conditions were 0.7-fold lower than in monoculture conditions (Supplementary Table 9). After normalization, the number of mutations observed under both conditions became nearly identical (67 vs 68.6).
Whole-genome sequence analysis of resistant mutants selected at sub-MIC of STR in the presence (co-culture) and absence (monoculture) of interspecies competition. (A) Mutations per resistant mutant in 18 mutants selected in the absence and presence of interspecies competition. (B) Venn diagram illustrating the overlap of genes mutated in two or more clones isolated from co-culture and monoculture conditions. Genes involved in the electron transport chain, which are known to contribute to streptomycin resistance, are highlighted in bold and underlined. (C) Principal coordinate analysis (PCoA) of mutational profiles in evolved S. enterica isolates is shown. Mutational profiles of evolved clones from monoculture and co-culture conditions were converted into a binary presence/absence matrix (1 = gene mutated, 0 = gene not mutated) and analyzed using PCoA based on Jaccard distances. Sample coordinates along the first two PCoA axes are shown, with points colored by growth condition (monoculture vs. co-culture). Group-level differences between co-culture and monoculture isolates were assessed using a permutation-based ANOVA.
Mutations were identified in genes with similar functions across both conditions, particularly those involved in the electron transport chain, motility, adhesion, and chemotaxis (Supplementary Table 10 and Fig. 4B). Because each clone contained multiple mutations, we performed a principal coordinate analysis (PCoA) using a binary presence/absence matrix of all detected mutations. This revealed a clear separation between co-culture-evolved and monoculture–evolved isolates (permutation-based ANOVA P = .001; Fig. 4C). Because this separation occurred primarily along PCoA axis 2 (PCoA2), we next examined the gene loadings that drove this pattern. The top contributors toward the co-culture–associated region were rfbJ, rfbP, barA, dxs, and cheR, whereas genes such as hemL, fliA, rfbF, lrhA, and treC contributed most strongly toward the monoculture-associated region. Mutations in dxs (deoxy-D-xylulose-5-phosphate synthase) and in the hem genes (hemA, glutamyl-tRNA reductase; hemL, glutamate aminomutase) have previously been implicated in aminoglycoside resistance.
We then asked whether whole-genome mutation profiles could also explain the different fitness epistasis categories observed among the isolates. However, PCoA based on mutation presence/absence did not reveal significant clustering by epistasis class (permutation-based ANOVA P = .94; Supplementary Fig. 8). Likewise, pairwise co-occurrence analysis of mutated genes using Fisher’s exact test with Benjamini–Hochberg FDR correction identified no significant gene–gene associations (all FDR > 0.05; Supplementary Table 11).
We identified mutations in genes previously reported in our earlier study, where the S. enterica–E. coli two-species system was evolved without streptomycin [35]. These mutations occurred in genes related to motility (motA), cellular adhesion (genes in the fim operon), and chemotaxis (cheR), and likely reflect adaptation to the antibiotic-free environment.
Discussion
In this study, we demonstrate how interspecies interactions can influence the physiology of competing species, thereby altering evolutionary outcomes. By examining the evolution of antibiotic resistance in the gut pathogen S. enterica under sub-MICs of antibiotics, we found that interspecific competition with E. coli significantly shaped the fitness landscape, altering selection dynamics and driving distinct resistance phenotypes. Specifically, interspecific competition resulted in the selection of resistance mutations with higher levels of resistance and greater growth trade-offs. To explain this observation, we generated three lines of evidence. First, we demonstrated that the gene aadA, which has been shown to contribute to high-level resistance at sub-MICs of STR [34], is upregulated in the presence of interspecific competition with E. coli. Second, several resistance mutations selected in the presence of interspecies competition interacted with the aadA gene to give a negative fitness epistatic effect. Lastly, the resistance levels and fitness effects of the resistance mutations selected in the presence and absence of interspecies interactions converged when the aadA gene was deleted, suggesting that the aadA gene mediates the observed differences between the resistant mutants isolated from the two conditions. Based on these observations, a conceivable explanation for the enrichment of higher-level resistance mutations in the presence of interspecies competition is that the increased expression of the aadA gene shifts the susceptible ancestor upward along the concave fitness-resistance curve. This curve reflects a fitness landscape characterized by negative fitness epistasis, a pattern previously described for rifampicin resistance mutations [45] (Fig. 5). Consequently, small changes in fitness would necessitate substantial increases in resistance, driving the selection of higher level resistance mutations. In contrast, starting from a lower point on this curve, as would be the case in the absence of interspecific competition, would result in mutations causing relatively smaller changes in resistance to still have a substantial fitness advantage (Fig. 5). A potential alternative explanation for these observations could also be that larger phenotypic changes might be needed to overcome challenges faced due to interspecific competition, e.g. resource competition. However, this explanation seems unlikely because the fitness of the resistant mutants isolated from co-culture conditions is lower than those isolated from monoculture conditions (Fig. 3C).
A fitness epistasis model illustrating how interspecies interaction alters evolutionary outcomes by modifying cellular physiology is shown. Interaction with E. coli shifts the ancestral S. enterica along the fitness-phenotype curve due to increased aadA gene expression (blue circles). This shift positions the ancestral strain at a point where larger phenotypic changes (resistance in this case) are required for small changes in fitness. The starting point can vary depending on the strength of interspecies interaction (represented as different shades of blue). This is in contrast to the starting point of the ancestral S. enterica strain in the absence of interspecific competition (black circle), where small changes in resistance can result in large changes in fitness. Susceptible bacteria are shown as circles, and resistant bacteria are shown as triangles.
We observed that the interaction between resistance mutations selected in the presence of interspecies competition and the aadA gene resulted in either negative fitness epistasis or additive effects. These patterns likely reflect the changing intensity of interspecific competition for resources during the experiment, which in turn led to different physiological states and altered evolutionary outcomes. Additionally, because the two species in our system are competing for shared resources, the physiological responses and subsequent evolutionary outcomes may, in some cases, resemble those observed in a resource-poor environment even in the absence of any interspecies interaction. However, our system differs fundamentally from a standard two-stressor setup involving antibiotics and low nutrient availability. Here, the second stressor is biotic, another species with its adaptive trajectory and the potential for reciprocal ecological and evolutionary interactions. This creates a dynamic and context-dependent selective environment in which fluctuating resource competition, interference, and niche partitioning can unpredictably influence stress-response pathways. We propose that the variability observed in how de novo mutations interact with aadA, producing a range of fitness effects, is a consequence of these shifting biotic interactions. Such complexity and evolutionary variability are unlikely to arise under static abiotic stress conditions, underscoring the importance of studying resistance evolution in ecologically realistic, multi-species systems. Previous work has shown that interspecific interactions can constrain resistance evolution in E. coli. For example, in natural microbial communities, competition increased the costs of resistance mutations and thus raised the minimal selective concentration of antibiotics, and also provided a protective effect for susceptible strains [46]. Other studies have shown that resident microbial communities can reduce the selective advantage of resistance mutations, thereby suppressing resistance evolution [6, 47]. In contrast to these inhibitory outcomes, our study demonstrates that interspecies competition can also act in the opposite direction by inducing increased antibiotic tolerance, which in turn shifts the selective landscape to favour the emergence of higher-level resistance mutations. Together, these results highlight that interspecific interactions can either suppress or potentiate resistance evolution, depending on the ecological and physiological context.
Our study and previous research [42] highlight how susceptibility to antibiotics can be affected due to interspecies competition. Our study has extended this understanding by determining the evolutionary implications of this observation, highlighting the importance of considering community-level interactions in resistance management strategies, as competitive dynamics within microbial communities can alter the fitness landscape and promote the emergence of different resistant phenotypes. These results also suggest that interspecific competition, by upregulating stringent- and stress-response pathways, can alter evolutionary outcomes for a broad range of phenotypes in these microbial communities.
A previous study found that interspecies interactions had only a minor effect on evolutionary outcomes, contributing just 2%–5% to the overall evolutionary response [48]. Although this seems to contrast with our findings, the two observations can be reconciled by considering the environmental context: our evolution experiments were conducted under stressful conditions (i.e. in the presence of sub-MICs of antibiotics), whereas the previous study was performed in a non-stressful setting. Given that interspecies interactions are known to especially influence stress-related pathways [14, 15, 17], it is plausible that their impact on evolutionary outcomes becomes more pronounced under stressful conditions. Although further controlled experiments are necessary to test this hypothesis, taken together, our study and the previous one underscore the importance of examining interspecies interactions in stress-inducing environments.
The genomic patterns observed here suggest that interspecies competition not only alters the number of mutations acquired during low-level streptomycin selection but also shifts the evolutionary trajectories of S. enterica. Co-culture-evolved isolates were primarily shaped by mutations in genes rfbJ, rfbP, barA, dxs, and cydD, which are linked to LPS biosynthesis, cell membrane potential, central metabolism, and global regulation. These genes pushed isolates toward the co-culture region in PCoA space, suggesting that alterations in cell surface properties, metabolic pathways, and specific aminoglycoside resistance mutations may improve competitive fitness or antibiotic tolerance in the presence of E. coli. Monoculture-evolved isolates were pushed toward the monoculture region by mutations in hemL, hemA, fliA, lrhA, rfbF, and treC, reflecting adaptations in heme biosynthesis, motility, and carbohydrate utilization. Among these, mutations in the heme biosynthesis genes have also been linked to aminoglycoside resistance through alteration of membrane potential. Whereas hemL contributed strongly to monoculture clustering, the presence of hemA exclusively in monoculture isolates further supports the role of heme biosynthesis adaptations. Overall, these patterns indicate that interspecies interactions can result in the selection of a different spectrum of resistance-conferring mutations.
Understanding fitness-phenotype landscapes is a cornerstone of evolutionary biology, as it dictates the adaptive dynamics and evolutionary trajectories of populations [49, 50]. Several theoretical and empirical studies have explored the characteristics of these landscapes to establish general principles dictating the trajectories of adaptive mutations [51–54]. Our work builds on this foundation by integrating the role of interspecies interactions in shaping these landscapes, offering a significant development to previous research. Although prior studies have highlighted fitness epistasis both for beneficial [23, 55] and deleterious mutations [56], the focus of these studies has generally been to understand changes in fitness. In contrast, our study examines the mechanisms linking fitness landscapes to shifts in phenotypic distributions. Our study emphasizes that interspecies interactions are not only pivotal in shaping evolutionary trajectories but also in driving phenotypic shifts among the interacting species.
Supplementary Material
Supplementary_materials_wrag014
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Barraclough TG . How do species interactions affect evolutionary dynamics across whole communities? Annu Rev Ecol Evol Syst 2015;46:25–48. Available from: 10.1146/annurev-ecolsys-112414-054030. · doi ↗
- 2Brown JH, Whitham TG, Morgan Ernest SK. et al. Complex species interactions and the dynamics of ecological systems: long-term experiments. Science 2001;293:643–50. Available from: 10.1126/science.293.5530.643.11474100 · doi ↗ · pubmed ↗
- 3Chesson P . Mechanisms of maintenance of species diversity. Annu Rev Ecol Syst 2000;31:343–66. Available from: 10.1146/annurev.ecolsys.31.1.343. · doi ↗
- 4Thompson JN . Variation IN interspecific interactions. Annu Rev Ecol Evol Syst 1988;19:65–87. Available from: 10.1146/annurev.es.19.110188.000433. · doi ↗
- 5Johansson J . Evolutionary responses to environmental changes: how does competition affect adaptation? Evolution 2008;62:421–35. Available from: https://academic.oup.com/evolut/article/62/2/421/6854571.18031306 10.1111/j.1558-5646.2007.00301.x · doi ↗ · pubmed ↗
- 6Nair RR, Andersson DI. Interspecies interaction reduces selection for antibiotic resistance in Escherichia coli. Commun Biol 2023;6:331.36973402 10.1038/s 42003-023-04716-2PMC 10043022 · doi ↗ · pubmed ↗
- 7Collins S . Competition limits adaptation and productivity in a photosynthetic alga at elevated CO 2. Proc R Soc B Biol Sci 2011;278:247–55. Available from: 10.1098/rspb.2010.1173. · doi ↗
- 8Van Doorslaer W, Stoks R, Swillen I. et al. Experimental thermal microevolution in community-embedded daphnia populations. Clim Res 2010;43:81–9. Available from: http://www.int-res.com/abstracts/cr/v 43/n 1-2/p 81-89/.