
==== Front
Microbiology (Reading)
Microbiology (Reading)
micro
mic
Microbiology
1350-0872
1465-2080
Microbiology Society

001489
10.1099/mic.0.001489
Research Article
Antimicrobials and AMR
Characterizing a stable five-species microbial community for use in experimental evolution and ecology
Castledine Meaghan 1m.castledine@exeter.ac.uk

Pennycook Joseph 2118226379@umail.ucc.ie

Newbury Arthur 1arn203@exeter.ac.uk

Lear Luke 1l.m.lear@exeter.ac.uk

Erdos Zoltan 1z.erdos@exeter.ac.uk

Lewis Rai 1r.bransgrove@exeter.ac.uk

Kay Suzanne 1s.m.kay@exeter.ac.uk

Sanders Dirk 1d.sanders@exeter.ac.uk

Sünderhauf David 1d.sunderhauf@exeter.ac.uk

Buckling Angus 1a.j.buckling@exeter.ac.uk

Hesse Elze 1e.hesse@exeter.ac.uk

http://orcid.org/0000-0001-6799-9670
Padfield Daniel 1*d.padfield@exeter.ac.uk

1 Centre for Ecology and Conservation, College of Life and Environmental Sciences, University of Exeter, Penryn, Cornwall, TR10 9FE, UK
2 University College Cork, Cork T12 K8AF, Ireland
DanielPadfield, d.padfield@exeter.ac.uk
There are no conflicts of interest.

Supplement: One supplementary figures and two supplementary tables are available with the online version of this article.

2024
19 9 2024
170 9 00148928 3 2024
09 8 2024
Copyright © 2024 The Authors
2024
https://creativecommons.org/licenses/by/4.0/ This is an open-access article distributed under the terms of the Creative Commons Attribution License. This article was made open access via a Publish and Read agreement between the Microbiology Society and the corresponding author’s institution.

Abstract

Model microbial communities are regularly used to test ecological and evolutionary theory as they are easy to manipulate and have fast generation times, allowing for large-scale, high-throughput experiments. A key assumption for most model microbial communities is that they stably coexist, but this is rarely tested experimentally. Here we report the (dis)assembly of a five-species microbial community from a metacommunity of soil microbes that can be used for future experiments. Using reciprocal invasion-from-rare experiments we show that all species can coexist and we demonstrate that the community is stable for a long time (~600 generations). Crucially for future work, we show that each species can be identified by their plate morphologies, even after >1 year in co-culture. We characterise pairwise species interactions and produce high-quality reference genomes for each species. This stable five-species community can be used to test key questions in microbial ecology and evolution.

coexistence
microbial community
model community
stability
http://dx.doi.org/10.13039/501100000270 Natural Environment Research Council NE/V012347/1, NE/S000771/1 Buckling Angus
==== Body
pmcData Availability

All data and code for recreating the analyses and plots from the processed datasets is available on GitHub: https://github.com/padpadpadpad/five_species_community_paper. The assembled genomes have been submitted to NCBI and are available under the following accessions: Achromobacter veterisilvae AB1 (CP148753.1), Ochrobactrum teleogrylli AB1 (JBBHKQ000000000.1), Pseudomonas fluorescens AB1 (CP148752), Stenotrophomonas sp. AB1(2024)(JBBHKO000000000.1), and Variovorax sp. AB1(2024)(JBBHKP000000000.1). The strains are also available for purchase from the NCIMB culture collection under the following accessions: Achromobacter veterisilvae AB1 (NCIMB 15570), Ochrobactrum teleogrylli AB1 (NCIMB 15569), Pseudomonas fluorescens AB1 (NCIMB 15567), Stenotrophomonas sp. AB1(2024)(NCIMB 15568), and Variovorax sp. AB1(2024)(NCIMB 15571).The raw sequencing of the original colony identification (Sanger sequencing) is available on the GitHub repository. The raw short- and long-read whole genome sequencing data can be made immediately available on request, as can the sequencing for the long-term colony identification (16s amplicon sequencing), and the code that creates the processed datasets.

Introduction

Controlled experiments on synthetic microbial communities have been used to test general ecological and evolutionary theories about - for instance - species coexistence [1], adaptive radiations [24], and ecosystem functioning [56]. The short generation times and large population sizes of bacteria make them ideal for testing ideas that are experimentally intractable in other model systems. Culturable bacteria are also easy to manipulate, allowing us to change species richness [57 8], abundances [910], and even interaction type [1112] and strength [13] in a systematic and increasingly high throughput way.

The majority of previous experimental work has been done on monocultures or in synthetic communities consisting of only a few species. However, an increasing number of studies have used more diverse microbial communities to study how species rich communities assemble and community dynamics through time to understand the maintenance of the high levels of microbial diversity seen in nature [8,1417]. Experiments on both simple and complex model communities are key to understanding the processes underpinning community assembly and coexistence, and there is an inherent trade-off between more detailed mechanistic understanding and complexity [18].

A key assumption of this previous work is that the community members stably coexist, meaning that the densities of species in the system do not show long-term trends [19]. Modern coexistence theory presents the ability to invade from rare (negative frequency dependence) which tests that each species is able to increase from low density when in the presence of the rest of the community, as a key test of coexistence [1920]. Experimentally establishing that model communities are stable and that species coexist is key to increasing the likelihood that the findings are relevant to natural communities (where high levels of diversity are maintained, composition is relatively stable through time [2123], and where patterns of diversity cannot be explained solely through neutral processes [2426]). Second, when a model community is not stable, it risks biasing our understanding of species interactions to strongly competitive and negative interactions. In contrast, natural communities are thought to be generally stable owing to ecological and evolutionary processes that diminish competition. Third, demonstrating stable coexistence and understanding the mechanisms of coexistence (e.g. equalising or stabilising [19]) would increase the generalisability of findings and allow future work to understand under what circumstances stability breaks down.

