Non-additive microbial community responses to environmental complexity

Selection and initial metabolic profiling of organisms

In order to maximize the chance of obtaining communities with diverse taxonomic profiles from different environmental compositions, the organisms selected were drawn from a number of bacterial taxa known to employ varying metabolic strategies. In addition, given the growing relevance of synthetic microbial communities to industrial and biotechnological applications^73,74,75,76, we chose to employ bacterial species that have previously been used as model organisms and have well-characterized metabolic capabilities. This criterion, paired with the availability of flux-balance models associated with a majority of these organisms, allows us to explore the metabolic mechanisms observed in our various experimental conditions with higher confidence. These selection principles resulted in a set of 15 candidate bacterial organisms (Acinetobacter baylyi, Bacillus licheniformis, Bacillus subtilis, Corynebacterium glutamicum, Escherichia coli, Lactococcus lactis, Methylobacterium extorquens, Pseudomonas aeruginosa, Pseudomonas fluorescens, Pseudomonas putida, Salmonella enterica, Streptomyces coelicolor, Shewanella oneidensis, Streptococcus thermophilus, and Vibrio natriegens) spanning three bacterial phyla (Actinobacteria, Firmicutes, and Proteobacteria, Supplementary Table 1, Supplementary Fig. 1a).

A microtiter plate-based phenotypic assay was used to assess the metabolic capabilities of each of the 15 candidate organisms. Each organism, stored in glycerol at −80 °C, was initially grown in 3 mL of Miller’s LB broth (Sigma–Aldrich, St. Louis, MO) for 18 h with shaking at 300 rpm at each organism’s recommended culturing temperature (Supplementary Table 1). To maximize oxygenation of the cultures and prevent biofilm formation, culture tubes were angled at 45° during this initial growth phase. Candidate organism Streptococcus thermophilus was found to have produced too little biomass in this time period and was grown for an additional 8 h. Each culture was then separately washed three times by centrifuging at 6000 × g for 2 min, removing the supernatant, suspending the pellet in 1 mL of M9 minimal medium with no carbon source, and vortexing or triturating to homogenize. The cultures were then diluted to OD₆₀₀ 0.5 ± 0.1 as read by a microplate reader (BioTek Instruments, Winooski, VT) and distributed into each well of three PM1 Phenotype MicroArray Plates (Biolog Inc., Hayward, CA) per organism at final OD₆₀₀ of 0.05 ± 0.01. The carbon sources in the PM1 plates (Supplementary Table 2, Supplementary Fig. 1b) were resuspended in 150 µl of M9 minimal media prepared from autoclaved M9 salts (BD, Franklin Lakes, NJ) and filter-sterilized MgSO₄ and CaCl₂ prior to inoculation. The cultures in each PM1 plate were incubated at each organism’s recommended culturing temperature with shaking at 300 rpm for 48 h. After this growing period, the OD₆₀₀ of each culture was measured by a microplate reader to quantify growth. To account for evaporation in the outer wells of the plates, which could yield in inflated OD readings, three ‘evaporation control’ plates with no carbon source were inoculated with bacteria at a final OD₆₀₀ of 0.05 and incubated at 30 °C for 48 h. The averaged OD₆₀₀ readings of these plates were subtracted from the readings of the bacterial growth plates to correct for evaporation. A one-tailed t-test was performed using these corrected OD₆₀₀ values to determine significance of growth above the value of the negative controls (p < 0.05). These final growth yields for the 15 candidate organisms are reported in Supplementary Fig. 2a, and aggregated analyses of the growth profiles of the organisms are reported in Supplementary Fig. 2b-d.

After this initial metabolic profiling, Streptococcus thermophilus was not included in any of the subsequent experiments as it displayed too low of a growth rate in the initial overnight growth phase and grew very minimally (no more than OD₆₀₀ 0.2) in fewer than 20% of the carbon sources in the PM1 plate. After inclusion in an initial mixed-culture experiment (com14, Supplementary Table 5), Salmonella enterica was also removed from future experiments due to its high levels of growth on all but one of the PM1 plate carbon sources. Its exclusion, meant to prevent its complete dominance in the subsequent mixed-culture experiments, resulted in a final set of 13 bacterial organisms.

