Formation of a constructed microbial community in a nutrient-rich environment indicates bacterial interspecific competition

The isolates Pseudomonas sp. GM17, Pantoea sp. YR343, and Sphingobium sp. AP49 were used to construct a synthetic microbial community. In previous work, these isolates were identified as the final surviving species after serially passaging a mixture initially composed of 10 bacteria in a nutrient-rich medium (7). Therefore, we hypothesized that a co-culture consisting of these three bacterial strains can form a stable community in the same medium environment. Meanwhile, competition may be the primary factor influencing the synthetic bacterial community, as all three community members could grow individually in the R2A medium without interdependency (7). An evaluation of pairwise interactions showed that the strain Pseudomonas sp. GM17 inhibits the growth of strain Sphingobium sp. AP49 (7). For validation, strains GM17, YR343, and AP49 were inoculated into R2A medium with equal proportions and were passaged for five dilution cycles. To assess the resulting populations of the community members, spread plating was used to assess colony-forming units of each species as they could be distinguished based on colony morphology. The succession and relative populations of the three-member bacterial community were also analyzed by quantitative PCR (qPCR). As assessed by either technique, the three strains co-exist throughout the serial dilution process.

The relative percentage for each species in the community was calculated to elucidate temporal variations in community structure at different passages. In the three-member community, Pseudomonas sp. GM17 dominates the cell population 24 h after inoculation and keeps its majority presence throughout the following five growth-dilution cycles (Fig. 1a and b). The community structure assessed by both plate counting (Fig. 1a) and the qPCR method (Fig. 1b) shows a similar relative distribution in the percentage of each organism upon community stabilization. When analyzed by qPCR, the relative percentage of strain Pantoea sp. YR343 in the stabilized three-member community appears to be larger than that of the result obtained by the plate counting method. Additionally, the qPCR result shows a large deviation from the plate counting method at time 0 h for the community with equal relative ratio of each species.

The growth abundance of each species in the initial passage (P0) was temporally probed to dissect growth behavior during the assembly process. The component organisms’ growth behavior when grown in monoculture, under the same conditions, was used for comparison. Both colony counting and qPCR results identified that the three species in co-culture reached different abundance values that spanned two orders of magnitude, by the end of the growth cycle (Fig. 1c). As the dominant strain in the co-culture, Pseudomonas sp. GM17 did not attain the same abundance level compared with its monoculture. When using the colony counting method, Pantoea sp. YR343 showed no statistical difference in growth capacity when compared with its growth in monoculture. In general, for the results obtained by qPCR, the genome equivalent of each species in the mixed culture was lower when compared with that of the corresponding monoculture. Among the three species, Sphingobium sp. AP49 displays a decrease in abundance of two orders of magnitude when compared with its growth in pure culture when assessed by both the colony counting and qPCR methods. The decrease in the abundance of each species in the constructed community, compared with that of the corresponding monocultures in the same medium environment, indicates competitive relationships within this three-member community. This competitive relationship is further supported by the growth curves of the total biomass in the 3-member community compared to that of the monocultures as assessed by optical density (Fig. S1). It was observed that the overall growth of the 3-member community was lower at the end of initial passage, in contrast to the growth of individual strains Pseudomonas sp. GM17 and Sphingobium sp. AP49.

Using bacterial abundance measurements obtained by plate counting and qPCR, the maximum growth rate (µ_max) was determined from growth curves for each species in either mixed or pure culture mode (Fig. S2 through S5). The µ_max of each strain measured by qPCR shows no statistical variation (P > 0.05) between growth in mixed and monoculture modes (Fig. 1d). This suggests that competition among these three species does not affect their exponential growth rates in a significant manner even though these bacteria do not attain the same growth abundances as their corresponding monocultures. When assessed by the plate counting method, Pseudomonas sp. GM17 and Sphingobium sp. AP49 also show no statistical variation in their maximum growth rates when comparing growth in mixed and pure culture modes. In contrast to the qPCR assessment, plate counting indicates that Pantoea sp. YR343 has a statistically significantly lower µ_max in the three-member community when compared with its growth rate in monoculture (Fig. 1d).