There are several proposed methods to measure invasion growth rate [27], but they generally involve testing whether each species (e.g. species A) has a positive growth rate after invading into a resident community (containing all other community members other than species A) at equilibrium [28]. However, this becomes more complicated to test as species diversity increases. For example, following the removal of a species, the resident community may destabilise or convert to an alternative stable state [29]. Explicit tests of coexistence in model microbial communities have been done when communities are simple and contain only strains of one [3032], two [33], or three species [7], but they remain scarce when species diversity goes beyond this [58 13 15 17], but see [34]). Most model microbial communities are performed in batch culture, where populations are regularly transferred into fresh microcosms [2], meaning that ‘equilibrium’ is never reached in a measurable sense. Consequently, microbiologists have used invasion-from-rare assays to test for negative frequency dependence and stable coexistence. In these assays, the community is inoculated into fresh media, but one species (the invader) is at 100-fold lower density than the other species (the resident community). If each species has a higher growth rate when rare, relative to the resident community (a relative invader fitness greater than 1) then the community shows negative frequency dependence and is likely to be stable [3135 36].

Here, we report the characterisation of a stable five species community that can be used for future experiments in ecology and evolution. We cultured a community pool of 46 isolates retrieved from the same soil sample over several weeks in lab medium (1/64 Tryptic Soy Broth (TSB)) and the final community composition was determined. As outcomes were highly repeatable, we took a single replicate community and conducted reciprocal invasion experiments to check coexistence. In addition, we measured species interactions using spent media and co-culture assays, tracked long-term community dynamics, and generated high quality reference genomes. Finally, we discuss potential questions that could be answered with this - or any - stable community when so much information about its constituent species is known.

Methods

Community disassembly experiment

We created a model microbial community consisting of 5 species that (dis)assembled for approximately 85 generations (Fig. 1a). Briefly, this community was formed by inoculating 46 bacterial strains – belonging to 27 genera (Table S1, available in the online version of this article) – that were isolated from an experimental compost community [37] into 6 mL growth medium (1/64 Tryptic Soy Broth (TSB) diluted with demineralised H₂O) in 25 ml glass vials with loosened plastic lids. Forty eight cultures were set up, incubated statically for 13 weeks at 28 °C. For the first 5 weeks no transfers were made, after which weekly transfers were done, passaging 1% (60 µL) of each community into fresh 6 mL microcosms every week. After this time, we plated all 48 communities onto KB (King’s Medium B) agar and found they converged on a similar composition of ~5 dominant colony morphotypes. We did 16S amplicon sequencing of these 48 end point communities (see Supplementary Information, Figure S1), and selected a single replicate of these to explore its stability.

Fig. 1. Community (dis)assembly experiment to form a stable community. (a) Schematic of the experimental design for the formation of the stable microbial community. (b) Five morphs dominated the community. (red) Achromobacter sp., the smallest colony, circular and uniform in size. (blue) Ochrobactrum sp. white/grey in colour and opaque. (green) Pseudomonas sp. is the largest colony, cream in colour and opaque. (purple) Stenotrophomonas sp. is typically larger than Achromobacter sp. although variable in size and has a stronger colour. (orange) Variovorax sp. has a spreading mat surrounding the colony and is variable in size. (c) Reliability of morphotypes for identifying the different species.

This community consisted of five species: Achromobacter veterisilvae AB1, Ochrobactrum teleogrylli AB1, Pseudomonas fluorescens AB1, Stenotrophomonas sp. AB1, and Variovorax sp. AB1 (see Methods and Results for sequencing and identification methods). Crucially, these species could be reliably identified by their colony morphology (Fig. 1b), allowing us to track each species in the community and re-isolate them for phenotyping and clone re-sequencing. A single clone of each member of this community was picked and grown in monoculture for 2 days in 6 mL of 1/64 TSB, shaking, in 25 mL glass vials with loosened plastic lids at 28 °C to achieve high cell densities. These samples were then cryogenically frozen at −70 °C in glycerol (final concentration: 25%). These clones were used to recreate the community for subsequent experiments, and are available from the NCIMB culture collection (see Data Availability Statement).

Sanger sequencing of 16S rRNA of each morphotype

To understand how well our characteristion of different morphotypes represent different species, we sequenced the 16S rRNA gene of multiple replicates of each morphotype. Of the five distinct morphotypes, we picked at least eight colonies across all 48 microcosms. Polymerase chain reactions (PCRs) were performed in 20 µL reactions containing 10 µL of DreamTaq Green PCR MasterMix (2X) (Thermo Scientific), 0.4 µL each of the 338F and 518R primers (10 µM) and 2 µL of 1 : 100 diluted culture. In total 51 samples were Sanger sequenced using the 515F primer (Eurofins Genomics).

The Sanger sequencing files were analysed in R, with R functions specified as Rpackage::function. Chromatograms for each sample were created using sangeranalyseR::secondary.peaks [38]. Low-quality bases were trimmed from the start and end of each file using sangeranalyseR::trimm.mott, which uses Mott’s modified trimming algorithm, with a cut-off of 40 to define a bad-quality score. Files were further filtered to retain only sequences that were longer than 100 base pairs long, had five or less secondary peaks (where there is more than one peak at a given position; a secondary peak was defined when the lower peak was at least a third as high as the higher peak), and had a mean quality score of >30. This left 43 trimmed and cleaned sequences that went forward to further analysis.

For each of the five identified morphs, we calculated a pairwise distance matrix from the DNA sequences using ape::dist.dna [39], with a phylogenetic tree being constructed and plotted using ape::njs and ape::plot.phylo. If all replicates of the same morph were attributed to the same level on the tree, a single consensus sequence for the morph was constructed using msa::msaConsensusSequence [40] after aligning multiple sequences of the same morph using msa::msa. In instances where there was more than a single branch of the tree, samples within a morph were split along branches to create multiple consensus sequences. Consensus sequences were aligned using msa::msa to look at similarities between morphs. Taxonomy was assigned to genus level using dada2::assignTaxonomy [41], which implements the RDP Naive Bayesian Classifier to assign taxonomy to the RDP database (v16) [42].

Invasion-from-rare assays

To explicitly test for coexistence, we conducted invasion-from-rare assays [4]. To assess whether all species were necessary for coexistence, we did assays where each species was inoculated from rare in every species by diversity combination. In total, this yielded six replicates of 75 treatments. Species were grown from frozen to achieve high cell densities as described above. Cell density (colony-forming units [c.f.u.]) was estimated from optical density (600 nm wavelength; OD600) using a species-specific calibration curve (Table S2 for equations) and normalised to 105 c.f.u. −1. Focal species were invaded at a 100-fold lower density (1.6 µL) than that of the resident community (160 µL in total). When there was more than one resident species, resident communities were created from equal ratios of each species. Cultures were incubated for 1 week, after which samples were cryogenically frozen at −70 °C in glycerol. Species abundances (c.f.u. mL−1) were estimated by plating cultures onto KB agar and incubating for 2 days at 28 °C.

