Modulation of microbial community dynamics by spatial partitioning.

Microbial communities inhabit spatial architectures that divide a global environment into isolated or semi-isolated local environments, which leads to the partitioning of a microbial community into a collection of local communities. Despite its ubiquity and great interest in related processes, how and to what extent spatial partitioning affects the structures and dynamics of microbial communities are poorly understood. Using modeling and quantitative experiments with simple and complex microbial communities, we demonstrate that spatial partitioning modulates the community dynamics by altering the local interaction types and global interaction strength. Partitioning promotes the persistence of populations with negative interactions but suppresses those with positive interactions. For a community consisting of populations with both positive and negative interactions, an intermediate level of partitioning maximizes the overall diversity of the community. Our results reveal a general mechanism underlying the maintenance of microbial diversity and have implications for natural and engineered communities.

Introduction

Microbial communities are critical to natural ecological processes, such as biogeochemical cycling[1], animal and human health[2,3], and engineering applications[4,5]. Microbial community structure, meaning species identities and their abundance, is a primary feature that defines the functioning of microbial communities[6]. Along with internal factors, such as growth rate, death rate, and interactions, external factors, such as ecological factors and chemical environments also modulate microbial community structures[7]. However, our knowledge is still limited regarding what factors impact microbial community structures in a scalable and general manner and how they operate.

Survey-based studies of complex microbial communities using sequencing technologies provide large amounts of high-quality data and empirical insights[8,9] but causal and mechanistic links are often missing between external factors and community structure[10]. In contrast, controlled assembly of a few species can provide mechanistic interpretations since specific variables related to community structure can be manipulated. These studies have investigated the contributions of different factors that are biological[10-13], chemical[14,15], or physical[16-18]. However, how the learned insights scale up to more complex communities, where diverse interaction types and higher-order interactions may be present, is difficult to test and remains unclear[10].

Among these factors, spatial partitioning is ubiquitous yet mostly overlooked for microbial communities. Spatial partitioning describes the physical separation of a community into local communities. For example, the physical architectures of the gut[19], plant root[20], and soil[21] all partition microbial communities into distinct local communities that are separated to different extents (Fig. 1a). Due to the complexity of the physical architecture of microbial communities, the partitioning can be mostly complete, such as the microbiota in two different animals or the local microbial communities in two separate droplets. It can also be partial, resulting from the cell mobility or dispersal across local environments[18] or diffusion of signaling molecules across local communities[16].

Fig. 1.

Spatial partitioning and simulation framework

a. Partitioning of microbial communities in natural habitats. The structures of natural habitats including colons, soil, and droplets prohibit global interactions by physically separating a global community into local communities that only contain subsets of the global populations.

b. Depiction of the partitioning concept. At t0, a community of interest is sampled into completely isolated local environments. Because of partitioning, many local communities are subsets of the global community, where some populations and interactions are missing. We use the number of local communities (N) to quantify the degree of partitioning. A partitioning of N = 5 is shown.

c. The impact of spatial partitioning on local and pooled communities. Communities in each local environment grow independently and the final local community structures are quantified at tf. Local communities at tf are pooled to determine the final global community structure.

d. Modeling growth of local communities. Each equation describes temporal dynamics a population (X) by accounting for its logistic growth, positive interactions (by removing stress), and negative interactions.

e. Three sets of model parameters. Stress (δ), positive interactions (β), and negative interactions (γ) are shown in matrix forms. Fitness is equivalent to 1 – δ.

f. A typical temporal dynamic of a local community simulated using our model. The simulated time courses are bounded for any parameter combinations that satisfy: non-negative initial condition (X ≥ 0 at t0), δ ≥ 0, , and . The time course shown is generated using a 10-population community, with an interaction matrix generated from: δ ∈ [0, 1.5], β = [0, 5], and γ ∈ [0, 0.5].

In the simplest case, where partitioning is complete, local environments each consists of only a subset of all members and partitioning restricts interactions within local communities. In general, spatial partitioning reduces the overall strength of interactions in the global microbial communities and lowers the number of interacting species for each individual member[22]. Moreover, the type of interactions experienced by a member can vary drastically depending on the random assembly of local environment[23]. In other words, spatial partitioning can modulate the dynamics of a microbial community by globally modulating the type and strength of interactions experienced by each member. This emphasis on interactions, derived from studying microbial communities, differs substantially from research in multicellular organisms, which places much greater emphasis on dispersal between local communities, abiotic factors, and neutral dynamics[24]. Focusing on interactions therefore has potential to contribute to the historically organismal-level study of spatial effects on local and global community diversity.

It is yet unclear whether the effect of spatial partitioning is highly system specific or whether it follows general rules. Beyond the challenges of distilling causal mechanisms and general rules, defining spatial partitioning in a relevant and quantitative manner is also challenging. To address this question and overcome these challenges, we first established a theoretical framework to explain the mechanisms by which spatial partitioning affects community structure. Based on the theoretical framework, we formed a hypothesis that spatial partitioning reduces biodiversity for negative interaction dominated community and increases biodiversity for positive interaction dominated community, and biodiversity peaks at an intermediate partitioning level for communities with both positive and negative interactions. We then tested our hypothesis using precisely controlled top-down experiments of simple communities and scaled up to complex natural communities. The ability to control microbial community structures through modulation of spatial partitioning can address a wide range of challenges we face with natural and engineered microbial communities for ecological, medical, and engineering purposes.

Results

A theoretical framework to model spatial partitioning

Consider a microbial community residing in a global environment. In the absence of partitioning, all interactions between members are retained. If the environment is divided into N equal-sized local environments, the community members will be allocated to these local environments. We simplify the process by assuming seeding follows a Poisson distribution where the Poisson parameter λ dictates the average number of cells in one local community (Fig. 1b). On average, the local communities have the same total number of cells subject to random variations: the relative variation in this number increases as the total cell number in the overall community decreases. The parameter N measures the level of partitioning: a larger N corresponds to higher partitioning. We assume local communities are completely isolated from each other; as such, no interactions, dispersal, or migration occur across different local communities. For a sufficiently large N, some local communities will only have a subset of the members in the overall community. As a result, some members will experience fewer interactions in comparison with when they reside in an unpartitioned environment. With extreme partitioning where each local environment contains at most one cell, all interactions between members are eliminated. Thus, when interactions are considered, spatial partitioning, at its core, blocks interactions across local communities.

After partitioning, each local community grows separately, and the ensemble of local communities captures the global community dynamics (Fig. 1c). We model the temporal dynamics of each local community using a set of ordinary differential equations (ODEs) (Fig. 1d), where δ describes the intrinsic death rate of each species and β and γ describe positive and negative effects of species on others, respectively (Fig. 1e). 1/β is the density of a partner species to reduce death rate (δ) of this strain by 50%. When β increases, the positive interaction strength is higher (i.e. lower density of partner species is required to halve the stress) and when γ increases, the negative interaction strength is higher (i.e. lower density of partner species is required to reach the same death rate). Thus, we use strength of interactions to indicate the magnitude of β and γ.