To further understand dynamic responses and the potential role of founder effects on the stable 3-member community, a set of communities inoculated at different initial species proportions were constructed using the serial dilution workflow. The absolute abundance of each species in the microbial community was analyzed by both plate counting and qPCR in 24-h intervals throughout the dilution process. Despite 1,000-fold changes in initial inoculum ratios, the bacterial community composition moved toward similar population ratios that were dominated by Pseudomonas sp. GM17 (Fig. 2). Although the final population ratio depends on whether plate counting or qPCR is used for the assessment, the final community composition shows no statistical difference across the different initial inoculum ratios when using the same analytical method. These results indicate that the final population ratio does not depend on the relative population level of the initial inoculum. Furthermore, regardless of which strain initially dominates the inoculum, a common population ratio results.

The process of community structure formation over the course of the serial dilution passages for three representative initial inoculum ratios is shown in Fig. 3. With varying species proportions in the inoculum, Pseudomonas sp. GM17 was always dominant after five passages. When strain Pseudomonas sp. GM17 was diluted 1,000-fold in the initial inoculum, it recovered to a dominant level in two growth-dilution cycles (Fig. 3a and b). When strain Pseudomonas sp. GM17 is the major species in the initial inoculum (1.0:0.001:0.001 of GM17:YR343:AP49), it maintained its dominance over the course of the dilution cycles. However, strains YR343 and AP49 maintained presence in the community and increased their relative population levels by more than an order of magnitude by passage No. 2 as determined by both plate counting and qPCR methods (Fig. 3c and d). When starting with an equal initial inoculum ratio, the final composition of the 3-member community was 1.0:0.058:0.015 (GM17:YR343:AP49) as determined by the plate counting method (Fig. 1a). This ratio was applied as an initial inoculum ratio for community assembly, and the community composition showed no variation in subsequent passages (Fig. 3e and f).

Other initial inoculum ratios including 1.0:0.001:1.0, 0.001:1.0:0.001, 1.0:1.0:0.001, and 0.001:0.001:1.0 of GM17:YR343:AP49 were also assessed for their impact on the resulting population structure of the bacterial community (Fig. S6 and S7). In all cases, when Pseudomonas sp. GM17 starts at a 1,000-fold dilution level, it recovers its dominance in the community by the end of passage 1. Similarly, regardless of whether Pantoea sp. YR343 and Sphingobium sp. AP49 start out as dominant or minority components, they return to similar relative ratios after just a couple passages. Both plate counting and qPCR assessment methods show the same population trajectories for the different members of the community. Meanwhile, no statistical differences among all evaluated initial inoculum ratios were found in the final population ratio. Additionally, the absolute abundances of the individual members in the communities were compared with their corresponding abundances when grown in monoculture (Fig. S8). As observed earlier when using a common starting inoculum ratio, decreased cell counts for the individual isolates result when grown in a community setting. This decrease could be attributed to competition among the community members. Negative interactions can also be discerned through a reduction in community productivity when compared to a model that represents the summation of the productivities of the individual species (27, 29). Our results indicate that the productivities of all three-member communities with various starting ratios were significantly lower than the predictions based on the null model (Fig. S9). Meanwhile, the measures of overall biomass for all three-member communities with different initial relative ratios were found to be lower at the end of the assembly process compared to the strain Sphingobium sp. AP49, which exhibited the highest growth in monoculture among all three community members in monoculture mode (Fig. S10). These observations are consistent with the competitive interactions detected when comparing the growth level of each specific species in either a community setting or as an individual culture. Collectively, these results indicate that a stable microbiome results when mixing together these three microbes and that the resulting population structure is independent of the relative ratio of the starting inoculum. Furthermore, competition within the members of the community decreases their final population levels.

The interspecies interaction coefficients for the gLV models fitted to colony-forming unit (CFU) counting results of community growth under different initial inoculum ratios imply that 85% of the pairwise couplings within the three-member communities are negative interactions. The gLV models derived from qPCR results of growth dynamics estimated that 77% of the pairwise interactions are negative (Fig. S11). These results suggest that the competitive relationship among the bacteria in the three-member community was not significantly changed by different initial inoculum ratios.