Relative invader growth rate of the focal species was calculated as the ratio of estimated Malthusian parameters (m), mfocal: mcommunity where mcommunity is the total density change of all the other populations combined. m=ln(N1/N0)/t where N1 is the final density, N0 is starting density, and t is the assay time (7 days) [4]. We performed independent one-sample t-tests on each combination of focal species by resident community combination, with the null-hypothesis being that the mean relative invader growth rate equals 1. A relative invader growth rate above 1 would indicate that the invader had a faster growth rate than the resident community, demonstrating negative frequency dependence. This resulted in 75 different statistical tests; P values were adjusted using the false discovery rate (FDR) method [43].

Medium-term persistence assays

To test whether diversity impacts coexistence and persistence, and to quantify pairwise and indirect (higher-order) interactions, each species was also inoculated at equal density in every possible community combination, from monoculture to five species. Species monocultures were prepared with 20 µL of each species inoculated (total inoculum varied across diversity levels) into fresh microcosms for each community combination. Twelve replicates of each unique community combination were then incubated for either 1, 2, 3 or 4 weeks (three replicates each) without resource replenishment. These assays allowed time for extinctions to occur, to measure the effect of resource scarcity and secondary metabolites for cultures left for 3 or 4 weeks, and for us to calculate species interactions of the community. At the end of each incubation period, individual species abundance (c.f.u. mL−1) was estimated using the methods described above. Abundances from week 1 were removed as it was visually clear that equilibrium density had not been reached for most species.

We calculated pairwise species interactions (e.g. effect of species 2 on species 1, w1|2) as the abundance of species 1 in the presence of species 2 relative to that of species 1 grown in monoculture (no interspecific interactions): w1|2=N1|2 /N1. If w is greater than 1, then the interaction is positive, whereas a w<1 demonstrates competitive interactions, with a value of 0 demonstrating competitive exclusion. To get a broad perspective of the nature of species interactions, we calculated the average (the mean of ~9 replicates) for each pairwise species interaction.

We estimated indirect (higher order) interactions from the persistence assays using a modified version of a framework originally developed to help quantify higher order antibiotic interactions on Escherichia coli [44]. As our assays cannot determine whether the presence of a third, fourth, or fifth species changes competitor density (an interaction chain) or changes per capita competitive effects (higher order interaction), our method calculates indirect interactions generally. This method overcomes the issues of the additive model when single species have large impacts on fitness, such that the additive expectation can never be achieved (i.e. a population size <0) [4546].

For each community combination with more than two species (e.g. species 1, 2, 3), the predicted effect on the focal species, ŵ1|23, is the product of the independent pairwise interactions: ŵ1|23=w1|2 × w1|3. The indirect interaction (ii) of any community combination on focal species 1, ii1|23, was calculated as the deviation of the observed abundance relative to monoculture, w1|23, from that predicted from pairwise interactions: ii1|23 = ŵ1|23 - w1|23. Using this definition, if an ii is the same sign as the original prediction (ŵ1|23), then it indicates synergism: the combined presence of species 2 and 3 on species 1 is greater than either on their own. If an ii is the opposite sign of the original prediction, then it indicates non-additivity and potential buffering (the combined presence of species 2 and 3 on species 1 is less than either on their own). To make all synergistic indirect interactions positive (and buffering interactions negative), we multiplied indirect interactions of predicted negative interactions by −1. We looked at whether the type of indirect interaction (buffering or synergism) was related to the strength of the expected pairwise interaction using a linear model with expected pairwise interaction as the response and type of indirect interaction as the response. This was only done on predicted negative interactions as there were only seven predicted positive interactions (all belonging to Variovorax sp.) in the dataset.

The original Tekin method rescales interaction estimates and quantifies indirect interactions after accounting for all interactions at lower species combinations (for example, in four species communities the indirect interaction of ii1|234 includes taking into account the estimates of ii1|23, ii1|24, and ii1|34). Recent work has discussed whether pairwise interactions alone can predict species abundances, so we chose to calculate all indirect interactions relative to the prediction of only the pairwise interaction estimates.

Supernatant assays

Cell-free supernatants were used to establish cross-feeding assays in 96-well plates. To obtain the supernatants, three replicate populations per species were established by inoculating a single clone in 6 mL of 1/64 TSB in 30 mL glass vials. Following 2 days of growth at 28 °C, we centrifuged individual cultures at 3000 r.p.m. for 10 min, and filtered these through 0.22 mm syringe filters (Merck, Millex PES). We confirmed supernatants were cell free by plating these onto KB agar. For each focal strain we inoculated ~5 x× 104 cells from an overnight culture in triplicate wells, each containing either 200 µL of fresh 1/64 TSB or cell-free supernatant of either the same or a different species (n=3 per pairwise combination). Following 48 h of growth at 28 °C, serial-diluted cultures were plated onto KB agar to obtain final cell densities.

Growth rate per day was calculated as ln(N1/N0)/t, where N1 is final density, N0 is the initial density, and t is two (the number of days of growth). If the growth rate was negative, we replaced the value with 0.001 to prevent negative values of relative growth required for comparisons between species. Relative growth rate was calculated by dividing the growth rate in supernatant by the growth rate in fresh media. In this way, if the value is above 1, the supernatant has a positive effect on growth, and if it is below 1 it has a negative effect on growth. This method gives a measure of the intraspecific effect and interspecific effect of each species by comparing the relative growth in each species’ own supernatant to those in other species’. No formal statistics were done on the interactions of individual species, choosing instead to only describe the types and strengths of interactions. These measurements may help explain why we see coexistence, and can act as a benchmark against which to understand how interactions change in future experiments.

We looked at whether estimates of species interactions from supernatant assays and co-culture experiments matched both qualitatively and quantitatively. Quantitatively, we used standardised major axis (SMA) regression to look at the relationship between the two estimates for each pairwise combination. SMA was useful here as we have no expectation about one variable being the response or the predictor, and allows us to test for an expected slope of 1. Qualitatively, we looked at whether the same sign of interaction (positive or negative) occurred when calculating pairwise interactions using co-culture and supernatant assays.

Long-term stability