In contrast to the classic general Lotka-Volterra (gLV) model formulations that can generate unbounded growth with some parameter ranges[25,26], our formulation generates bounded dynamics with the entire space of defined parameter domain (Fig. 1f, Supplementary Fig. 1). Our model only incorporate competition in the pairwise γ term, instead of assuming all species have background negative interactions. Thus, our model formulation accounts for different types of negative interactions (including competition) explicitly. Note that our simulation condition implies that the interaction length scale is larger than or similar to the scale of the local environment and that the local communities are well-mixed. Another assumption is that the interaction logic and strengths are preserved regardless of the initial partitioning.

Through its effect on local community membership at t = 0, the spatial partitioning level N modulates the growth of populations and the pooled global community structure after growth at t = t. This seeding and growth process capture the critical aspects of the temporal evolution of microbial communities in nature and in the laboratory setting. Examples include the inoculation and growth of communities in germ-free animal models[27], infant guts[28], cheese microbiome[29], and so on. The pooling process is analogous to the mixing of samples when quantifying microbial communities in natural habitats.

Emergence of biphasic dependence of biodiversity

By varying N, we focus on examining how varying spatial partitioning affects biodiversity (Fig. 2a), which is a key parameter that influences community stability[30], function[31], and evolution[32,33]. Here we primarily use the inverse Simpson index (referred as Simpson index in the following texts) as the metric to serve as the proxy for effective number of members. where p is the relative abundance of species i and the community has M number of species. The index reaches the maximum, which is M, when all species have the same relative abundance. Whereas the index reaches the minimum, which is 1, when only one species persists.

Fig. 2.

An emerging biphasic dependence of biodiversity on partitioning level

a. Overall framework of spatial partitioning. The number of local communities (N) changes while the total initial cell number and total volume keeps constant.

b. Simulations show that the partitioning level has no impact on the biodiversity of two populations with no interaction. δ1 = 0.1 and δ2 = 0.8 and all elements in the interaction matrix are zero. Data are presented as mean values +/− SD and the same applies to all the following panels and figures. Simulations were repeated 10 times (n=10) to calculate the mean and SD and the same with panel e-i.

c. Simulations of a pairwise one-directional negative interaction show that increased partitioning reduces the negative impact from A on B that leads to higher relative abundance of B and an increased biodiversity.

d. Simulations of a pairwise one-directional positive interaction show that increased partitioning reduces the positive impact from A on B and leads to decreased biodiversity.

e. Biphasic response emerges with both positive and negative interactions in 20-member communities with 100% connectedness. The panel titles indicate the fraction of negative interactions among all interactions. δ, β, and γ are sampled from uniform distributions of interval (0, 1.5), (0, 3), and (0, 0.8). The following panels share the same parameter setting while assuming 1:1 ratio of positive to negative interaction and varying the corresponding parameter and its specific value is indicated as panel titles. In panels e-i, grey dots represent the responses of a random interaction network, and the black dots are average response.

f. Increasing number of species increases the amplitude of biphasic response.

g. Increasing connectedness increases the amplitude of biphasic response.

h. Increasing interaction strengths increases dynamic range of the response. Interaction strength is used as a multiplier to scale the range of interaction parameters in panel e.

i. Biphasic dependence also emerges with intermediate ratio of positive and negative interaction strengths. The total interaction strength is held constant at 2. The ratio of maximum positive to maximum negative interaction strength ranges from 0.1 to 10.

Indeed, the same type of natural microbial community can reside in environments with different partitioning levels. Starting with two-member communities, when there is no interspecific interaction, spatial partitioning has no impact on final global community diversity and composition (Fig. 2b). However, increasing partitioning promotes the biodiversity of a pairwise community where one member suppresses the other, by shielding the suppressed strain from its suppressor (Fig. 2c, Extended Data Fig 1a). In contrast, increasing partitioning decreases the biodiversity a pairwise community when one member promotes the growth of the other. This is due to the partitioning of the dependent strain from its helper strain, which reduces its growth (Fig 2d, Extended Data Fig. 1b). This observation extends to communities with more populations that are dominated by either positive interactions or negative interactions (Extended Data Fig. 1c). Even for a pairwise community, the degree by which partitioning affects biodiversity also depends on the strength of interactions (Supplementary Fig. 2).

Extended Data Fig. 1

Intuition of the dependence of community diversity on partitioning level

a. Local community behaviors for one-directional negative interaction. Low partitioning allows the interaction between the two populations, leading to the collapse of the green population. High partitioning separates the two populations and enables the growth of the green population.

b. Local community behaviors for one-directional positive interaction. Low partitioning levels enable the growth of the green population, whereas the green population cannot grow well at high partitioning levels in absence of the other population.

c. The same principle applies to 12-member communities with all negative interactions (red) or all positive interactions (blue). Increasing partitioning increases the diversity of communities with negative interactions but decreases the diversity of communities with positive interactions. Data are represented as mean values +/− SD, with n=10, error bar represents standard deviation.

d. Steady state is not required to generate biphasic response. The simulation results of a fully connected 15-member community at different t values: 50 (dark grey), 100 (grey), and 200 (light grey). The community has 50% of negative interactions and 50% of positive interactions, with a maximum δ of 1.5, a maximum γ value of 1, and a maximum β of 5.

e. Biphasic dependence of diversity of local communities, which is quantified as the count of local communities that have unique combinations of members. The plot is generated with a 10-member community; each dot represents one randomization of the initial seeding. The solid trace represents the average number +/− SD (error bars) and n=10 and the same as panel f.

f. Biphasic dependence of diversity of local communities containing a population. The plot is generated with the same 10-member community in panel e, with one set of initial seeding. Each dot represents the diversity of local communities containing a population at a partitioning level. The solid trace represents the average number across 10 members.

g. More unique types of local communities lead to higher diversity. For communities with both positive and negative interactions, types of local communities correlate with final community diversity. The lighter the trace, the more members there are in a community.

When communities have both competition and cooperation, simulations reveal a biphasic dependence of biodiversity on the partitioning level for a large range of relative prevalence between cooperation and competition (Fig. 2e). When there is a balance of positive and negative interactions, biphasic dependence emerges with increasing number of species in the community, increasing overall interaction strength, and increasing interaction connectedness (Fig. 2f–i). Reaching the steady state is not a requirement for biphasic dependence to occur, though the biphasic dependence becomes more pronounced when approaching the steady state (Extended Data Fig. 1d).

The biphasic dependency can be explained from two aspects. One is the suboptimal biodiversity at both low or high partitioning levels, where the suppressed members or the dependent members have reduced growth or are unable to persist, respectively. The other aspect is that for each species, only a subset of initial local community compositions can enable its persistence (by excluding its competitors or providing its cooperators). Thus, increasing the number of unique local communities increases the chance for each species to persist in at least one local community. Since the highest count of unique local communities at t0 (Extended Data Fig. 1e) and the count of unique local community containing any species (Extended Data Fig. 1f) both peak at an intermediate partitioning level, the intermediate level can maintain the growth of the most species. Consistent with this notion, we find that indeed, the highest biodiversity overlaps with the highest diversity of local communities (Extended Data Fig. 1g). However, the biphasic dependency cannot be explained by local community biodiversity because it decreases with increasing partitioning levels regardless of interaction types (Extended Data Fig. 2).

Extended Data Fig. 2

Distribution of the biodiversity of local communities.