dFBA accommodates dynamic effects of the extracellular environment on microbial metabolism. Stoichiometric metabolic models of Pseudomonas sp. GM17, Pantoea sp. YR343, and Sphingobium sp. AP49 were used as compartments in the community model. The dFBA predictions for biomass growth of individual cultures generally agreed with the corresponding growth levels obtained from experimental data (Fig. S12). The final composition of the three-member community as measured by the qPCR method (1.0:0.126:0.013 of GM17:YR343:AP49) was used as the reference structure for community dFBA simulations. By changing the maximum uptake rates of limiting nutrients to a lower level, the three-member community dFBA models converged to the same composition by the end of passage 5, even when initial abundances of individual species differed by three orders of magnitude (Fig. 4; Fig. S13). Experimentally, although convergence happens at an earlier passage number, these results suggest that the community dFBA models are able to recapitulate the experimental observations by manipulating the uptake rates of limiting nutrients.

In order to suggest patterns of metabolic interactions in a nutrient-rich environment, the metabolite exchange fluxes between different species in community dFBA models over six-passage time intervals were captured. The simulations of the three-member communities with equal initial inoculum ratio predict that Pseudomonas sp. GM17 secretes xanthine (XAN), octadecanoic acid (ocdca), and putrescine that are consumed by Pantoea sp. YR343. In return, strain YR343 is proposed to excrete glycerol, which is consumed by strain GM17. Sphingobium sp. AP49 is predicted to uptake putrescine produced by strain GM17 (Fig. 5). The complex medium environment in dFBA models does not contain these metabolites, which are thus considered as potential interspecies metabolic exchanges in the three-member community. When modeled by dFBA in R2A medium as a monoculture, strain YR343 secretes glycerol as a byproduct of its growth, which is consistent with the model prediction in mixed cultures. Simulation of strain GM17 in monoculture also shows that it secretes XAN, ocdca, and putrescine as predicted by community dFBA models. Therefore, the secretions of the metabolites involved in interspecies exchanges by strain YR343 and GM17 are not emergent metabolic fluxes in the three-member community dFBA model under nutrient-rich conditions.

In the community model, Pseudomonas sp. GM17 is predicted to secrete four amino acids including alanine, methionine, cysteine and proline, while Pantoea sp. YR343 is predicted to excrete alanine and proline into the medium environment. Although Sphingobium sp. AP49 is predicted to consume these amino acids, they were not considered as absolute interspecies exchanges since the medium environment also contains these amino acids as nutrients for bacterial growth. Besides amino acids, strain GM17 is predicted to secrete thymidine and H⁺. Strain YR343 is predicted to secrete uracil, NH₄⁺, and H⁺, and strain AP49 is predicted to secrete niacin, inosine, and NH₄⁺ to the medium environment. Although different community members are predicted to have uptake fluxes of these compounds and inorganic ions, the medium environment contains all of them as unlimited substrates for bacterial growth. Several other metabolic byproducts were predicted to be secreted into the medium environment but without any uptake by other community members (Fig. 5).

The proteomes of the three bacterial species in the community at 48 h (Table S2) were compared with the proteomic results of the monocultures (Table S3) by performing Student’s t-tests to identify the differentially abundant proteins (Fig. 6; Table S4). In the pairwise comparisons between the proteins identified in the monocultures versus community samples, several enzymes involved in biosynthesis of putrescine, glycerol, and other amino acids as predicted by the dFBA model were identified to be significantly more abundant in the community compared to their monocultures. Although the majority of the identified proteins overlapped between the community and monoculture proteomes, a large number of proteins were uniquely identified in the community for strain GM17 and YR343. An opposite trend was observed for strain AP49, where the majority of the proteins were uniquely identified in the monoculture, which could be due to the relatively lower number of proteins detected for Sphingobium sp. AP49 than those of Pseudomonas sp. GM17 and Pantoea sp. YR343 in the community (Fig. S14; Tables S5 and S6).

