Pig faecal community

Pig faeces were collected from four Cornish Black pigs without previous exposure to antibiotics in April 2016 on Healey’s Cornish Cyder farm (Penhallow, Cornwall, United Kingdom). Two hundred grams of faeces from each pig were pooled, mixed with 400 mL each of sterile glycerol and 1.8 g/L NaCl solution. The mixture was homogenized for 3 min in a Retsch Knife mill Gm300 (Retsch GmbH, Haan, Germany) at 2000 rotations per minute (rpm), filtered through a sieve (mesh size ~1 mm2), centrifuged at 500 rpm for 60 s at 4 °C and the liquid supernatant fraction was collected and frozen at −80 °C as the inoculum.

Pig faecal extract

Two hundred grams of faeces from each pig were pooled, mixed with 800 mL of sterile 0.9 g/L NaCl solution. The mixture was homogenized for 3 minutes in a Retsch Knife mill Gm300 (Retsch GmbH, Haan, Germany), at 2000 rpm, filtered through a sieve (mesh size ~1 mm2) and the liquid fraction was collected. The extract was then centrifuged (3500 rpm, 20 min, 4 °C), the supernatant collected and autoclaved (121 °C, 20 min). The autoclaved extract was centrifuged again (3500 rpm, 20 min, 4 °C) and the supernatant collected and used as a nutrient supplement.

Strains

The focal species, E. coli MG1655, was chromosomally tagged with a Tn7 gene cassette encoding constitutive red fluorescence, expressed by the mCherry gene [20] to ensure that E. coli can be detected and distinguished from other community members after competition based on red fluorescence. The Kn resistant, red fluorescent variant containing resistance gene aph(3′)-IIb encoding an aminoglycoside 3′-phosphotransferase was created previously [21, 22].

To create the Gm resistant mutant the strain was further tagged through electroporation with the pBAM delivery plasmid containing the mini-Tn5 delivery system [23, 24] for Gm resistance gene aacC1 encoding a Gm 3′-N-acetyltransferase [25]. Successful clones were screened for Gm resistance (30 μg/mL) and for the chosen clone a single strain growth curve in lysogeny broth (LB) was measured to ensure that the cost of the resistance gene was lower than 10% compared with the susceptible strain to ensure competitive ability.

Competition experiments

Competition experiments as well as initial growth of focal species strains were performed in 25 mL serum flasks with butyl rubber stoppers. As growth medium 10 mL of sterile LB medium supplemented with 0.1% pig faecal extract, 50 mg/L Cysteine-HCl as an oxygen scrubber and 1 mg/L resazurin as a redox indicator to ensure anaerobic conditions [26], was added to each reactor, heated in a water bath to 80 °C and bubbled with oxygen-free N 2 gas until the oxygen indicator resazurin turned colourless. After cooling down to 37 °C the appropriate concentration of antibiotic was added from a 1000× anaerobic stock solution.

Two isogenic pairs of the focal species, the susceptible, red fluorescent E. coli strain with either its Gm or Kn resistant counterpart, were competed across a gradient of six antibiotics concentrations (Gm [μg/mL]: 0, 0.01, 0.1, 1, 10, 100; Kn [μg/mL]: 0, 0.02, 0.2, 2, 20, 200). Strains as well as the community (100 μL of frozen stock) were grown separately under anaerobic conditions in triplicate reactors, replicates were combined, harvested through centrifugation, washed twice in 0.9% anaerobic NaCl solution and finally resuspended in 0.9% NaCl solution, adjusted to OD 600 0.1 (~107 bacteria/mL) and subsequently used in competition experiments. While the community was grown as an inoculum from the same frozen, homogenized stock, both subsampling and cultivation bias, inherent when growing an environmental community under laboratory conditions led to differences in original composition of the model community (Figs. S1, S2). When growing the community in isolation in the absence of antibiotics carrying capacity was reached after 18 h based on OD 600 readings in a spectrophotometer.