For experiments involving a subset of the 13 organisms, the organisms were chosen to ensure they could be differentiated via agar plating. In the 3-species community experiment involving E. coli, M. extorquens, and S. oneidensis in combinations of 5 carbon sources (com3a), the organisms were selected based on their easily differentiable colony morphologies (Supplementary Fig. 21). In the second 3- and 4-species community experiments (B. subtilis, M. extorquens, P. aeruginosa, and S. oneidensis (com3 and com4)), selection was informed by differentiable colony morphology and additional metabolic criteria based on generalist-specialist relationships. Absolute growth yield on the Biolog PM1 plates was also considered, with the goal of including both high- and low- yielding organisms. Therefore, P. aeruginosa (high-yield generalist), B. subtilis (low-yield specialist), M. extorquens (low-yield generalist), and S. oneidensis (high-yield specialist) were selected.

Selection and combination of carbon sources

The carbon sources used in all experiments were selected from the 95 carbon sources contained in the Biolog PM1 Phenotype MicroArray Plate. This plate contains a variety of molecule types such as mono- and disaccharides, sugar alcohols, organic acids, amino acids, and more complex polymers (Supplementary Table 2, Supplementary Fig. 1b). Using the metabolic profiling experiments for each individual organism as a basis (Supplementary Fig. 2a), different criteria were established to choose the carbon sources used in each experiment depending on the desired complexity of the environment. An overarching criterion was that each experiment contain at least one sugar, one organic acid, and one amino acid to increase the possibility of synergistic interactions between carbon sources and resource-use orthogonality between the organisms.

For the communities grown in 5 carbon sources (i.e., com3a and com13a in D-Glucose, pyruvate, D-glcNAc, L-proline, and L-threonine), the following criteria were applied: D-glucose was selected as it resulted in the highest yield of each of the individual organisms, pyruvate was an organic acid with relatively high yields, D-glcNAc was a more complex sugar that resulted in varying individual growth yields, and L-proline and L-threonine were amino acids that resulted in generally high and low individual species yields, respectively. Communities were grown in all combinations of these five carbon sources (5 conditions of 1 carbon source, 10 of 2, 10 of 3, 5 of 4, and 1 of 5) for a total of 31 unique environmental compositions (Supplementary Table 4).

The carbon sources for the 32-carbon source experiments were selected based on the following criteria, in decreasing order of importance: carbon sources in which generalists individually displayed low levels of growth but favored at least one specialist (3 carbon sources), carbon sources that resulted in high-variance in growth yields across organisms (5 carbon sources), and carbon sources that resulted in low-variance in growth yields across organisms (7 carbon sources). These criteria were meant to increase the probability of observing more taxonomically diverse communities. The remaining 21 carbon sources were selected based on the total organism-specific yields they conferred (Supplementary Fig. 2a), with higher-yielding carbon sources being prioritized. Communities were grown in selected combinations of these 32 carbon sources (32 conditions of 1 carbon source, 16 of 2, 8 of 4, 4 of 8, 2 of 16, and 1 of 32) for a total of 63 unique environmental compositions (Supplementary Table 3). The selected combinations were chosen based on the Biolog growth yields of the 13 organisms under each carbon source, with the lowest-yielding carbon source (D-sorbitol) being paired with the highest (D-glucose) followed by the second-lowest and second-highest, etc.

Growth media conditions were assembled by resuspending each carbon source in distilled water to stock concentrations of 1.25 mol C/L and filter sterilizing using 0.2 µm membrane filter units (Nalgene, Rochester, NY). Carbon source stock solutions were stored at 4 °C for no longer than 30 days. A liquid-handling system (Eppendorf, Hamburg, Germany) was used to distribute the individual carbon source stocks in the appropriate combinations in 96-well plates. These prepared carbon source stocks were then sterilized using filter plates (Pall Corporation, Port Washington, NY) via centrifugation at 1500 × g for 2 min. These were then combined with M9 minimal medium (containing M9 salts, MgSO₄, CaCl₂, and no carbon) and filter-sterilized water to final working concentrations of 50 mM C in 96 deep-well plates (USA Scientific, Ocala, FL) for a total volume of 300 µl. This working concentration was selected such that all organisms would not grow beyond the linear range of OD₆₀₀ for biomass measurements.

Culturing of microbial monocultures and communities