The patterns of the distribution of local patch biodiversity are similar regardless of the nature of interactions. Intermediate partitioning levels lead to bimodal distributions. The black lines indicate the biodiversity of the global community; the light purple lines show the average of local community biodiversity and error bars show the standard deviation of the local community diversity with N = # of partitions. The histogram shows the distribution of local community biodiversity corresponding to each portioning level. Normalized count is the count of local communities in each bin normalized by the total number of local communities in the partitioning level.

a. 20-species communities with all negative interactions. The interaction networks have 100% connectedness. δ is generated by a uniform distribution with minimum equals to 0 and maximum equals to 2. β is generated by a uniform distribution with minimum equals to 0 and maximum equals to 3. γ is generated by a uniform distribution with minimum equals to 0 and maximum equals to 0.8.

b. 20-species communities with 1:1 count of negative to positive interactions. All parameter settings are the same as panel a, except for the split between negative and positive interactions.

c. 20-species communities with all positive interactions. Intermediate partitioning levels lead to bimodal distributions. All parameter settings are the same as panel a, except for the split between negative and positive interactions.

Robustness of the biphasic dependence

So far, we have found that biphasic dependence arises when there is a balance of positive and negative interactions, including their magnitudes (Fig. 2). This behavior is widely applicable to a broad range of parameter settings (Supplementary Fig. 3). In nature, microbial communities often undergo successive mixing and partitioning. Our simulations indicate that, in the presence of intermittent mixing, the general impact of partitioning on biodiversity is similar to the case without mixing (Supplementary Fig. 4). That is, even when communities go through multiple cycles of partitioning, growth, and mixing, our conclusion holds.

We then examined the robustness of our conclusion with respect to the model formulation. We tested a modified gLV model formulation. Compared with our original formulation, the standard gLV model can generate unbounded growth of the community. To avoid this situation, we introduced a carrying capacity to cap each species’ growth. The results are consistent with results from our base model formulation (Supplementary Fig. 5).

To test the effects of temporal stochasticity, we used a stochastic differential equation model. With a moderate level of noise, the stochastic model generates qualitatively the same results as simulations without temporal stochasticity (Supplementary Fig. 6). To simulate the variation of size and the level of nutrient of local communities, we implemented a randomization of local community carrying capacity within each partitioning level. Even when the local community carrying capacity has a normal distribution of sigma equals to 50% of mean, the results are qualitatively similar to the base simulation results (Supplementary Fig. 7). Taken together, these additional simulations demonstrate the robustness of the general conclusion as revealed by the base model simulation (Fig 2).

Biphasic dependence with simple communities

We use microtiter plates, each with 6, 24, 96, 384, or 1536 isolated wells, to implement various levels of spatial partitioning in experiments (Fig. 3a). In each experiment, we calibrated the initial density of each community such that the average number of cells in each local environment is ~0.5 at the highest partitioning level. We then allocated the same total volume of the same mixture of microbial community into different wells in each plate. After culturing for 30 hours, we pooled all the wells in each plate to measure the global community structure.

Fig. 3.

Experimental demonstration of the predicted principle with simple communities

a. Spatial partitioning implemented with microtiter plates. The same starting cell culture was dispensed into wells of 6, 24, 96, 384, and 1536 plates. The total volume was constant across all plates and evenly distributed into all wells within the plate. After culturing, each plate was pooled to measure the overall community structure.

b-c. Schematic of the one-directional negative (strain 1 and 2) and positive (strain 1 and 3) pairwise interaction.

d. Increased biodiversity for pairwise negative interactions. The final community structures of are quantified by selection plating and CFU counting. Data are represented as mean values +/− SD and n=16. The same method applies to panel e. 10 nM of aTc and 100 mM of IPTG were used to control circuit functions, and the same with panel e.

e. Decreased biodiversity for pairwise positive interactions. The QS signal from strain 1 promotes the growth of the receiver strain 3 by inducing the production of CcdA, which is the antitoxin to CcdB. Data are represented as mean values +/− SD and n=16. The slight increase in at low partitioning level was likely due to background-level negative interactions when the two strains compete for nutrition and space.

f. Simulated results of experimental pairwise negative community response to partitioning. Both simulations have been run 10 times and the mean (black open circle and grey trace) and standard deviation (grey error bars) are shown. The negative interaction has δ = [0,0], γ = [0,0;0.7,0], and β = [0,0;0,0]. Data represented as mean values +/− SD with n=10 (the same applies to panel g)

g. Simulated results of experimental pairwise positive community response to partitioning. The positive interaction is simulated using δ = [0,1.1], γ = [0,0.15;0.15,0], and β = [0,0;0.8,0].

We first investigated two pairwise synthetic communities, which use quorum sensing (QS) to mediate one-directional positive interaction or one-directional negative interaction. To implement the interactions, we used an engineered strain 1 that produces 3OC6HSL, which is a QS signal (Extended Data Fig. 3a). 3OC6HSL induces the expression of CcdB in strain 2 that results in cell death, forming the negative interaction (Extended Data Fig. 3b). 3OC6HSL induces the expression of CcdA which reduces the toxicity of CcdB that strain 3 produces, forming the positive interaction (Extended Data Fig. 3c). We designed each strain to have a different antibiotic resistance profile to use selection plating to quantify community composition. The circuit functions were validated using monocultures (Extended Data Fig. 3d) and interactions on agar plates (Extended Data Fig. 3e). Consistent with model predictions, increasing spatial partitioning promoted the biodiversity of the negative-interaction pair (Fig. 3b, Extended Data Fig. 3f) but suppressed that of the positive-interaction pair (Fig. 3c, Extended Data Fig. 3g). Further details of these circuits are provided in Supplementary Table 1.

Extended Data Fig. 3

Characterizations of the synthetic strains

a.-c. Circuit diagrams of strains 1, 2, and 3

d. Monoculture response to inducers and QS signals. Growth of strain 1 was not strongly impacted by inducers and QS signals. Growth of strain 2 was inhibited by 3OC6HSL, as well as the induction of circuit by aTc and IPTG. The growth of strain 3 was inhibited by IPTG that induces CcdB but rescued by addition of 3OC6HSL. The experiments were done with 1000-fold dilution of overnight monocultures. The strains were cultured in M9 medium at 30°C in a plate reader. “-” indicates no inducers or QS signals were added. “a” indicates the addition of [aTc] of 10nM. “I” indicates the addition of [IPTG] of 1mM. “C6” indicates the addition of [3OC6HSL] of 10nM.

e. Response of strain 2 and 3 to supernatant of 1. The supernatant of strain 1 (introduced in the center) inhibited growth of strain 2 (initially spread on the entire plate) around the center. The rescue strain 3 by strain 1 was confirmed by the elevated growth at the center of the plate. The agar plates were made with LB medium, 1.5% agar. 1mM IPTG and 100nM of aTc were also added to induce circuit functions. The overnight culture of strain 1 was induced by 1mM of IPTG and 100nM of aTc to produce 3OC6HSL.

f. Community structure of the pair with negative interaction. Strain 2 had an increased relative abundance with increasing partitioning level. The community structures were measured by selective plating. Data are represented as mean values +/− SD and n=16.

g. Community structure of the pair with positive interaction. Strain 3 has a reduced relative abundance with increasing partitioning level. The slight reduction of biodiversity at the low partitioning was due to the slight decrease of relative abundance of strain 3 possibly driven by background competition between the two strains. Data are represented as mean values +/− SD and n=16.