Isogenic strains were mixed at 1:1 ratio (community absent treatment), and that mix further added at 10% ratio to 90% of the faecal community (community present treatment). Approximately 106 bacteria of either mix were transferred to six replicate reactors of each of the antibiotic concentrations and grown at 37 °C with 120 rpm shaking for 24 h, which allowed growth up to carrying capacity. As a consequence of normalizing the total inoculum size the resulting inoculum size of the focal species in absence (~106 bacteria) and presence (~105 bacteria = 10% of total inoculum) of the community differed. A volume of 100 μL of each reactor was then transferred to a fresh bioreactor, grown for 24 h, transferred again for a final 24 h growth cycle and finally harvested for subsequent analysis.

Fitness assay

From each reactor after 3 days (T 3 ), as well as the inocula (T 0 ), a dilution series in sterile 0.9% NaCl solution was prepared and plated on LB and LB + AB (30 μg/mL Gm or 75 μg/mL Kn). For appropriate dilutions total and resistant red fluorescent E. coli colonies were counted under the fluorescence microscope. Plating of the susceptible strain on LB + AB plates further did not lead to any growth of spontaneous mutants. The relative fitness (ρ) of the resistant (r) compared with the susceptible strain (s) strain was subsequently calculated based on their individual growth rate (γ) throughout the competition experiment:

$$\rho = \frac{{\gamma _r}}{{\gamma _s}} = \frac{{{\mathrm{log}}\left( {10^6 \times n_r^{T_3}/n_r^{T_0}} \right)}}{{{\mathrm{log}}\left( {10^6 \times n_s^{T_3}/n_s^{T_0}} \right)}} \\ = \frac{{{\mathrm{log}}\left( {10^6 \times n_r^{T_3}/n_r^{T_0}} \right)}}{{{\mathrm{log}}\left( {10^6 \times \left( {n_{total}^{T_3} - n_r^{T_3}} \right)/\left( {n_{total}^{T_0} - n_r^{T_0}} \right)} \right)}}$$

Statistical significance testing (n = 6) was performed using a one-tailed t-test against neutral selection (ρ = 1) and ANOVA corrected for multiple testing to compare the relative fitness of different samples.

MIC assay

To assess the MIC of the susceptible and the resistant focal strain individually in the presence and absence of the microbial community reactors were inoculated with 105 of the focal bacteria and 106 bacteria from the community for the community present treatment. Triplicate reactors were grown overnight across a gradient of antibiotics. Concentrations were increased by 1 μg/mL (Gm susceptible), 2 μg/mL (Km susceptible), 25 μg/mL (Gm resistant) and 50 μg/mL (Km resistant), respectively. Reactors were then harvested and plated out on LB agar. Positive growth was scored as more than fourfold bacterial colonies growing on the plates compared with plating of the inoculum. The MIC was defined as the first concentration at which no positive growth was observed.

DNA extraction and sequencing

Bacteria from each reactor, as well as inoculum and original pig faecal community were harvested through centrifugation of 2 mL of liquid, followed by DNA extraction using the Qiagen PowerSoil kit as per the manufacturer’s instructions. The quality and quantity of the extractions was confirmed by 1% agarose gel electrophoresis and dsDNA BR (Qubit), respectively.

16S rRNA gene libraries were constructed using multiplex primers designed to amplify the V4 region [27]. Amplicons were generated using a high-fidelity polymerase (Kapa 2G Robust), purified with the Agencourt AMPure XP PCR purification system and quantified using a fluorometer (Qubit, Life Technologies, Carlsbad, CA, USA). The purified amplicons were pooled in equimolar concentrations based on Qubit quantification. The resulting amplicon library pool was diluted to 2 nM with sodium hydroxide and 5 mL were transferred into 995 mL HT1 (Illumina) to give a final concentration of 10 pM. Six hundred millilitres of the diluted library pool was spiked with 10% PhiXControl v3 and placed on ice before loading into Illumina MiSeq cartridge following the manufacturer’s instructions. The sequencing chemistry utilized was MiSeq Reagent Kit v2 (500 cycles) with run metrics of 250 cycles for each paired end read using MiSeq Control Software 2.2.0 and RTA 1.17.28.