After selecting the set of 32 carbon sources above, each of the 13 bacterial organisms was cultured independently in each individual carbon source prepared from stock solutions. Each organism was first inoculated from a glycerol stock stored at −80 °C into 3 mL of LB broth and incubated at 30 °C with shaking at 300 rpm for 18 h. The culture tubes were angled at 45° to prevent biofilm formation and to enable oxygenation of the cultures. The overnight cultures were then washed three times by centrifuging at 6000 × g for 2 min, resuspending in M9 medium without carbon, and vortexing and triturating if necessary to homogenize. The individual cultures were then inoculated in biological triplicate into the prepared media plates at final concentrations in 300 µl of OD₆₀₀ 0.05 ± 0.01 as measured by a microplate reader (BioTek-Synergy HTX). Additionally, a control plate was assembled containing one well inoculated with each individual organism with no carbon source to assess the decay/evaporation of the initial inocula, and one uninoculated well with 50 mM C of D-glucose to control for contamination. These monocultures were grown at 30 °C with shaking at 200 rpm for 48 h, after which their growth yields were quantified by transferring 150 µl to clear 96-well plates (Corning, Corning, NY) and reading absorbance values (OD₆₀₀). Biomass quantities are reported as the difference between the raw OD₆₀₀ readings of each sample and the corresponding OD₆₀₀ value of the negative control wells. Outlying OD₆₀₀ readings were removed by calculating Z-scores M for each individual measurement x_i using the median absolute deviation (MAD):

where (widetilde x) is the median across three biological replicates and 0.6745 represents the upper quartile of the normal standard distribution, to which the MAD converges. If the Z-score of an individual measurement exceeded 3.5, it was considered an outlier and removed. This process resulted in the elimination of 71 data points (out of 1248) across all organism monocultures. As with our Biolog phenotypic assay, a one-tailed t-test was performed using these corrected OD₆₀₀ values to determine significance of growth above the value of the negative controls (p < 0.05). These final growth phenotypes for the 13 organisms in stock solutions are reported in Supplementary Fig. 9 along with a comparison to their growth in the Biolog phenotypic assays.

For community experiments in combinatorial media, all communities were assembled using a bottom-up approach with each organism initially grown separately and diluted to the same starting concentrations before being combined. All individual organisms were prepared as described in the above paragraph, then combined at equal proportions and inoculated in biological triplicate into the prepared combinatorial media plates at final concentrations of OD₆₀₀ 0.05 ± 0.01 in 300 µl. Each community growth experiment additionally contained three control wells: one uninoculated well with 50 mM C of D-glucose to control for contamination, and two inoculated wells with no carbon source to assess the decay of the initial inocula. The communities were grown at 30 °C with shaking at 300 rpm for periods of 24 or 48 h before each passage. At each passaging step, the cultures were triturated 10 times to ensure the communities were homogenized and 10 µl were transferred to 290 µl of fresh media for the subsequent growth period. Yields of the cultures were quantified as described above, and processed using the aforementioned outlier removal procedure. This process resulted in the elimination of 8 data points (out of 192) for com3, 15 for com4, 10 for com13, and 17 (out of 768) across the four monocultures. A summary of the organisms, carbon sources, and culturing conditions for each community experiment is found in Supplementary Table 5.

Communities to be sequenced were centrifuged at 1500 × g for 10 min and the supernatant was removed. Cell pellets were stored at −20 °C until DNA collection was performed using a 96-well genomic DNA purification kit (Invitrogen). To harvest the DNA, each cell pellet was resuspended in 180 µl lysis buffer containing 25 mM Tris-HCl, 2.5 mM EDTA, 1% Triton X-100, and 20 mg/ml Lysozyme (Sigma–Aldrich). The samples were mixed by vortexing and incubated at 37 °C for 30 min, after which 20 mg/ml of RNase A (Invitrogen, Carlsbad, CA) and 20 mg/ml of Proteinase K (Invitrogen) with PureLink Pro 96 Genomic Lysis/Binding Buffer (Invitrogen) were added. The samples were mixed by vortexing and centrifuged after each reagent was added. The samples were incubated at 55 °C for 30 min, after which 200 µl of 100% ethanol (Sigma–Aldrich) were added. DNA from the resulting lysates was purified using a vacuum manifold according to the purification kit protocol (Invitrogen). Purified DNA was normalized to 15 ng/µl using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Waltham, MA). Library preparation was performed based on a paired-end approach developed by Preheim et al.⁷⁷, which targets the V4 region of the 16S rRNA gene with the primers U515F (5’-GTGCCAGCMGCCGCGGTAA) and E786R (5’-GGACTACHVGGGTWTCTAAT). Libraries were sequenced using an Illumina MiSeq platform at either the MIT BioMicroCenter, Cambridge, MA (com13a) or at QuintaraBio, Boston, MA (com13). QIIME2⁷⁸ was used to demultiplex raw files and produce FASTQ files for forward and reverse reads for each sample. DADA2⁷⁹ was used to infer sequence variants from each sample and a naïve Bayes classifier was used to assign taxonomic identities using a custom reference database with a 95% confidence cutoff. Reads not able to be assigned to a genus in our database were marked as ‘Unassigned’ and comprised 0.19% and 0.01% of all reads across all samples for com13 and com13a, respectively.