Many studies have tested the assembly of a group of species by coculturing all assemblages of a certain number of species out of a larger set of species. This experimental setup provides an alternative implementation of spatial partitioning where the physical separation is controlled through controlled initial seeding (Extended Data Fig. 4a). Consistent with our simulation and experimental results, simple experimental communities dominated with negative interactions demonstrate increasing biodiversity with increasing spatial partitioning, whereas communities dominated by positive interactions show an opposite trend. This is true for all 6 types of pairwise interaction diversities (Extended Data Fig. 4b) and for larger communities (Extended Data Fig. 4c). Communities with both negative and positive interactions also reveal biphasic dependence of diversity on partitioning level (Extended Data Fig. 4d). Previous studies have demonstrated that partitioning decreases diversity for mutualistic pairs[11,34]; partitioning increases diversity for two-directional negative interactions[35]; partitioning increases[35] or decreases the diversity of a pair with both positive and negative interactions. Previously published results of multi-member microbial communities also follow the same general principle (Extended Data Fig. 5).

Extended Data Fig. 4

Spatial partitioning by controlled seeding

a. Partitioning based on controlled seeding. Number of species in local communities (group size) decreases with increasing partitioning. All combinations of the same group size were tested and pooled after growth to determine the richness of the pooled community -- the total number of populations. Due to the exponential nature of combination, controlled seeding is difficult to implement in an exhaustive manner for communities with large number of populations.

b. Simulations of all 6 major types of interactions in monocultures and cocultures. A check (cross) mark indicates that a population can (cannot) survive by itself. Overall, partitioning promotes coexistence for negative interaction-dominated pairs and impedes that for positive interaction-dominated pairs. The predator-prey interaction shows both trends, depending on whether the populations survive by themselves. Previously published data are consistent with these simulation results.

c. Spatial partitioning by controlled seeding in large communities with only negative or positive interactions. Increasing partitioning increases diversity for communities with only negative interactions and decreases for communities with only positive interactions. Each dot represents a randomly generated interaction network, and 10 networks are generated for negative interaction (red trace) and positive interaction (blue trace) networks. The open circles represent the mean across the 10 random interaction networks and error bars represent the standard deviations.

d. Robust biphasic dependence is observed with controlled seeding. Partitioning implemented by controlled seeding also generates robust biphasic dependence for communities with both negative and positive interactions. Simulations are done on 10 randomly generated interaction networks of 10 populations with 1:1 ratio of positive versus negative interactions (grey dots). Open circles represent the mean and error bars represent the standard deviation.

Extended Data Fig. 5

Controlled seeding as an alternative implementation of spatial partitioning.

Local group size refers to the number of populations that are seeded into a local community at t0. Reducing group size is effectively increasing partitioning level.

a. An 8-member community that is dominated by negative interactions reaches higher diversity with increasing. The previous study has collected eight soil bacterial species and found that the final community compositions are primarily driven by competitive exclusion. We have used the published data that are generated by co-culturing all possible combinations of a specific number of species to test our theory.

b. A 14-member auxotroph community reaches maximum diversity at an intermediate partitioning level, with a sharp drop of diversity from local group size of 2 to 1. This previous publication investigated 14 E. coli auxotrophs where none can grow as monocultures but some pairs grow collectively in cocultures. When all 14 auxotrophs are cocultured, only a few coexist whereas others are competed out due to competition.

c. The biodiversity of a 12-member synthetic human gut microbiome consortia that has both positive and negative interactions follow a biphasic dependence on partitioning. This previous study has cultured the single species, all pairwise assemblages, all single-species dropout communities, and all 12-member community. For each community, the abundance is measured using 16S rRNA gene sequencing.

Biphasic dependence with complex communities

We next examined the applicability of our insights (Fig. 2) to much more complex experimental communities. To this end, we generated a collection of plasmid-barcoded Keio strains[36] to enable quantification of the community dynamics by next-generation sequencing (Extended Data Fig. 6). We constructed a community of 47 Keio strains that are auxotrophs (Supplementary Table 2). The auxotrophic strains were selected based on their final density in the minimum medium (MOPS) versus rich medium (LB) presented in the original Keio strain publication and we obtained 47 strains that were successfully barcoded.

Extended Data Fig. 6

Construction of barcoded Keio strains and sequencing quantification.

a. The backbone of barcoded plasmids. b. Calibration experiment. Samples prepared at known fixed concentrations were generated and prepared through the entire workflow of NGS library preparation, sequencing and data analysis, as described in the methods. The image represents the layout of sample prepared at known ratios for sequencing where daker shading indicates higher relative strain concentrations in mixture. Samples were prepared based on 2-fold dilution between each group.

c. Identifying outliers. A box and whisker plot of 94 samples sequences at equal ratio was used to identify barcodes that amplified poorly or over-amplified as compared to other barcode sequences. The center of the box and whisker plot is 762. Outlier barcodes were defined as those giving sequencing counts with a distance greater than 1.5 times below the 1st quartile’s interquartile range (distance between the 1st and 3rd quartile) or 1.5 times above the 3rd quartile’s interquartile range. Eight barcodes were removed from future analysis.

d. Overall calibration results. Normalized relative abundance of sequencing counts obtain for each barcode plotted versus the expected sample ratio shows a good correlation between expected and actual barcode abundances. Each data point represents a single barcoded Keio strain.

e. Correlation between NGS measurements of technical replicate experiments. Sequencing technical replicate of three different samples sequences independently.

f. Correlation between replicate NGS measurements on the same biological samples. Sample preparation and sequencing replicate for three sample starting from the same template DNA that were processed independently through sample preparation and sequencing protocols.

These auxotrophs compete for the nutrient components other than the amino acids they provide to each other. Their positive interactions only emerge when these amino acids are absent. Therefore, increasing the concentrations of these amino acids attenuates the positive interactions, elevating the relative contribution of negative interactions (Fig 4a). This modulation of the relative magnitude of positive and negative interactions by adjusting amino acid concentrations has been demonstrated in previous auxotrophic communities[14,15]. In our system, the positive interaction in the absence of exogenously added amino acids was verified by collective survival of the community (Extended Data Fig. 7a). We also measured the distribution of OD in each plate to verify that sufficient level of partitioning was achieved (initial density of ~0.5 cells per well in 1536 plate), and an increase of OD was observed with increasing amino acid concentration (Extended Data Fig. 7b). Sequencing results revealed an increasing diversity at high casamino acid concentration and a decreasing diversity at low concentration, and a biphasic dependence at an intermediate concentration (Fig. 4b, c).

Fig. 4.

Biphasic dependence observed with complex communities.

a. Schematic of a 47-member auxotroph community. The interactions in the community are modulated by the concentration of casamino acid. When the casamino acid concentration is low, the primary interaction is cooperation; when the concentration is high, the primary interaction is competition.

b. Partitioning experiment using an auxotrophic complex community. A selection of 47 auxotrophic strains were identified from the Keio collection and barcoded for NGS quantification. Complex communities were developed and grown at five levels of partitioning and varying amino acid concentration for 30 hours. Decreasing amino acid concentration shifts communities from negative interaction dominated to positive interaction dominated while transitioning through a state where both negative and positive interactions are strong. The panel presents representative relative abundance of each strain at the end of the experiment, as measured by NGS, for a single biological replicate.