Metagenomic libraries were created using the KAPA high throughout Library Prep Kit (Part no: KK8234) optimized for 1 μg of input DNA with a size selection and performed with Beckman Coulter XP beads (Part no: A63880). Samples were sheared with a Covaris S2 sonicator (available from Covaris and Life Technologies) to a size of 350 bp. The ends of the samples were repaired, the 3′–5′ exonuclease activity removed the 3′ overhangs and the polymerase activity filled in the 5′ overhangs creating blunt ends. A single ‘A’ nucleotide was added to the 3′ ends of the blunt fragments to prevent them from ligating to one another during the adapter ligation reaction. A corresponding single ‘T’ nucleotide on the 3′ end of the adapter provided a complementary overhang for ligating the adapter to the fragment ensuring a low rate of chimera formation. Indexing adapters were ligated to the ends of the DNA fragments for hybridisation on a flow cell. The ligated product underwent size selection using the XP beads detailed above, thus removing the majority of un-ligated or hybridized adapters. Prior to hybridisation the samples underwent six cycles of PCR to selectively enrich those DNA fragments with adapter molecules on both ends and to amplify the amount of DNA in the library. The PCR was performed with a PCR primer cocktail that anneals to the ends of the adapter. The insert size of the libraries was verified by running an aliquot of the DNA library on a PerkinElmer GX using the High Sensitivity DNA chip (Part no: 5067–4626) and the concentration was determined by using a High Sensitivity Qubit assay. All raw sequencing data have been submitted to ENA under study accession number PRJEB29924.

16S rRNA gene analysis

Sequence analysis was carried out using mothur v.1.32.1 [28] and the MiSeq SOP [27] as accessed on 07.08.2017 on http://www.mothur.org/wiki/MiSeq_SOP. Sequences were classified based on the RDP classifier [29]. Diversity was assessed based on observed OTUs at 97% sequence similarity. All sequences of the focal species E. coli were removed based on ≥99% sequence similarity. No sequences with this degree of similarity to the focal species were detected in the original faecal community. NMDS plots for the community were created based on the Bray–Curtis dissimilarity metric [30].

Further sample similarity was tested using analysis of molecular variance (AMOVA) a nonparametric analogue of traditional ANOVA testing. AMOVA is commonly used in population genetics to test the hypothesis that genetic diversity between two or more populations is not significantly different from a community created from stochastically pooling these populations [31, 32].

Metagenomic analysis

Metagenomic samples, as well as a reference genome for the focal species E. coli MG1655, were analysed using the ARG-OAP pipeline for antibiotic resistance genes detection from metagenomic data using an integrated structured antibiotic resistance gene database [33]. This resulted in the abundance of different resistance gene classes and subtypes within these groups normalized by 16S rRNA gene copy number. Antibiotics resistance genes detected in the E. coli reference genome were subtracted from the total number of hits per 16S rRNA gene copy based on the abundance of E. coli 16S rRNA gene/total 16S rRNA gene. Further, all antibiotics resistance gene numbers were normalized to the amount of pig faecal community 16S rRNA gene per total 16S rRNA gene copy.

Mathematical model

In order to illustrate possible mechanisms underlying the data for bacterial fitness in the presence/absence of the community for varying concentrations of Gm and Kn, we described our experimental setup mathematically. For this we first developed a discrete-time mathematical model for the growth of the susceptible and drug-resistant bacteria, s and r, respectively, in the presence or absence of the community, c.

Bacterial growth

The discrete-time model describing the growth of the bacteria i, i = s, r, c, is governed by the following iterative model