Communities to be assayed by agar plating were diluted by a factor of 10⁴ and spread on LB agar plates using autoclaved glass beads. Plates were prepared by autoclaving and distributing 18 mL of LB agar (Sigma–Aldrich) into petri dishes using a peristaltic pump (TriTech, Los Angeles, CA). Inoculated plates were incubated at 30 °C and imaged after 72 h using a flatbed scanner for colony counting. Colony counts for com3 and com4 were adjusted based on a standard dilution of the community members at equal concentrations measured by OD.

Significance between growth yields under differing environments was determined using a one-sided two-sample t-test with significance cutoffs of 0.05, 0.01, and 0.001. Species richness (S) is defined as the number of different organisms detected in a particular environment. Shannon entropy (H) is defined as follows:

where p_i is the relative abundance of organism i in a sample.

For hierarchical clustering analysis of communities, Spearman correlation coefficients were computed either for pairs of environments or pairs of organisms based on normalized vectors of species abundances. Hierarchical clustering was performed on the correlation coefficients using the ‘clustergram’ function in MATLAB, which calculated distances between clusters using the UPGMA method based on Euclidean distance.

Computation of epistasis scores

To quantify the non-additivity in how taxonomic diversity and balance could change in incrementally more complex environments, we first established definitions of expected values of species richness S and Shannon entropy H based on their values in lower-complexity conditions. Let a combined set of carbon sources AB be defined as the union of carbon source sets A and B. For carbon source sets A, B, and AB, the vectors of relative species abundances in each set are defined as V_A, V_B, and V_AB, respectively. The species richness values S for each set are therefore simply the number of positive species abundance values in each vector. Based on the organisms that survived in sets A and B, we establish the naïve assumption that at least as many organisms as survived in either environment will also survive in set AB. We therefore define S_{expected, AB} as max (S_A, S_B). This value is then compared to the experimentally observed species richness value S_AB using our epistasis score for species richness E_{S, AB}:

We use a similar expectation to calculate Shannon entropy epistasis E_H. Here, we first calculate the observed H for carbon source set AB as (H_{AB} = – {sum} {V_{AB}log _2V_{AB}}), and the expected H for carbon source set AB based on A and B as H_{expected, AB}=max(H_A, H_B). The epistasis score for Shannon entropy E_{H, AB} is therefore:

Consumer resource modeling

We employed a dynamical modeling framework to simulate the yields of arbitrary communities in increasingly complex environments, as well as the relative abundances of com4. The model builds upon Robert MacArthur’s consumer resource model^80,81,82 and a subsequent modification by Marsland et al.⁴³, which simulates the abundances of organisms over time as a function of resource availability, metabolic preferences, and exchange of secreted metabolites.

We define the individual species abundances as N_i for i = 1,…,S, and the resource abundances as R_α for α = 1,…,M. The key variable in calculating the abundances N_i is a stoichiometric resource utilization matrix C_iα, which defines the uptake rate per unit concentration of resource α by organism i (Fig. 1e). To calculate the growth of each organism on each resource, we multiply this matrix by a Monod consumption term R_α/(k_i,a + R_α) that simulates concentration-dependent resource depletion. Each consumed resource type α with abundance R_α is therefore consumed by organism i at a rate C_iαR_α/(k_i,a + R_α). These resources α are then transformed into other resources β by the organisms via a species-specific normalized stoichiometric matrix D_αβi. A fraction l of the resultant metabolic flux is returned to the environment to be made available to other organisms, while the rest is used for growth. In addition to these resource consumption terms, the species abundances are also defined by (i) a species-specific conversion factor from energy uptake to growth rate g_i, (ii) a scaling term representing the energy content of each resource w_α, and (iii) a quantity representing the minimal energy uptake required for maintenance of each species m_i. These terms are further defined in Supplementary Table 7. Taken together, the species abundances N_i over time are defined by:

The initial resource abundances R_α,0 and the initial organism abundances N_i,0 are first defined; then each resource is consumed in a manner dependent on the matrix C_iα, and converted to other resources based on the stoichiometric matrix D_αβi:

We selected the parameters for our equations based on experimental observations and quantities obtained from the literature (Supplementary Table 7). Our model contains no terms for resource replenishment or culture dilution, as all species were diluted every 48 h at a proportion defined by our experiments (10 µl of culture passaged into a total of 300 µl of fresh uninoculated media). Kinetic growth curves for com14 (Supplementary Table 5, Supplementary Fig. 22a) were used to estimate the orders of magnitude for the remaining parameters, based on the community reaching an average of approximately 2.4 × 10⁸ CFU/mL (OD₆₀₀ 0.3) within 20 h in 50 mM C of D-glucose. The conversion factor from energy uptake to growth rate g_i, as well as the energy content of each resource w_α, were set to 1 and 1 × 10⁸, respectively, in order to approximate this magnitude of growth. The Monod resource consumption half-velocity constant was set to 1 × 10⁴ g/mL for all resources in order to approximate the experimentally observed growth timeframe. Lastly, the minimum energy requirements for all organisms m_i were informed by the community yields at steady state and the leakage fraction l_α was set to 0.8 based on community simulations in Marsland et al.⁴³. These quantities are summarized in Supplementary Table 7.

We used our model to simulate the growth of communities containing S = 13, S = 3, and S = 4 arbitrary organisms. Here, values in C could only be nonzero if a randomly defined probability (P_{{mathrm{util}}}^{ialpha }) of an organism i being able to utilize a particular resource α was below a given threshold value θ_i. This threshold, representative of the fraction of resources usable by each organism, ranges from 0 to 1 and can be decreased to make an organism more of a resource specialist. Each community was simulated 50 times in each environment, so that the resource preference matrices C and the resource utilization probabilities (P_{{mathrm{util}}}^{ialpha }) could be randomly repopulated. This process allowed us to more effectively sample the large space of possible resource utilization matrices and obtain a clearer indication of how mean community yields changed in response to increasing environmental complexity. The environments were generated from a set of M = 32 arbitrary resources, which were combined in a scheme similar to that of com3, com4, and com13: 32 conditions with one resource, 16 conditions with two resources, and so on, up until one condition with all 32 resources.

We initially simulated our communities with a θ_i of 1 for all organisms i, as an initial null model for assessing how growth yields could vary with increasing environmental complexity (Supplementary Fig. 5a, e, i). These communities, which had high degrees of niche overlap (Supplementary Fig. 10) and distributions of yield epistasis E_Y centered at zero, were used as a baseline to quantify the yield epistasis distributions of our in vitro communities (Fig. 2d). We next carried out simulations with decreasing values of θ_i (and correspondingly decreasing degrees of niche overlap), which recapitulated the increases in yield we observed experimentally for com3 and com4 (Supplementary Fig. 5f-g, j-k). We then carried out simulations in which we defined each value θ_i as the proportion of nutrients able to be consumed by each organism in our monoculture assays, θ_mc (Supplementary Fig. 9a, Supplementary Fig. 5d, h, l). Two of these simulation results (that of a 13-species community with θ_i = 0.5 for all i, and of a 13-species community with θ = θ_mc) were used as CRM-A and CRM-B, respectively, in our main analysis of taxonomic diversity with generalists and specialists.

To test the effects of metabolic exchange, our CRM simulations also contained between 1 and 10 unique secreted metabolites, which could also be consumed by organisms according to randomly defined preferences in C. For a more realistic representation of metabolic conversion, these byproducts were matched with primary resources in conversion matrix D, which was randomly populated according to a transition probability of 0.25, meaning that a given metabolic byproduct had a 25% chance of being converted from a given primary resource. This matrix was normalized across each primary resource to ensure conservation of mass. Our results for yield, species richness, and Shannon entropy are presented as the average across all quantities of metabolic byproducts. The initial species and resource abundances were set to 6 × 10⁶ CFU/mL and 1.5 g/mL, respectively, to approximate the initial OD₆₀₀ of 0.05 and the initial resource concentration of 50 mM C of glucose used in our experiments. All communities were simulated over the course of 288 h with dilutions every 48 h (based on com3, com4, and com13 culturing timescale) with a timestep of 0.01 h.

Our calculation for the degree of niche overlap ρ of our simulated communities was based on a previous formulation of the metric designed for consumer resource models⁸³. This metric, bounded between 0 and 1, is 0 if the organisms do not compete for resources and is 1 if the organisms’ resource utilization profiles completely overlap. It is defined using the mean µ and variance σ² of the community resource utilization matrix C as:

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com