c. Validation of the principle by a 47-member auxotroph community. Simpson index of communities demonstrates the effects of partitioning and amino acid concentration on population biodiversity. Each dark open circle represents the average among the 2 sequencing runs of 3 biological replicates (grey dots). Error bars represent standard deviation of the 6 replicates (N=6) and error bar centers represent the means. One-sided t test was used to test against a null hypothesis of a 0 slope. Degrees of freedom of the three-piece models are 10, 16, and 10 and all other panels are 28. The p values are: 0.1%: 0.1078, 0.02%: 0.0199, 0.005%: 0.0060 (left), 0.4980 (middle), 0.0168 (right); 0.001%: 0.0132, 0.0002%: 0.2371. The asterisks indicate that the trends of data points enclosed by the corresponding brackets have p-values less than 0.05.

Extended Data Fig. 7

Characterizations of the Keio auxotrophic (a & b) and non-auxotrophic (c) communities.

a. Community growth response to amino acid concentrations and initial density. Collective growth demonstrates the presence of positive interactions among 47 Keio auxotroph strains. With no casamino acids added (the darkest trace) and 0.0002% of casamino acid, the auxotroph Keio community shows no significant growth at 48 hours in M9 medium at 30°C despite having an initial density above 0. However, with higher initial cell densities, the community began to grow, which is a typical behavior of cooperative communities. The data are represented as mean with error bars representing standard deviations across 7 replicate wells (N = 7). From the darkest to lightest, the lines casamino acid concentrations of 0%, 0.0002%, 0.001%, 0.005%, 0.02%, and 0.1%.

b. The distribution of OD600 of each well at different casamino acid concentrations and partitioning levels. As expected, there is systematic increase in OD600 with increasing casamino acid concentrations (labeled at each row). Based on the initial cell density and Poisson distribution, 96 and 384-well plates are almost always are seeded with at least one cell in each well. However, at low casamino concentrations (0.0002% and 0.001%), bimodal distributions occur at 96 and 384-well plates indicating that the initial community compositions and their interactions play a crucial role in the growth of a strain. At the 1536 partitioning level, high casamino acid concentrations (0.1% and 0.02%) also show bi-modal distributions, which is expected due to a probability of ~30% wells being empty at a λ = 1.2.

c. Increased diversity with a 94-member community dominated by negative interactions. The strains were randomly selected from the Keio collection with auxotrophic strains excluded. Since all strains share the same genetic background, strong negative interactions are expected due to competition for nutrients and space.

We also constructed one community consisting of 94 non-auxotrophic strains (Supplementary Table 3), where none of the strains overlap with the 47 auxotrophic strains. The culture medium was supplemented with 0.1% casamino acid, which is a typical casamino acid concentration for M9 medium. The 0.1% [CA] also provides a direct comparison between the 94-nonauxotrophic community and the first condition of the auxotrophic community. Since Keio strains are all derived from Escherichia coli K-12 strain BW25113 with single gene deletions, the strains share similar genetic backgrounds and metabolic pathways, leading to primarily competition for nutrition in the community. Indeed, we have found that increasing partitioning level overall promotes the diversity of the non-auxotroph community (Extended Data Fig. 7c).

A previously-published dataset on natural groundwater bacterial communities also exhibits a biphasic dependence of biodiversity on partitioning[37]. The community underwent serial dilutions in microtiter plates where one well represents a local community and there are 96 local communities for each dilution level (Extended Data Fig. 8a). The samples were then cultured in the presence and absence of oxygen. There are three variables that change across serial dilution: partitioning level, average initial density, and down-sampling of the natural community. By accounting for both changes in initial density (Extended Data Fig. 8b) and down-sampling (Extended Data Fig. 8c), we find that biodiversity indeed follows a biphasic dependence on partitioning level (Extended Data Fig. 8d). The biodiversity curves for anaerobic and aerobic conditions suggest that aerobic condition creates an environment that favors competition, which is consistent with the conclusions in the original publication (Extended Data Fig. 8e). Water content is another implementation of spatial partitioning that shows biphasic dependency of biodiversity for soil microbiomes at microscale level[21].

Extended Data Fig. 8

Deducing the effect of spatial partitioning on a groundwater community

a. The experimental setup of the groundwater study. Partitioning increases with higher dilution due to the decrease of local community sizes and reduced number of interactions. Serial dilution introduces two additional variables: the decrease of initial cell densities and subsampling of the initial community that occurs with increasing levels of dilution.

b. Negligible impact of initial cell density on final community diversity. 10 randomly generated 1000-member community were simulated with increasing dilution rate, which is equivalent to decreasing density of the initial community (grey dots). Data are represented as mean values +/− SD with N=10.

c. Estimated global initial community structure. The population abundance distribution loosely follows power law, which approximates the distribution pooled final communities of different dilution levels, especially at 10X and 102X. The estimated abundance distribution of the original sample is based on power-law distribution and an estimate of 5000 total OTUs, which is feasible compared with 399 OTUs that are sequenced in the final communities that account for both 10X dilution and OTUs that were not able to grow.

d. Comparison between the measured richness at each dilution level and the estimated richness of initial communities. The baseline is the mean sampled from the estimated community composition of the original sample by simulating serial dilution. The error bars of the baseline are the standard deviation across 10 simulated samplings (n=10). The # of OTU present was calculated from the published data. Dilution of 10 in anaerobic condition has no OTU present and it is thus not shown.

e. Biphasic dependence of a groundwater community. The community was analyzed in both anaerobic (NO3) and aerobic (O2) culture conditions. Down sampling due to dilution was accounted for to estimate the number of OTUs sampled at t0. Y-axes indicate the percentage of OTUs that are present at t (after culturing) out of the number of sampled OTUs at t0 (see SI for detailed method). Data are represented as mean values +/− SD with N=10.

Responses of archetypes and multilevel partitioning

Within each community, member species can be categorized in 8 archetypes based on fitness level and the interactions they receive (Extended Data Fig. 9a). Simulation of a 160-species community (20 species of each archetype) with a random interaction network shows how each archetype responds to various partitioning levels (Extended Data Fig. 9b). There are four types of responses (Extended Data Fig. 9c): flat, positive, negative, and biphasic. Flat curve corresponds to species receiving no interactions. Positive and negative curves correspond to species receiving negative and positive interactions, respectively, whereas biphasic curve corresponds to species receiving both types of interactions. Note the parallel between the response of individual archetypes and how the biodiversity of communities of one archetype changes across partitioning levels.

Extended Data Fig. 9

Response of members of different architypes to spatial partitioning

a. Eight archetypes of species based on fitness, competition, and cooperation, using three metrics: fitness (1 – δ), positive interaction it receives (β), and negative interaction it receives (γ). The color of each archetype is determined by the level of fitness, the strength of positive interaction received, and the strength of positive interaction received, which correspond to the level of blue, red, and green. Each archetype is named to reflect their characteristics.

b. A 160-member community comprised of all eight archetypes. The same simulation framework described in Figure 2 was used to simulate how each archetype responds to partitioning. 20 populations were defined for each archetype. The interaction matrix was randomized while constraining the strength of interactions each population receives according to its archetype. The size of each dot represents the final population density.