$$n_i^{t + 1} = n_i^t\left( {1 + \phi _i\left( {1 - g_i} \right)\left( {1 - f_i} \right)} \right),$$

where \(n_i^{t + 1}\) is the size of the population of strain i at time t + 1, and ϕ i is the maximum growth rate in the absence of competition and drug pressure. The reduction in growth due to density-dependent regulation/resource limitation, given as

$$g_i = \frac{{\mathop {\sum }

olimits_j e_{ij}n_j}}{{k_d}},$$

with k d as the carrying capacity and e ij being the competition coefficient, describing how much the presence of an allospecific strain j impacts the competitive fitness of strain i. The reduction in bacterial growth due to drug pressure, f i , is governed by a generalised logistic function

$$f_i = {\mathrm{min}}\left( {f_{{\mathrm{max}}},\frac{1}{{1 + e^{\alpha _i - \beta _i\ln c}}}} \right),$$

where c is the drug concentration (in μg/mL), α i and β i are the parameters describing the dose-response relationship for strain i, and f max = 0.9 is the maximum growth inhibition.

Model simulation and relative fitness calculation

Starting from an initially small number of bacteria in fresh medium, we ran the model for 30 generations, at which point the bacterial population had reached carrying capacity, and diluted the population accordingly. The bacteria were again allowed to grow for 30 generations before being diluted and grown for a final 30 generations. At this point we calculated the relative fitness of the resistant strain as

$$\rho = \frac{{\gamma _r}}{{\gamma _s}} = \frac{{{\mathrm{log}}\left( {10^6 \times n_r^{90}/n_r^0} \right)}}{{{\mathrm{log}}\left( {10^6 \times n_s^{90}/n_s^0} \right)}}.$$

Community-dependent change in drug resistance/susceptibility

The Kn data seem to suggest that the benefit of the drug resistant bacteria is reduced in the presence of the community at medium to high drug concentrations pointing towards a decrease in the susceptibility of the susceptible strain in a community context. We captured this scenario by making the dose-response parameters α s,r and β s,r explicitly dependent on the density of the community by increasing the resistance of susceptible strain, s, i.e.

$$\alpha _s\left( t \right) = \alpha _{s,0}\left( {1 + \frac{{1.3\,n_c^t}}{{n_c^t + 10^3}}} \right),$$

$$\beta _s\left( t \right) = \beta _{s,0}\left( {1 + \frac{{0.35\,n_c^t}}{{n_c^t + 10^3}}} \right),$$

where α i,0 and β i,0 are the time-independent dose-response parameters. The effect of density dependence is further illustrated (Fig. S3).

Parameter estimations

For each drug (Gm and Kn) we obtained a set of parameter values that resulted in a good overall fit between the model simulations and the data, where the data comprised the observed relative fitness for both sets of experiments (i.e. bacteria grown in the presence and absence of the community) for six different drug concentrations. To allow for logarithmic regression the non-antibiotic control was assumed as one order of magnitude lower than the lowest concentration used in the experiment. The parameter values were determined by minimising the root-mean-square error using an optimisation algorithm akin to simulated annealing [34]. The aim here was not to perform rigorous parameter estimation but rather to find a set of parameters that, given specific model constraints and assumptions, resulted in model behaviours that qualitatively agreed with both the observed dynamics over the repeated growth cycles and the empirically determined fitness values. In fact, our method failed to find a unique set of values that consistently gave the best fitting model, which suggests that the available data was insufficient to determine the global maximum. However, the qualitative relationships between individual parameters and between the parameters comparing the two antimicrobials were fairly consistent between model runs. Tables 1, 2 list the sets of parameters estimated for the two different antibiotics. These parameters allow changes in community and focal species densities to be estimated throughout the individual competition experiments based on start- and end-point measurements (Fig. S4).

Table 1 Model parameter values for gentamicin selection curves Full size table