To understand long-term stability of the five-species community, 12 cultures were transferred weekly for 60 weeks. To start the cultures, monocultures of species were grown from the freezer stocks to achieve high cell densities as described above. Populations were set-up using 20 µL inoculum of each species’ culture normalised using the species-specific calibration curves and incubated at 28 °C [47]. Serial 100-fold dilutions (60 µL culture into 6 mL 1/64 TSB as described earlier) took place every week for over a year. Cultures were plated every 2 weeks to check for contamination. After 60 weeks (at least 400 generations), our plate counts and identification of morphotypes indicated that all five species remained present in all cultures. To confirm this, we sequenced the 16S rRNA gene of multiple replicates of each morphotype (species) across multiple replicates. We picked at least 30 colonies of each morphotype across all 12 populations, with 169 samples in total. DNA was extracted by following Qiagen’s standard DNeasy UltraClean Microbial protocol, and the concentration of gDNA in each sample was determined using QuBit. A conserved fragment from the hypervariable region of the 16 s rRNA gene was targeted using N501F and N806R primers and amplicon sequencing was undertaken by the Centre of Genomic Research (Liverpool, UK) using an Illumina MiSeq to create 2×250 bp paired end reads.

We processed and analysed the sequence data in R using the packages dada2 and phyloseq [48]. Following the standard full stack workflow, we estimated error rates, inferred and merged sequences, constructed a sequence table, removed chimeric sequences and assigned taxonomy. During processing, forward and reverse reads were truncated between 25 and 250 nucleotide positions due to poor quality scores. Assembled amplicon sequence variants (ASVs) were assigned taxonomy using the RDP database. We then aggregated OTUs at the genus level as we knew each of the five species belonged to a different genus. One sample was removed as it only had 135 reads, and one sample was removed because the dominant ASV (the identity of the clone) was <0.95, indicating that the colony likely consisted of multiple species.

Whole genome sequencing of each morphotype

To better understand the genomic features of each species, we did de novo hybrid assembly using both short- and long-read sequencing. DNA was extracted using a NEB monarch genomic DNA kit following the manufacturer’s instructions, apart from for Ochrobactrum spp. samples, which were vortexed with beads and had an extra proteinase-K treatment to aid lysis. Short read Illumina sequencing was done by the Exeter Sequencing Service on a HiSeq 2500 ran in standard mode. Long read sequencing was done using the PacBio HiFi platform on a Sequell II using one SMRT cell by the Centre of Genomic Research (Liverpool, UK) and quality processed by them.

Hybrid assembly was done using HybridSPAdes [49] with the option ‘-isolate’ specified. The assemblies were polished with the trimmed Illumina data using Pilon [50] under default settings. Circlator [51] was used to circularise the assemblies, which utilised the PacBio reads which were corrected using Canu [52] with the option ‘genomeSize’ set to the approximate size for each isolate (estimated from the blast results of each isolates 16S sanger sequence). All contigs<500 bp were removed, and a contig that was identified as the phage phiX - used as spike in for the short-read sequencing - was removed. Assemblies were annotated using prokka [53] with default settings. We ran breseq [54] to map the short reads back to each assembly to identify potential errors in the assembly and to estimate depth. Taxonomy was estimated using gtdbtk [55] which uses FastANI [56] to do whole genome alignment, and by using dada2 to assign taxonomy using 16S genes in the assemblies. The quality of the assemblies was assessed using CheckM2 [57].

Results

Final (dis)assembled community

From plating all 48 communities, we identified five distinct morphotypes that were present in all communities. Based on Sanger sequencing, the five colony morphotypes were identified as Achromobacter sp., Ochrobactrum sp., Pseudomonas sp., Stenotrophomonas sp., and Variovorax sp. (Fig. 1b). Overall we have high confidence that these morphs represent different species. For Pseudomonas sp., Stenotrophomonas sp., and Variovorax sp., all morphotype replicates were assigned to the same consensus sequence. For Ochrobactrum sp., 11/12 morphotype replicates were aligned to a single consensus sequence, with the remaining replicate assigned to the closely related Rhizobium sp. As both sequences are on the same branch of the global tree, this morph is unlikely to have been misidentified to any other of the four morphs. For isolate 1 (likely Achromobacter sp.), 4/8 replicates were assigned to Achromobacter sp., 3/8 were assigned to Stenotrophomonas sp. and one replicate was assigned to Pseudomonas sp. (Fig. 1C). This may lead to some underreporting of Achromobacter sp. and overreporting of Stenotrophomonas sp., but it is unlikely to qualitatively alter results of experiments that use this community. We chose a single replicate community and a single isolate of each colony morphology from which to do all future experiments.

Reciprocal invasion-from-rare experiments

We determined whether all species were key for coexistence by invading each species from rare in all species combinations. Species could invade, irrespective of resident community composition or diversity, in 71 of the 75 combinations (one-sample t-test: Padj <0.05, Fig. 2). Moreover, there was considerable variation in relative invader growth rate, ranging from a minimum of 1.12 to a maximum of 4.34. Of the non-significant combinations, three were when Achromobacter sp. was the invading species (including when it was invading into the other four species together; Fig. 2a), the other when Stenotrophomonas sp. was invading into a community of Achromobacter sp., Variovorax sp. and Ochrobactrum sp. However, although not significant, mean estimates of relative invader growth rate were still above 1.

Fig. 2. Ability of species to invade-from-rare in all combinations of species and diversity. Panels represent the invader species. Colours of pie points reflect the mix of the resident community. Black points are observed values of relative invader growth rate (n~6), pie points are the estimate of relative invader growth rate from each one-sample t-test, and error bars represent lower and upper 95% confidence intervals. Shaded regions separate different levels of resident community diversity; dashed, vertical line drawn where relative invader growth rate=1, values above which indicate the invading species can successfully re-establish from rare.

Interactions based on co-culture and spent media

From doing co-culture assays in all combinations of species diversity, we were able to explore the nature of species interactions between species pairs in the model community and also how these changed at higher levels of diversity, allowing us to estimate indirect interactions. Four of the five species (Achromobacter sp., Ochrobactrum sp., Pseudomonas sp., Stenotrophomonas sp.) had competitive interactions, growing worse in co-culture with other species compared to monoculture (Fig. 3b), with three of the species negatively impacted by all other species. Although Variovorax sp. negatively impacted the growth of all other species, it grew better in the presence of three of the other species, demonstrating asymmetric positive interactions (Fig. 3b). The only species which had a symmetrically competitive interaction with Variovorax sp. was Stenotrophomas sp.