The bacterial strains Pseudomonas sp. GM17, Pantoea sp. YR343, and Sphingobium sp. AP49 were previously isolated from the rhizosphere of Populus deltoides, and draft genome sequences were reported (7, 45). Bacterial monocultures were inoculated from R2A agar plates and grown in 10 mL of R2A medium in test tubes overnight at 30°C with shaking at 200 rpm for preparing seed cultures. The R2A complex medium was prepared as described in reference (46). For community construction, each species was normalized to create a bacterial population of 10⁷ CFU/mL, using a calibration curve between OD₆₀₀ and CFU/mL established for each strain (Fig. S15). An initial inoculum was prepared by mixing 300 µL of overnight cultures of the three strains at equal ratios (100 µL of each strain) in 9.7 mL of R2A medium and cultivated at 30°C with shaking at 200 rpm. Every 48 h, 10% (vol/vol) of the culture was transferred to a new test tube carrying 9-mL fresh R2A medium. Triplicate test tubes were subjected to serial growth experiments. At the 24- and 48-h time points of each growth cycle, both spread plate-based CFU counting on R2A and qPCR were used to assay the population of each species. The relative percentage of each species was computed as the read counts or copies for each strain divided by the total number of reads or copies at each time point.

To obtain improved quality genome sequences, from Pseudomonas sp. GM17, Pantoea sp. YR343, and Sphingobium sp. AP49, gDNA PacBio SMRT long-read sequencing was performed at SNPsaurus (Eugene, OR). Genomic DNA was extracted from 4 mL of overnight cultures grown in R2A medium. DNA extractions were performed using Qiagen DNeasy kit (Qiagen, Germantown, MD) following the manufacturer’s instructions. Flye v.2.8 was used to assemble and polish the sequenced genomes (47). The assembly quality was assessed using BUSCO v.4 (48). An average genome coverage of approximately 270-fold was obtained for each genome. The assemblies were annotated using Prokka with default parameters (49).

Bacterial cells were collected by centrifugation at 12,000 rpm for 15 min. Total DNA was extracted using a Qiagen DNA Isolation kit (Qiagen) as per the manufacturer’s instructions. The extracted DNA was used as a template for the qPCR process. The concentration of DNA templates was measured using a Nanodrop system (Thermo Scientific, Wilmington, DE). Each qPCR mixture (20 µL) was composed of 10 µL of SYBR Green PCR Master Mix (Bio-Rad, Hercules, CA), 4 µL of primer mixture (5 mM for forward and reverse primers), 5 µL of H₂O, and 1 µL of DNA template. The qPCR conditions were as follows: initial denaturation at 95°C for 5 min, followed by 45 cycles of 95°C for 30 s and 60° for 45 s. Quantitative PCR amplification was performed using a CFX96 Real-Time System (Bio-Rad) and monitored with the Bio-Rad CFX manager v.3.1 software. The primers used for qPCR analysis in this work are shown in Table S1. The standard curves were generated using triplicate 10-fold dilutions of the DNA extracted from each species. Species abundance was expressed as the number of genome equivalents calculated from qPCR amplification of genomic DNA.

The tri-culture community system was constructed under eight different initial inoculum ratios. In addition to an initially equal proportion, other initial inoculum ratios include 1.0:1.0:0.001, 1.0:0.001:1.0, 0.001:1.0:1.0, 1.0:0.001:0.001, 0.001:1.0:0.001, 0.001:0.001:1.0, and 1:0.058:0.015 of Pseudomonas sp. GM17:Pantoea sp. YR343:Sphingobium sp. AP49. The initial ratio of 1.0:0.058:0.015 is the observed, stabilized community ratio when starting with equal amounts of the three species in the initial inoculation. To prepare an initial inoculum, an overnight culture of each bacterium was normalized to 10⁷ CFU/mL by calibration curve between OD₆₀₀ and CFU/mL, and then one or two of the corresponding strains were diluted 1,000-fold to 10⁴ CFU/mL. One hundred microliters of each component was combined with 9.7-mL R2A liquid medium to attain a final volume of 10 mL. Triplicate samples were subjected to serial growth experiments. The population densities for each of the three strains were tracked over six growth-dilution cycles, and the community structures were determined by both plate counting and qPCR methods.