c. Relative abundance of the eight archetypes. Each trace is the average relative abundance of all 20 species that belong to the same archetype. Higher partitioning level enriches species receiving negative interactions (traces with colors that have red components) whereas lower partitioning level enriches species receiving positive interactions (traces with colors that have red components). Instead of white, we plotted the relative abundance of “Drama Queens” in grey. Data are represented as mean values +/− SD with N=20.

d. Mixed partitioning level provides a robust way to maintain community biodiversity. From the top left to bottom right panel, proportion of negative interaction change from 0 to 1.0 for a 40-species community. Even though the biodiversity response changes from negative monotonic dependence to biphasic dependence, and to positive monotonic dependence, mixed partitioning, which in this set of simulation is even volume mix across partitioning levels, creates a robust way to maintain biodiversity of the community (bar on the right). Data are represented as mean values +/− SD with N=10.

The diversity of responses by different archetypes suggests that a single partitioning level cannot accommodate all species. Although there exists a level of partitioning that optimizes the biodiversity of microbial communities, this optimal biodiversity may not ensure the persistence of all species (Fig. 5a). Furthermore, without a priori characterization of the interactions within the community, the optimal partitioning level is unknown. To this end, we reasoned that a mixed partitioning level could ensure the maximum possibility for all species of any archetype to persist. Simulation results show that when more partitioning level are included, more species will be able to persist (Fig. 5b) and mixed partitioning level provides a robust strategy to maintain high community diversity (Extended Data Fig. 9d).

Fig. 5.

Multilevel partitioning is a robust strategy to maintain community biodiversity

a. Sub-optimal maintenance of diversity with a homogeneous partitioning level. We simulated the impact of spatial partitioning on a random 20-member community with both positive and negative interactions and non-archetypical species. The populations persist at different partitioning levels. Although an intermediate partitioning level maintains the highest diversity, it does not maintain all 20 species.

b. Robust maintenance of diversity through mixed level of partitioning. With more partitioning levels included, the persistence of more species can be maintained. When all 5 partitioning levels are included, all 20 species can be maintained. Instead of using any single partitioning level, a mixed level of partitioning allows most populations to persist.

c. The cross-section of the cellulose sponge and experimental procedure. The pore size has a wide distribution that naturally creates multi-level partitioning for the global community. The sponge is then fit into an autoclavable culture tube. After autoclaving, the same initial cell culture as used in Fig. 4 is added to the tube and saturates the sponge.

d. Multilevel partitioning to maintain biodiversity. Without prior knowledge of the interactions within a community, it is challenging to choose a single partitioning level that best maintain the community diversity. Multilevel partitioning can promote the chance for any archetype to thrive in the partitioning level it prefers. We cultured the Keio collection auxotroph community either as a homogeneous liquid culture (−) or with a sponge (+) that has varying pore sizes that provide a multilevel partitioning. Three amino acid concentrations (x axis labels) were tested, and all have shown an increase in community diversity with multilevel partitioning compared with liquid culture (no partitioning).

We tested this notion using cellulose sponges to create multi-level partitioning due to the wide distribution of pore sizes inside a sponge[38] (Fig. 5c). We used the same microbial community shown in Fig. 4 to test if submerging a sponge in the liquid culture with low initial cell density will increase the biodiversity of the final community compared with a shaken liquid culture. In addition, we used three casamino acid concentrations to test how well the sponge can maintain a higher biodiversity with different overall strengths of positive and negative interactions. Indeed, our results show that the global diversity of communities cultured in sponge was overall significantly higher than those grown in liquid culture (Fig. 5d).

In nature, the physical architectures of microbial habitats indeed impose multiple levels of partitioning through interconnected local communities, varying sizes of local environments, and self-assembled local patches. Our results demonstrate a robust mechanism that can explain the high biodiversity observed in natural microbial communities. Following this mechanism, multi-level partitioning also facilitates the emergence of spatial heterogeneity and niche differentiation, which are key mechanisms that maintain biodiversity[39]. Since spatial partitioning is important to maintaining the biodiversity of natural microbial communities, its disruption in laboratory cultures most likely reduces biodiversity of the community of interest. Gut and soil microbiome diversity may be better maintained with culturing methods implementing multi-level spatial partitioning.

Discussion

Our study reveals a simple principle that dictates how spatial partitioning modulates microbial community structure through global modulation of interactions. Although spatial partitioning modulates the growth of populations receiving negative versus positive interactions in opposite directions, it is robust that complex communities with both negative and positive interactions reach the highest biodiversity at an intermediate partitioning level. Previous studies that use precisely controlled assembly primarily focus on characterizing pairwise interactions, which limits the scalability of the experiments[23]. In contrast, we leveraged seeding stochasticity and community-level interaction characteristics to demonstrate the robustness and scalability of this principle. Further, this community-level principle is the result of the collective response of single species in the community.

For experimentalists, our study suggests the importance of the explicit consideration of physical arrangement of communities in experimental designs. For example, serial dilutions not only modulate the initial density of communities but also increase the spatial partitioning level of the community. Our study also suggests that to maintain the highest biodiversity of a natural community, beyond the design of chemical environment, the design of the physical environment is also crucial. For natural microbial communities that arise from complex spatial partitioning environments, spatial design of lab cultures can be especially important for maintaining their structure and biodiversity.

Engineering of habitat by spatial partitioning can be an effective strategy to modulate and control microbial community structures. Beyond using microtiter plates, other engineering methods can also be used to impose spatial partitioning that has not been fully tested yet, such as encapsulation[40] and inkjet printing[41]. The modulation of community structure can happen at three levels. First the relative abundance of a single population can be modulated based on archetype. Second, the proportions of interactions can be modulated through increasing or decreasing partitioning. Third, the biodiversity can be modulated or maintained through the modulation of spatial partitioning level or a mixed level of partitioning. Beyond being a general and robust modulator, spatial partitioning, as a physical factor, is orthogonal to chemical and biological factors and can be used in parallel.

The role of space in maintaining biodiversity is a central question beyond microbial ecology, as exemplified by classic Island Biogeography Theory[42], debate regarding the optimal design of biodiversity reserves[43,44], and the subsequent flowering of metacommunity ecology[24,45,46]. Many specific aspects of partitioning, such as nestedness[47], stochastic extinction[48], and migration rate[49] have been investigated. However, most studies largely overlook interactions, in particular cooperation. Our general principle of spatial partitioning accounts for both types of interactions, addressing a critical gap in our understanding of spatial mechanisms for the maintenance and promotion of biological diversity. This helps to clarify the role of spatial mechanisms alongside temporal mechanisms of biodiversity maintenance and spatiotemporal mechanisms, such as the intermediate disturbance hypothesis[50]. Finally, the insights presented here and offers a fresh perspective for interpreting, screening, and controlling microbial community structures and the relative abundance of individual population.

Methods

Data availability

Experimental data generated for this manuscript are provided as source data files, and deposited to GitHub at: https://github.com/youlab/partitioning_NCB2021.

Code availability

The simulation and data-analysis codes used in this the study are deposited to GitHub at https://github.com/youlab/partitioning_NCB2021.

QS-based and synthetic E. coli strains

The strains are constructed from previously published strains and their antibiotic resistance are verified using antibiotic plates. Strain 1 is previously published as QS-CAT strain in[51]. Strain 2 is identical to the circuit of prey strain except that the CcdBs is replaced with the wild-type CcdB gene[52]. Strain 3 is previously published in[11] as M2 in the synthetic mutualistic pair.