Fig. 3. Estimated pairwise interactions of the model community from supernatant and co-culture assays. (a) Growth of each species was in the supernatant of each species (including itself) and fresh media. (b) Growth of each species was measured in co-culture and in monoculture. (c) Relationship between estimates of interactions estimated from supernatant assays and co-culture. In all plots, values below one indicate negative interactions. In (a) large white points is the estimate of intraspecific competition. In (a) and (b) large coloured points are the mean value of all interspecific interactions with 95% confidence intervals, and small coloured points represent individual interspecific interactions (the median value of three technical replicates). In (c) the solid fitted line is the best fit of a standardized major axis regression, the dashed lines demonstrate where interactions go from being negative to positive, and the inset plot demonstrates how often the interactions are qualitatively the same type using co-culture and supernatant assays.

In combinations of more than two species, we saw large variation in estimates of indirect interactions (Fig. 4a). Overall, Ochrobactrum sp., Pseudomonas sp. and Stenotrophomas sp. had on average negative indirect interactions, which given they had negative pairwise species interactions, indicates interactions in multispecies culture are more synergistic than would be expected by pairwise interactions alone. The opposite is true for Achromobacter sp., where we see on average buffering indirect interactions where interactions in multispecies culture are less than would be expected by pairwise interactions alone. Variovorax sp. has huge variation in indirect effects, indicating that interactions from co-culture are bad at estimating how it will do in multispecies culture.

Fig. 4. Estimated indirect interactions of the model community from co-culture assays. (a) Total indirect interactions across levels of diversity show lots of variation across all focal species. (b) Impact of the strength of negative and positive pairwise interactions on whether the indirect interaction is synergistic (stronger than expected) or buffering (less than expected). As negative and positive pairwise interactions get stronger, there is a tendency for the indirect interaction to be buffering. In (a) and (b) large coloured points are the mean value of all interspecific interactions with 95% confidence intervals, and small coloured points represent individual indirect interactions (the median value of three technical replicates).

When looking at the indirect interactions of all species together (Fig. 4), we found that 41% were synergistic (23 out of 55), with 59% being buffering, when compared to the expected effect of pairwise interactions. For three- and four-species communities, buffering was found 40% of the time, and synergistic 60% of the time, with this pattern switching when all five species were present (albeit only with five data points). It is important to note that the definition of synergism and buffering of indirect interactions is sensitive to the null model we are comparing against, which we set here simply as the impact expected from pairwise interactions alone [46]. If we included the indirect interactions of three species combinations into the estimated indirect interactions of four species, we may have seen more buffering [44]. We chose not to do this for simplicity in explanation and interpretation.

No focal species consistently had buffering or synergistic indirect interactions irrespective of species combination they were co-cultured with (Fig. 4a). Across all indirect interactions, we found that as pairwise interactions intensify, indirect interactions are more likely to be buffering (lower effect than expected), with synergistic indirect interactions more likely when the expected pairwise interaction is weaker (Fig. 4b). This pattern is statistically significant for expected negative interactions (ANOVA between models with and without indirect interaction type as a predictor: F1,47=4.51, P=0.039), and the pattern is the same for expected positive interactions, albeit a lack of expected positive interactions prevents us from doing formal statistical analysis.

The supernatant assays allowed us to get an estimate of intraspecific competition (white points in Fig. 3a), as well as interspecific competition, with coexistence expected when intraspecific competition is greater than interspecific competition. In four of the five species (Pseudomonas sp. being the exception, where no difference was detected) intraspecific competition was stronger than the interspecific effect (Fig. 3a). In 3 of the species the strength of intraspecific competition was at least twice that of the estimate of the lower 95% confidence interval of interspecific competition, and in Variovorax sp. intraspecific competition was twelve times stronger than interspecific competition.

Both co-culture and supernatant assays are commonly used to estimate species interactions, but little work has compared the two methods. We find a positive relationship between estimates of interactions from supernatant and co-culture (slope=1.46, 95%CI=0.86–2.45, Fig. 3c), but the model only explains 15% of the variation in the data, meaning that the fit between the supernatant and co-culture estimates is poor. However, qualitatively, both approaches characterized the interaction as the same sign (positive or negative) in 18 of the 20 interaction estimates (Fig. 3c, inset). In the two instances where they disagreed, wo|a and ws|o, one of the estimates was more or less neutral and therefore difficult to define as either positive or negative.

Long-term coexistence

To examine whether the five-species community was stable long-term, we re-sequenced and visually identified multiple clones of each morphotype across six communities that had been cultured for 60 weeks. Crucially, this will also tell us whether morphotype differences still exist after long-term culturing which is key for the use of model communities in experimental evolution studies.

Of the 167 samples that passed quality control, 160 (~96%) were correctly identified from colony morphology (the genus assigned by the 16S sequencing was the same as that assigned by the assessor). Of the seven that were mis-assigned, 2 Achromobacter sp. were identified as Ochrobactrum sp., 2 Pseudomonas sp. were identified as Variovorax sp., 2 Variovorax sp. were identified as Pseudomonas sp., and 1 Achromobacter sp. was identified as Pseudomonas sp. These mis-assignments do not match reliably with those from the initial Sanger sequencing but remain extremely rare so are again unlikely to qualitatively alter results of experiments that use this community. This confirms that the community has long-term stability and that the individual species can still be identified by their morphotype after this time.

Genomic composition

Using short- and long-read sequencing, we created high quality reference genomes for each of the five species. All of the genomes have very high completion (>99.9%) and very low contamination (<2%), with the worst assembly having only three contigs (Table 1). Reassuringly, taxonomic assignment from GTDBtk corresponds to the genus-level taxon assignment from the Sanger sequencing, but for 3 of the five isolates a species name is now available (Table 1). Using GTDBtk to assign taxonomy using genome similarity, 4 of the species show >98% similarity to previously sequenced genomes (Table 1). The isolate genomes have now been submitted to NCBI and have been given strain names: Achromobacter veterisilvae AB1, Ochrobactrum teleogrylli AB1, Pseudomonas fluorescens AB1, Stenotrophomonas sp. AB1, and Variovorax sp. AB1. Genomic elements associated with plasmids, prophages, secondary metabolites, defence systems, and antibiotic resistance are present in all species. These genomes can be used as high quality references to look at genetic changes in future experimental evolution experiments.