The maximum growth rate of each species in the mixed culture was determined from growth curves derived from species-specific CFU counting and qPCR protocols. For determining the maximum growth rate of each species in monoculture, 0.1 mL of each strain was normalized to 10⁷ CFU/mL and inoculated into 9.9 mL R2A medium. The number of cells was monitored by both CFU counting and qPCR methods to generate the growth curve. The maximum growth rate of each bacterial strain was calculated as the slope of the plot between ln(X/X₀) and the time during the exponential growth phase, where X is the cell mass at t and X₀ = initial cell mass. The increase in population size between t = 0 and 48 h, that is, N_t=48/N_t=0, was calculated using both CFU counting and qPCR results for cell abundance.

The gLV model was simulated using a time-series data set of microbial abundance of community members during serial dilution processes by Metagenomic Microbial Interaction Simulator software (50). The mean abundance of species at each time point over three parallel experiments was applied as input for gLV modeling, and the interaction coefficient calculated by gLV model was used to infer a positive or negative interaction type.

The genomes of the three strains were used to reconstruct sequence-specific metabolic networks using the Department of Energy Systems Biology Knowledgebase (KBase) “Annotate Genome/Assembly with RASTtk” and “Build Metabolic Model” applications (51–53). The metabolic models were all gap-filled in the R2A medium environment. The KBase narrative containing the genomes and metabolic models can be found at https://kbase.us/n/117918/17/, https://kbase.us/n/117919/16/, and https://kbase.us/n/117923/15/.

The different species were incorporated into COMETS in Python for community dFBA simulations by importing the corresponding metabolic models from KBase pipelines. The dFBA simulations for the tri-species communities with different initial inoculum ratios were all performed using the parsimonious algorithm (24, 54). Glucose and amino acids were considered as the limiting carbon sources in the R2A complex medium environment. The initial molar amount of glucose was calculated based on the glucose content in the R2A medium composition. The amino acids were considered as equally distributed in peptone, casamino acids, and yeast extract by weight, and they were exhibited as molar weight for their initial composition in the R2A medium. All other inorganic ions, gases, and growth factors were set as unlimited throughout all the simulations. The serial transfer of the community was simulated by COMETS using batch dilution function, and the dilution factor was set as 0.1, which is the same as the serial dilution experiment for the community assembly process. The initial biomass of each strain in the in silico community was defined to match the initial inoculum amount and ratios used in community assembly experiments.

The dFBA framework for modeling the three-member bacterial communities was calibrated to experimental data for the final community composition obtained by the qPCR method. To fit the simulated, stabilized community structure to the experimental results, the uptake rates of the carbon sources for each species in the simulations of communities were adjusted. The uptake rates of glucose and all amino acids for a community member were changed simultaneously to the same value. The limiting carbon source (glucose and amino acids) uptake rates of all three species in the community were first set at 10 mmol/gDW/h. The limiting carbon source uptake rates of strain Pantoea sp. YR343 and Sphingobium sp. AP49 decreased simultaneously from 10 to 1 mmol/gDW/h with a decrement of 1 mmol/gDW/h until reaching 1 mmol/gDW/h, and then decreased with a decrement of 0.1 mmol/gDW/h. The uptake and secretion values for each community model are listed in Tables S7 through S9. If the final biomass percentage of a species in the simulation was lower than 95% of its final biomass percentage after six passages in the experiment, its limiting carbon source uptake rates were returned to the last value and decreased again with a decrement of 0.01 mmol/gDW/h or 0.001 mmol/gDW/h. When the final biomass percentage of every species in the community model was in the range of 95%–105% of the biomass percentage of corresponding species in the stabilized community after six passages in the experiment, the simulation was then considered as successfully fitting to the experimental data. The default V_max of each strain in the community model for all other unlimited nutrients was set to 10 mmol/gDW/h. The initial cell dry mass was 3.9 × 10⁻⁷ g, and if a species was 1,000-fold diluted at the starting point, its initial cell dry mass was set as 3.9 × 10⁻¹⁰ g. The fluxes of metabolites involved in interspecies exchanges were extracted using the method in reference (24). The growth dynamics of the monoculture of each species was also simulated using the dFBA model, and the metabolite exchange fluxes between the bacterial cell and the medium environment were extracted for further comparison with those exchange fluxes in the community models.