Growth conditions

M9 medium contained standard M9 salt (6.8 g/L Na2HPO4, 3 g/L KH2PO4, 1 g/L NH4Cl, and 0.5/L NaCl), 0.5% (w/v) glucose, 1 mM MgSO4, 0.1 mM CaCl2, 1 μg/ml thiamine, and 0.1% (w/v) casamino acid. It was buffered with 0.1M of MOPS and adjusted to a pH of 7.0 with NaOH. The antibiotic concentrations used in this study to maintain the plasmids during overnight culture are: 100 μg/ml chloramphenicol, 50 μg/ml kanamycin, 100 μg/ml ampicillin, and 50 μg/ml spectinomycin. The agar plates were made with 1.5% agar and standard Luria-Bertani (LB; Miller) broth, containing the same antibiotic concentrations.

We grew the strains with M9 medium at 30°C for 8~12 hours from glycerol stocks to create overnight cultures. We then calibrate the OD of the overnight culture to 0.1 and mix each strain by equal ratio followed by 2×106fold dilution. This calibration and dilution procedure ensures a starting total cell concentration of 50~200 cells/ml. Each well contains 2.56 ml, 640 μl, 160 μl, 40 μl, and 10 μl for 6, 24, 96, 384, and 1536 well plates, respectively. QS-based synthetic systems were cultured using M9 medium supplemented with aTc and IPTG but no antibiotics. The Keio strains were cultured with 100 μg/ml chloramphenicol to maintain the barcoded plasmids. The Keio auxotroph strains were cultured in the same M9 medium except for the concentration of added casamino acid. The strains did not lose antibiotic resistances over the experiments. The plates were sealed using both AeraSeal (Excel Scientific) and Breathe‐Easy membrane (Diversified Biotech) and are covered by plastic lids. The microtiter plates were shaken at 225 rpm at 30 °C. for 24 hours for the QS-based strains. We cultured the regular Keio strains for 30 hours and auxotrophic Keio strains for 48 hours to reach sufficient cell density.

Quantifying community structure of QS-based communities

After culturing, we pooled all wells within each microtiter plate and measured the OD of each pooled community. Each sample was diluted and plated on multiple types of antibiotic plates with 4 replicate plates. The antibiotic plates were cultured at 37°C overnight before CFU counting. Since each strain has a unique antibiotic resistance profile, all strains can be measured using the correct type of antibiotic plates. For each pairwise community, the replicate plating gave 16 pairs (4×4) of measurements, which were used to determine community structure and biodiversity index.

Construction of barcoded Keio strains

The 96 non-auxotrophic strains were selected randomly from the Keio collection. We identified autotrophs by screen all Keio strains with OD in LB greater than 0.45 and OD in MOPS below 0.05. Among the total of 76 auxotroph strains that were screened, we randomly picked 47 to construct the auxotrophic community.

A library of barcoded plasmids was generated using Gibson Assembly cloning. Briefly, a plasmid vector backbone was generated by linearizing a plasmid vector (p15A-GFP-Carbenicillin, Extended Data Fig 6a) using PCR amplification (P1 and 2, Q5 MasterMix). Using lower amount of template DNA (as low as 0.1 ng of DNA) reduced the background of unmodified vectors downstream. The vectors were gel purified and underwent an overnight restriction enzyme digested using NOTI and PVUI, based on manufacturer’s instructions, to generate end that overlap two synthesized DNA fragments for Gibson Assembly.

Sequences of the synthesized DNA fragments are provided below. Each part contains a 15–20 base pair overlap with the plasmid vector (yellow) and with each other (red) for Gibson Assembly, one of the Illumina® adapter sequences in green (part 1 contains the i5 adapter and part 2 contains the i7 adapter sequence) and 18 random base pairs making up the barcode regions.

Part 1:

5’ -GCCTCAGGGCCCGATAGTACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNNNNNNNNNNNNNCCTCAGGGTCACTAGG -3’

Part 2:

5’ -CTCAGGGTCACTAGGNNNNNNNNNNNNNNNNNNAGATCGGAAGAGCACACGTCTGAACTCCAGTCACGTGGGCCGCTTAATTAATTAATC -3’

Following Gibson Assembly cloning, the barcode DNA fragment was assembled into the vector backbone in the following structure:

5’…GCCTCAGGGCCCGATAGTACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNNNNNNNNNNNNNCCTCAGGGTCACTAGGGNNNNNNNNNNNNNNNNNNAGATCGGAAGAGCACACGTCTGAACTCCAGTCACGTGGGCCGCTTAATTAATTAATC…3’

A library of barcoded plasmids was generated by transforming the Gibson Assembly reaction into electrically competent cells, allowing an overnight outgrowth of the transformed cells and harvesting the plasmid library by plasmid midiprep kit (Qiagen), based on the manufacturer’s instructions. Plasmids were then be transformed into target strains by chemical transformation[53]. Barcoded strains were validated by colony PCR (Supplementary Table 4, Primers 3 and 4, Q5 MasterMix). Following enzymatic clean-up, PCR products were sent for Sanger Sequencing to validate the barcode sequence of each strain.

A collection of strains from the Keio collection were uniquely barcoded, including more than 94 non-auxotroph strains and 47 auxotrophs. All barcoded strains were validated by sanger sequencing to contain a single unique barcode sequence. Sanger sequencing validation of all barcoded strains allowed for simplification of downstream data processing as all barcode sequences and associated strain identifies are known. The sequencing also allowed us to remove strains carrying duplicate or multiple barcodes as a result of double transformation.

NGS library preparation

Plasmid DNA was extracted from experimental samples by miniprep kit (Zymo Research), using the manufacturer’s instructions, or by boiling lysis of bacteria (95°C for 10 m in nuclease-free water). Extracted DNA was stored at −20°C until use in downstream NGS library preparation.

Libraries for NGS sequencing were prepared using a two-step PCR protocol using NEBNext® Ultra™ II Q5® Master Mix. The first PCR step is composed of 2 PCR cycles and is used to extract barcode sequences from the vector backbones using primers that bind the conserved Illumina adapter sequences flanking either end of unique barcode sequences (Green in sequences above). In addition, unique molecular identifiers (UMIs) and dual sample indexes are introduced during this PCR protocol. Primer sequences are provided in Supplementary Table 4, where the 10-base pair UMI sequence is represented as Y’s and the 8-base pair index sequences are represented as X’s in red. UMIs are used to uniquely label each template DNA molecular in the initial sample and allow for downstream correction for sequencing errors introduce by PCR amplification.

PCR 1 conditions:

Initial denature – 30 sec– 98°C

2x PCR cycles

10 sec – 98°C

30 sec – 67°C

20 sec – 72°C

Final extension – 5 min 72 C

Infinite Hold 4°C

Following PCR 1, PCR clean-up and size selection were achieved using SPRIselect magnetic beads (Beckman Coulter). Two-step size selection was used (0.95X → supernatant saved →0.8X → DNA eluted from bead) based on the manufacturer’s instructions. The cleaned-up PCR products are then pooled together and amplified using a second round of PCR using the conditions described below.

PCR 2 conditions:

Initial denature - 30s– 98°C

20x PCR cycles

10 sec – 98°C

30 sec – 66°C

20 sec – 72°C

Final extension – 5 min 72 C

Infinite Hold 4°C