Table 1. Genome characteristics and taxonomic assignment of the five species genomes

Isolate	NCBI submission name	GTDBtk assignment	16 s copy no.	Genome size (bp)	No. of contigs	GC content	Total coding sequences	Completeness	Contamination	
 1	 Achromobacter veterisilvae AB1	 Achromobacter veterisilvae	 4	 65 49 044	 1	 0.67	 5991	 100.00	 0.98	
 2	 Ochrobactrum teleogrylli AB1	 Ochrobactrum_B teleogrylli	 1	 50 64 753	 3	 0.57	 4714	 99.97	 1.41	
 3	 Pseudomonas fluorescens AB1	 Pseudomonas_E fluorescens_BH	 2	 67 59 727	 1	 0.60	 6053	 100.00	 0.15	
 4	 Stenotrophomonas sp. AB1 (2024)	 Variovorax sp003019815	 1	 66 11 336	 2	 0.68	 6041	 99.99	 0.37	
 5	 Variovorax sp. AB1 (2024)	 Stenotrophomonas spp.	 1	 39 26 254	 2	 0.66	 3384	 100.00	 0.04	

Discussion

We characterized a model five species community natural community that can be used in experiments in microbial ecology and evolution. The community contains species with mostly competitive interactions, but Variovorax sp. demonstrates exploitative (potentially cross-feeding) interactions with three of the four other species. We found that estimates of interactions from supernatant assays and co-culture assays agreed qualitatively, but not quantitatively, with estimates of interactions in co-culture being stronger on average than those in supernatant. This suggests that contact dependent interactions (e.g. T6SS and biofilm formation) may also be important, but it also may just reflect methodological differences in the approaches. We found evidence of indirect interactions, but they became weaker (and more likely to be buffering) as expected pairwise interaction increased. This meant that in the full five species community, indirect interactions were relatively weak.

This community has numerous advantages as a model system for ecological and evolutionary studies. Sequencing has confirmed that each species has a distinct colony morphology that means it can be tracked and re-isolated from the other members of the community, and that these colony morphologies remain distinct through a year of long-term culture. This makes this community ideal for experimental evolution. Furthermore, we used invasion-from-rare assays to experimentally demonstrate species coexistence, and sequencing showed that the community is stable through long-term (>1 year) batch culture. We encourage - where possible - more explicit tests of coexistence in microbial community experiments.

Our model community is composed of diverse bacterial species all isolated from the same soil sample. Even so, this community is not designed to mimic a soil community. Instead the aim was to create a community that we could show empirically to be stable that could then be used to investigate general questions about ecology and evolution. By using the complex media Tryptic Soy Broth (TSB) at a (relatively) low concentration (1/64), we likely created conditions that promote coexistence through both stabilising (increased niche differences meaning intraspecific competition was greater than interspecific competition (Fig. 3a)) and equalising (diluted media reduced the fitness differences between species) mechanisms. The long-term stability of the community and the unique ability to re-isolate each species to test for evolutionary responses in phenotype make the community perfect for doing long-term experimental evolution experiments, a multi-species version of Lenski’s pioneering LTEE. The genome assemblies provide a reference database to track genetic changes through evolutionary time.

A recent editorial considered what features would be desirable for a model community to be widely adopted. This included having a functional output for the community, balancing complexity with feasibility, reproducibility, the ability to study the community over time, and ensuring the community in some way represents a natural community [58]. Our community achieves most of these, but we believe coexistence is a key desirable feature for a model community, and is essential to many of the other features such as including a temporal aspect and reproducibility. Moreover, demonstrating stability means model communities can more easily be used to test ecological and evolutionary theory where the assumption is often that the species coexist and the community is stable. In addition, the ability to manipulate a model community and re-isolate individuals is an invaluable feature for measuring evolutionary changes.

The utility of this five species community is demonstrated in the current and previous work by us, and its potential to test future ideas. For example, we have demonstrated how species become (mal)adapted after coevolutionary time [47] and how invaders and disturbances interact synergistically in their effect on resident diversity [59]. While the five species community as described only contains bacteria, we have also discovered three phages that each can infect one of the five species [60] and combined the model community with a broad host range plasmid which can be carried by all five species [61]. Moreover, we have used this community to demonstrate that fitness effects of plasmids shape the structure of bacteria–plasmid interaction networks [61], and that the cost of plasmid carriage can depend on community context [62] . Given its long-term stability and the distinct morphotypes of each species, this community is ideal to understand how microbial communities evolve and how this alters stability in the face of environmental stressors. Furthermore, its stability means this community is an ideal one to use to measure invasibility, or as a way to validate and improve Genome Scale Models (GSMs) that use genomic information to predict metabolic functioning and species interactions [63]. Environmental conditions such as primary substrate, pH, or temperature could be changed to see how well GSMs do under different scenarios.

supplementary material

10.1099/mic.0.001489 Uncited Supplementary Material 1.

Funding: This work was funded by NERC. D.P. was also supported by a NERC Independent Research Fellowship (NE/W008890/1). E.H. was supported by a UKRI Future Leaders Fellowship (MR/V022482/1). D.S. was supported in part by grant MR/N0137941/1 for the GW4 BIOMED MRC DTP, awarded from the Medical Research Council. A.B. was supported by NERC grants NE/V012347/1 and NE/S000771/1.

Author contributions: J.P., E.H. and A.B. conceived the study and designed the experimental work. J.P., E.H., D.P. and M.C. conducted the experiments. R.L., and S.K. did DNA extractions and Q.C. D.P., and M.C. analysed the data. D.P. did the bioinformatic analyses. D.P. and M.C. wrote the first version of the manuscript and all authors contributed to editing subsequent versions.
==== Refs
References