The product from PCR 2 was run on a 2% agarose gel for PCR clean-up and final validation of the library size. DNA was extracted and purified from the gel using the Zymoclean Gel DNA Recovery Kit, based on the manufacturer’s instructions. The cleaned PCR product is the final sequencing amplicon which is compatible for sequencing on standard Illumina sequencing platforms. The DNA libraries were sequenced using 151 base pair, paired-end reads either through a sequencing facility on an Illumina MiSeq or in house using an Illumina MiniSeq. The data collection was performed in either MiSeq Software Suite or MiniSeq Software Suite. For in house sequencing, DNA libraries were denatured, diluted and mixed with a Phi-X spike-in of 30% based on standard Illumina protocols for library preparation of 16S Library on the Illumina Miniseq Platform.

Final NGS Library:

AATGATACGGCGACCACCGAGATCTACACXXXXXXXXYYYYYYYYYYACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNNNNNNNNNNNNNCCTCAGGGTCACTAGGGNNNNNNNNNNNNNNNNNNGATCGGAAGAGCACACGTCTGAACTCCAGTCACAXXXXXXXXATCTCGTATGCCGTCTTCTGCTTG

Where N’s represent barcode sequences, X’s represent sample indexes and Y’s represent UMIs.

Data analysis pipeline

Reads from NGS sequencing were analyzed using tools available on the open source, web-based platform, Galaxy (Galaxy version 20.05). Data analysis was simplified due to the fact that all sequencing reads are the same size, all barcode sequences are known a priori and forward and reverse reads are fully overlapping. Sequencing reads take the following format:

Read 1:

NNNNNNNNNNNNNNNNNNCCTCAGGGTCACTAGGNNNNNNNNNNNNNNNNNNAGATCGGAAGAGCACACGTCTGAACTCCAGTCACXXXXXXXXATCGCGGATGCCGGCTTATGGTTGGGAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGG

Read 2:

NNNNNNNNNNNNNNNNNNCCTAGTGACCCTGAGGNNNNNNNNNNNNNNNNNNAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTTYYYYYYYYYYXXXXXXXXTTTAGAAATACGTTGTACCCCTAACAATAAAAAAAAAAAAAGGGGGG

Where N’s represent barcode sequences, X’s represent sample indexes and Y’s represent UMIs.

Following quality control using FastQC software (Galaxy Version .72+galaxy1)[54], pooled paired end reads were demultipelexed into different experimental conditions (Barcode Splitter Galaxy Version 1.0.1[55], Trim Sequences Galaxy Version 1.0.2+galaxy0[55], and FASTQ joiner Galaxy Version 2.0.1.1+galaxy0)[56].

Next, paired end reads were merged using the software BBMerge. A high merge rate of 95% or above was achieved for all reads since forward and reverse reads are completely overlapping. Correction for PCR amplification error was performed on the Galaxy Platform using the following selection of tools: Sort Collection (Galaxy Version 1.0.0), FASTQ joiner (Galaxy Version .0.1.1+galaxy0)[56], UMI-Tools extract and deduplicate (Galaxy Version 0.5.5.1)[57], Bowtie2 (Galaxy Version .3.4.3+galaxy0)[58,59], and Samtools fastx (Galaxy Version 1.9+galaxy1)[60].

Finally, strain-specific barcodes were counting in each sample using a custom python script. Up to two mismatched base pairs were allow per barcode during barcode counting. This was selected as a reasonable number of mismatched bases that may result from PCR amplification error using Q5 high fidelity master mix and NGS sequencing error and considering that the minimum hamming distance between any two barcode sequences was 17 base pairs.

Calibration and replicates

Calibration samples were generated with barcoded strains prepared at known concentrations. As such, cultures were grown to a similar, low OD (0.3–0.4) and mixed at three different known ratios as well as one sample where all strains were mixed at an equal ratio (Extended Data Fig 6b). Samples were prepared through the experimental workflow, NGS sequencing and data analysis. The resulting sequencing reads correlated with the expected relative sample concentrations (Extended Data Fig 6d) demonstrating the reliability of the barcode sequencing approach. Normalizing the number of sequencing reads in each of the mixed ratio samples by the number of reads counts obtain in the equal ratio sample for each respective barcode improved the correlation between known and expected barcode counts. This normalization likely corrects for differences in actual strain abundance at the time of sample preparation, (i.e. variations in ODs at sampling time and differences in OD values and actual cell numbers for samples with similar OD may vary depending on cell size, etc.). There was one outlier data point in the ‘ratio 2’ sample. This strain did not appear as an outlier in other samples and was present in technical replicate sequencing of the same biological sample. As such, this error was most likely a result of human error in sample and was not associated to a problem with PCR amplification of the species barcode sequence.

Keio collection experiments

The characterization of Keio auxotrophic community (Extended Data Fig 7a–b) was generated with an array of initial densities and varying the concentrations of casamino acids.

The overnight monocultures of the 47 auxotrophic Keio strains were mixed by equal volume and calibrated to a cell density of 100 cells/ml. The medium is the same as the M9 medium described in section IV.2, except we adjusted the concentration of amino acid to 0%, 0.0002%, 0.001%, 0.005%, 0.02%, and 0.1%. The 3 replicates were generated in two separate days. The distribution of OD of the end points are shown in Extended Data Fig 7b.

The overnight monocultures of the 94 regular Keio strains were mixed by equal volume and calibrated to a cell density of 120 cells/ml, which is verified using CFU counting. M9 medium was used (see section IV.2). 6, 24, 96, 384, and 1536 plates were used to create 5 levels of partitioning. All microtiter plates were shaken at 225 r.p.m. at 30°C for 30 hours. The OD of the pooled samples were 0.375, 0.320, 0.220, 0.294, 0.179 for 6, 24, 96, 384, 1536 well plates. The measurements of the biodiversity are shown in Extended Data Fig 7c.

Sponge experiment

We used kitchen sponge from Scotch Brite and cut the sponge into 3 pieces with dimensions of 1.2 cm by 1.2 cm by 4cm. The sponges were sterilized in 15ml falcon tubes by autoclaving. We checked the sterilization of the sponges by adding LB growth medium to the sponge and shaken at 225rpm at 30 °C for 48 hours in a 15ml falcon tube. No cell growth was detected based on OD measurement.

We mixed the 47 Keio auxotroph strains by equal volume in M9-glucose medium and diluted the mixture to a cell density of around 300/ml, which is verified by CFU counting. M9-glucose medium with no casamino acid contains standard M9 salt (6.8 g/L Na2HPO4, 3 g/L KH2PO4, 1 g/L NH4Cl, and 0.5/L NaCl), 0.5% (w/v) glucose, 1 mM MgSO4, 0.1 mM CaCl2, 1 μg/ml thiamine. The M9-glucose medium is buffered with 0.1M of MOPS and adjusted to a pH of 7.0 with NaOH. We then created three variations of the M9 medium by adding 0.001%, 0.005%, and 0.02% casamino acid. 7ml of cell culture was added to 15ml falcon tubes with or without sponge. The cell cultures were shaken for 20 hours at 30°C to enable the cells to reach a higher cell density, then the cell cultures were cultured at 30°C for another 14 hours. We first applied physical pressure to the sponge multiple times to mix the cell culture and then we collected the supernatant for downstream processing.