1. Barber JN Nicholson LC Woods LC Judd LM Sezmis AL et al Species interactions constrain adaptation and preserve ecological stability in an experimental microbial community ISME J 2022 16 1442 1452 10.1038/s41396-022-01191-1 35066567
2. Buckling A Craig Maclean R Brockhurst MA Colegrave N The Beagle in a bottle Nature 2009 457 824 829 10.1038/nature07892 19212400
3. Jousset A Eisenhauer N Merker M Mouquet N Scheu S High functional diversity stimulates diversification in experimental microbial communities Sci Adv 2016 2 e1600124 10.1126/sciadv.1600124 27386573
4. Lenski RE Rose MR Simpson SC Tadler SC Long-term experimental evolution in Escherichia coli. I. adaptation and divergence during 2,000 generations. Am Nat 1991 138 1315 1341 10.1086/285289
5. García FC Bestion E Warfield R Yvon-Durocher G Changes in temperature alter the relationship between biodiversity and ecosystem functioning Proc Natl Acad Sci U S A 2018 115 10989 10994 10.1073/pnas.1805518115 30297403
6. Sanchez-Gorostiaga A Bajić D Osborne ML Poyatos JF Sanchez A High-order interactions distort the functional landscape of microbial consortia PLOS Biol 2019 17 e3000550 10.1371/journal.pbio.3000550 31830028
7. Friedman J Higgins LM Gore J Community structure follows simple assembly rules in microbial microcosms Nat Ecol Evol 2017 1 1 7 10.1038/s41559-017-0109 28812620
8. Hu J Amor DR Barbier M Bunin G Gore J Emergent phases of ecological diversity and dynamics mapped in microcosms Science 2022 378 85 89 10.1126/science.abm7841 36201585
9. Acosta F Zamor RM Najar FZ Roe BA Hambright KD Dynamics of an experimental microbial invasion Proc Natl Acad Sci U S A 2015 112 11594 11599 10.1073/pnas.1505204112 26324928
10. Jones ML Ramoneda J Rivett DW Bell T Biotic resistance shapes the influence of propagule pressure on invasion success in bacterial communities Ecology 2017 98 1743 1749 10.1002/ecy.1852 28397255
11. Hammarlund SP Gedeon T Carlson RP Harcombe WR Limitation by a shared mutualist promotes coexistence of multiple competing partners Nat Commun 2021 12 619 10.1038/s41467-021-20922-0 33504808
12. Harcombe W Novel cooperation experimentally evolved between species Evolution 2010 64 2166 2172 10.1111/j.1558-5646.2010.00959.x 20100214
13. Ratzke C Barrere J Gore J Strength of species interactions determines biodiversity and stability in microbial communities Nat Ecol Evol 2020 4 376 383 10.1038/s41559-020-1099-4 32042124
14. Bittleston LS Gralka M Leventhal GE Mizrahi I Cordero OX Context-dependent dynamics lead to the assembly of functionally distinct microbial communities Nat Commun 2020 11 1440 10.1038/s41467-020-15169-0 32188849
15. Estrela S Sanchez-Gorostiaga A Vila JC Sanchez A Nutrient dominance governs the assembly of microbial communities in mixed nutrient environments elife 2021 10 e65948 10.7554/eLife.65948 33877964
16. Estrela S Sánchez Á Rebolleda-Gómez M Multi-replicated enrichment communities as a model system in microbial ecology Front Microbiol 2021 12 657467 10.3389/fmicb.2021.657467 33897672
17. Goldford JE Lu N Bajić D Estrela S Tikhonov M et al Emergent simplicity in microbial community assembly Science 2018 361 469 474 10.1126/science.aat1168 30072533
18. Blasche S Kim Y Oliveira AP Patil KR Model microbial communities for ecosystems biology Curr Opin Syst Biol 2017 6 51 57 10.1016/j.coisb.2017.09.002
19. Chesson P Mechanisms of maintenance of species diversity Annu Rev Ecol Syst 2000 31 343 366 10.1146/annurev.ecolsys.31.1.343
20. Grainger TN Levine JM Gilbert B The invasion criterion: a common currency for ecological research Trends Ecol Evol 2019 34 925 935 10.1016/j.tree.2019.05.007 31300168
21. Costello EK Lauber CL Hamady M Fierer N Gordon JI et al Bacterial community variation in human body habitats across space and time Science 2009 326 1694 1697 10.1126/science.1177486 19892944
22. Delgado-Baquerizo M Oliverio AM Brewer TE Benavent-González A Eldridge DJ et al A global atlas of the dominant bacteria found in soil Science 2018 359 320 325 10.1126/science.aap9516 29348236
23. Fierer N Jackson RB The diversity and biogeography of soil bacterial communities Proc Natl Acad Sci U S A 2006 103 626 631 10.1073/pnas.0507535103 16407148
24. Dumbrell AJ Nelson M Helgason T Dytham C Fitter AH Relative roles of niche and neutral processes in structuring a soil microbial community ISME J 2010 4 337 345 10.1038/ismej.2009.122 19924158
25. Liao J Cao X Zhao L Wang J Gao Z et al The importance of neutral and niche processes for bacterial community assembly differs between habitat generalists and specialists FEMS Microbiol Ecol 2016 92 fiw174 10.1093/femsec/fiw174 27543321
26. Ofiteru ID Lunn M Curtis TP Wells GF Criddle CS et al Combined niche and neutral effects in a microbial wastewater treatment community Proc Natl Acad Sci U S A 2010 107 15345 15350 10.1073/pnas.1000604107 20705897
27. Godwin CM Chang F Cardinale BJ An empiricist’s guide to modern coexistence theory for competitive communities Oikos 2020 129 1109 1127 10.1111/oik.06957
28. Spaak JW Schreiber SJ Building modern coexistence theory from the ground up: the role of community assembly Ecol Lett 2023 26 1840 1861 10.1111/ele.14302 37747362
29. Chesson P Updates on mechanisms of maintenance of species diversity J Ecol 2018 106 1773 1794 10.1111/1365-2745.13035
30. Pande S Merker H Bohl K Reichelt M Schuster S et al Fitness and stability of obligate cross-feeding interactions that emerge upon gene loss in bacteria ISME J 2014 8 953 962 10.1038/ismej.2013.211 24285359
31. Preussger D Giri S Muhsal LK Oña L Kost C Reciprocal fitness feedbacks promote the evolution of mutualistic cooperation Curr Biol 2020 30 3580 3590 10.1016/j.cub.2020.06.100 32707067
32. Zhao L Zhang Q-G Zhang D-Y Evolution alters ecological mechanisms of coexistence in experimental microcosms Funct Ecol 2016 30 1440 1446 10.1111/1365-2435.12611
33. Hammarlund SP Chacón JM Harcombe WR A shared limiting resource leads to competitive exclusion in a cross-feeding system Environ Microbiol 2019 21 759 771 10.1111/1462-2920.14493 30507059
34. Chang C-Y Bajić D Vila JCC Estrela S Sanchez A Emergent coexistence in multispecies microbial communities Science 2023 381 343 348 10.1126/science.adg0727 37471535
35. Le Gac M Plucain J Hindré T Lenski RE Schneider D Ecological and evolutionary dynamics of coexisting lineages during a long-term experiment with Escherichia coli . Proc Natl Acad Sci U S A 2012 109 9487 9492 10.1073/pnas.1207091109 22645336
36. Turner PE Souza V Lenski RE Tests of ecological mechanisms promoting the stable coexistence of two bacterial genotypes Ecology 1996 77 2119 2129 10.2307/2265706
37. Hesse E O’Brien S Tromas N Bayer F Luján AM et al Ecological selection of siderophore-producing microbial taxa in response to heavy metal contamination Ecol Lett 2018 21 117 127 10.1111/ele.12878 29161760
38. Chao K-H Barton K Palmer S Lanfear R sangeranalyseR: simple and interactive processing of sanger sequencing data in R Genome Biol Evol 2021 13 evab028 10.1093/gbe/evab028 33591316
39. Paradis E Schliep K ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R Bioinformatics 2019 35 526 528 10.1093/bioinformatics/bty633 30016406
40. Bodenhofer U Bonatesta E Horejš-Kainrath C Hochreiter S msa: an R package for multiple sequence alignment Bioinformatics 2015 31 3997 3999 10.1093/bioinformatics/btv494 26315911
41. Callahan BJ McMurdie PJ Rosen MJ Han AW Johnson AJA et al DADA2: high-resolution sample inference from Illumina amplicon data Nat Methods 2016 13 581 583 10.1038/nmeth.3869 27214047
42. Maidak BL Cole JR Lilburn TG Parker CT Saxman PR et al The RDP (ribosomal database project) continues Nucleic Acids Res 2000 28 173 174 10.1093/nar/28.1.173 10592216
43. Benjamini Y Hochberg Y Controlling the false discovery rate: a practical and powerful approach to multiple testing J R Stat Soc Ser B Methodol 1995 57 289 300 10.1111/j.2517-6161.1995.tb02031.x
44. Tekin E White C Kang TM Singh N Cruz-Loya M et al Prevalence and patterns of higher-order drug interactions in Escherichia coli NPJ Syst Biol Appl 2018 4 1 10 10.1038/s41540-018-0069-9 29263797
45. Tekin E Savage VM Yeh PJ Measuring higher-order drug interactions: a review of recent approaches Curr Opini Syst Biol 2017 4 16 23 10.1016/j.coisb.2017.05.015
46. Tekin E Yeh PJ Savage VM General form for interaction measures and framework for deriving higher-order emergent effects. Front Ecol Evol 2018 6 166 10.3389/fevo.2018.00166
47. Castledine M Padfield D Buckling A Experimental (co)evolution in a multi-species microbial community results in local maladaptation Ecol Lett 2020 23 1673 1681 10.1111/ele.13599 32893477
48. McMurdie PJ Holmes S phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data PLoS One 2013 8 e61217 10.1371/journal.pone.0061217 23630581
49. Antipov D Korobeynikov A McLean JS Pevzner PA hybridSPAdes: an algorithm for hybrid assembly of short and long reads Bioinformatics 2016 32 1009 1015 10.1093/bioinformatics/btv688 26589280
50. Walker BJ Abeel T Shea T Priest M Abouelliel A et al Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement PLoS One 2014 9 e112963 10.1371/journal.pone.0112963 25409509
51. Hunt M Silva ND Otto TD Parkhill J Keane JA et al Circlator: automated circularization of genome assemblies using long sequencing reads Genome Biol 2015 16 294 10.1186/s13059-015-0849-0 26714481
52. Koren S Walenz BP Berlin K Miller JR Bergman NH et al Canu: scalable and accurate long-read assembly via adaptivek-mer weighting and repeat separation Genome Res 2017 27 722 736 10.1101/gr.215087.116 28298431
53. Seemann T Prokka: rapid prokaryotic genome annotation Bioinformatics 2014 30 2068 2069 10.1093/bioinformatics/btu153 24642063
54. Deatherage DE Barrick JE. Identification of Mutations in Laboratory-Evolved Microbes from Next-Generation Sequencing Data Using breseq Sun L Shou W Engineering and Analyzing Multicellular Systems New York, NY Springer New York 165 188
55. Chaumeil P-A Mussig AJ Hugenholtz P Parks DH GTDB-Tk v2: memory friendly classification with the genome taxonomy database Bioinformatics 2022 38 5315 5316 10.1093/bioinformatics/btac672 36218463
56. Jain C Rodriguez-R LM Phillippy AM Konstantinidis KT Aluru S High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries Nat Commun 2018 9 5114 10.1038/s41467-018-07641-9 30504855
57. Chklovski A Parks DH Woodcroft BJ Tyson GW CheckM2: a rapid, scalable and accurate tool for assessing microbial genome quality using machine learning Nat Methods 2023 20 1203 1212 10.1038/s41592-023-01940-w 37500759
58. O’Toole GA We have a community problem J Bacteriol 2024 206 e00073-24 10.1128/jb.00073-24 38529952
59. Lear L Inamine H Shea K Buckling A Diversity loss from multiple interacting disturbances is regime-dependent Ecol Lett 2023 26 2056 2065 10.1111/ele.14325 37847646
60. Castledine M Newbury A Lewis R Hacker C Meaden S et al Antagonistic mobile genetic elements can counteract each other’s effects on microbial community composition mBio 2023 14 e00460-23 10.1128/mbio.00460-23 37022158
61. Newbury A Dawson B Klümper U Hesse E Castledine M et al Fitness effects of plasmids shape the structure of bacteria-plasmid interaction networks Proc Natl Acad Sci U S A 2022 119 e2118361119 10.1073/pnas.2118361119 35613058
62. Sünderhauf D Klümper U Gaze WH Westra ER van Houte S Interspecific competition can drive plasmid loss from a focal species in a microbial community ISME J 2023 17 1765 1773 10.1038/s41396-023-01487-w 37558861
63. Frioux C Singh D Korcsmaros T Hildebrand F From bag-of-genes to bag-of-genomes: metabolic modelling of communities in the era of metagenome-assembled genomes Comput Struct Biotechnol J 2020 18 1722 1734 10.1016/j.csbj.2020.06.028 32670511
