Open access
Editor's Pick
Research Article
13 October 2020

Higher-Order Interactions Dampen Pairwise Competition in the Zebrafish Gut Microbiome


The microbial communities resident in animal intestines are composed of multiple species that together play important roles in host development, health, and disease. Due to the complexity of these communities and the difficulty of characterizing them in situ, the determinants of microbial composition remain largely unknown. Further, it is unclear for many multispecies consortia whether their species-level makeup can be predicted based on an understanding of pairwise species interactions or whether higher-order interactions are needed to explain emergent compositions. To address this, we examine commensal intestinal microbes in larval zebrafish, initially raised germfree, to allow the introduction of controlled combinations of bacterial species. Using a dissection and plating assay, we demonstrate the construction of communities of one to five bacterial species and show that the outcomes from the two-species competitions fail to predict species abundances in more complex communities. With multiple species present, interbacterial interactions become weaker, suggesting that higher-order interactions in the vertebrate gut stabilize complex communities.
IMPORTANCE Understanding the rules governing the composition of the diverse microbial communities that reside in the vertebrate gut environment will enhance our ability to manipulate such communities for therapeutic ends. Synthetic microbial communities, assembled from specific combinations of microbial species in germfree animals, allow investigation of the fundamental question of whether multispecies community composition can be predicted solely based on the combined effects of interactions between pairs of species. If so, such predictability would enable the construction of communities with desired species from the bottom up. If not, the apparent higher-order interactions imply that emergent community-level characteristics are crucial. Our findings using up to five coexisting native bacterial species in larval zebrafish, a model vertebrate, provide experimental evidence for higher-order interactions and, moreover, show that these interactions promote the coexistence of microbial species in the gut.


Intestinal microbes exist in complex and heterogeneous communities of interacting, taxonomically diverse species. The composition of these communities varies across individuals and is crucial to the health of the host, having been shown in humans and other animals to be correlated with dietary fat uptake (1, 2), organ development (3, 4), immune regulation (510), and a wide range of diseases (1120).
Despite the importance of intestinal communities, the determinants of their composition remain largely unknown. A growing number of studies map the effects of external perturbations, such as antibiotic drugs (21, 22) and dietary fiber (23) and fat (24, 25), on the relative abundance of gut microbial species. Intrinsic intermicrobial interactions, however, are especially challenging to measure and are important not only for shaping community composition in the absence of perturbations but also for propagating species-specific perturbations to the rest of the intestinal ecosystem.
The considerable majority of studies of the gut microbiota have been performed on naturally assembled microbiomes by sequencing DNA extracted from fecal samples, an approach that provides information about the microbial species and genes present in the gut but that imposes several limitations on the inference of interspecies interactions. The high diversity of natural intestinal communities and, therefore, the low abundance of any given species among the multitude of its fellow residents implies that stochastic fluctuations in each species’ abundance will be large, easily masking true biological interactions. The accuracy of inference is considerably worse if only relative, rather than absolute, abundance data are available (2629), as is typically the case in sequencing-based studies. Finally, we note that fecal sampling assesses only the microbes that have exited the host, which may not be representative of the intestinal community (30).
An alternative approach to using DNA sequencing and naturally assembled host-microbiota systems is to build such systems from the bottom up using model organisms. This is accomplished by using techniques for generating initially germfree animals, and well-defined sets of small numbers of microbial species, and then measuring the populations of these species resident in the intestine. Recent work along these lines has been performed using the nematode Caenorhabditis elegans (31) and the fruit fly Drosophila melanogaster (32, 33). However, as described further below, these studies imply different principles are at play in the different systems. Moreover, it is unclear whether conclusions from either model platform translate to a vertebrate gut, which has both greater anatomical complexity and more specific microbial selection (34). To address this, we measure bacterial interactions in larval zebrafish (Fig. 1A), a model vertebrate organism amenable to gnotobiotic techniques (3538), which enabled, in earlier works that investigated pairs of bacterial species, the discovery of specific interbacterial competition mechanisms related to intestinal transport (39, 40). The experiments described here involve several hundred fish, each with 1 to 5 resident bacterial species, enabling robust inference of interspecies interactions.
FIG 1 Five chosen commensal species are robust colonizers of the larval zebrafish intestine. (A) A 7-dpf larval zebrafish, with a dotted line outlining the intestine. Scale bar, 500 μm. (B) Chromogenic agar plate showing colonies of all five candidate species (A. calcoaceticus, milky opaque; Aeromonas sp. strain ZOR0001, reddish purple; Enterobacter sp. strain ZOR0014, blue; Plesiomonas sp. strain ZOR0011, dark purple; P. mendocina, colorless translucent). (C) The abundance per zebrafish gut of each of the five bacterial species when colonized in monoassociation with the host, assessed as number of CFU from plated gut contents. Each circular data point is a CFU value from an individual fish (N = 13, 17, 15, 8, and 10, from left to right), with the mean and standard deviation indicated by the square markers and error bars. AC, Acinetobacter calcoaceticus; AE, Aeromonas sp. strain ZOR0001; EN, Enterobacter sp. strain ZOR0014; PL, Plesiomonas sp. strain ZOR0011; PS, Pseudomonas mendocina.
The ability to quantify species abundance and to manipulate it by controlled addition or subtraction of species is commonplace in macroscopic ecological investigations. Its implementation here enables connections between intestinal microbiome research and a large literature on ecosystem dynamics. An issue whose importance has been realized for decades is the extent to which interspecies interactions are pairwise additive or whether higher-order (often called indirect) interactions are necessary to explain community structure (41, 42). The term “higher-order interactions” has been defined in various ways in the ecological literature (42, 43), in some cases referring specifically to nonadditive changes in a species’ growth rate given the presence of additional species, or to changes in the nature of the interaction between two species induced by additional species. In other cases it refers more generally to any interaction that cannot be captured by a pairwise model. We adopt the latter, commonly used definition, which is agnostic to underlying mechanisms (44). In our analyses below, we consider various pairwise models and assess their ability to describe data from multispecies communities; mismatch is indicative of the existence of higher-order interactions. Pairwise additivity, if dominant, simplifies the prediction of ecosystem composition, which would be desirable for therapeutic applications of microbiome engineering. Higher-order interactions may stabilize multispecies communities according to several recent theoretical models described further in Discussion (4548), implying that quantifying and controlling indirect effects is necessary for reshaping gut microbiomes.
Whether host associated or not, microbial communities have shown a variety of interaction types. A classic study involving cultured protozoan species found good agreement between the dynamics of four-species consortia and predictions derived from measurements of pairs of species (49). Similarly, Friedman and colleagues showed that the outcomes of competitions among three-species communities of soil-derived bacteria could be predicted simply from the outcomes of pairwise combinations (50). In contrast, experiments based on the cheese rind microbiome found significant differences in the genes required for a nonnative Escherichia coli species to persist in a multispecies bacterial community compared to predictions from pairwise coexistence with community members (51). A closed ecosystem consisting of one species each of algae, bacteria, and ciliate exhibited a strong nonpairwise interaction, in which the bacteria is abundant in the presence of each of the algae or ciliate alone but is subject to strong predation in the three-species system (52).
Within animals, the interaction types observed in the few studies to date that make use of controlled microbial communities in gnotobiotic hosts are also disparate. Competitive outcomes of three-species communities from subsets of 11 different bacterial species in the gut of the nematode C. elegans could be predicted from the outcomes of two-species experiments, with indirect effects found to be weak (31). In contrast, work using well-defined bacterial assemblies of up to five species in the fruit fly D. melanogaster found strong higher-order interactions governing microbe-dependent effects on host traits, such as life span (32).
To our knowledge, there have been no quantitative assessments of interbacterial interactions using controlled combinations of microbial species in a vertebrate host, leaving open the question of whether higher-order interactions are strong or whether pairwise characterizations suffice to predict intestinal community structure. Therefore, we examined larval zebrafish, inoculating initially germfree animals with specific subsets of five different species of zebrafish-derived bacteria and assessing their subsequent absolute abundances. Although the number of species is considerably smaller than the hundreds that may be present in a normal zebrafish intestine, it is large enough to sample a range of higher-order interactions yet small enough that the number of permutations of species is tractable.
As detailed below, we find strong pairwise interactions between certain bacterial species. However, we find weaker interactions and a greater than expected level of coexistence in fish colonized by four or five bacterial species. This suggests that measurements of pairwise intermicrobial interactions are insufficient to predict the composition of multispecies gut communities and that higher-order interactions dampen strong competition and facilitate diversity in a vertebrate intestine.


Zebrafish (Fig. 1A) were derived to be germfree and then were inoculated at 5 days postfertilization (dpf) with the desired combination of microbial species by addition of bacteria to the flasks housing the fish. Approximately 48 h later, fish were euthanized and their intestines removed by dissection. Intestines and their contents were homogenized, diluted, and plated onto chromogenic agar (see Materials and Methods). Secreted enzymes from each of the five candidate bacterial species generate particular colors due to substrates in the chromogenic medium, allowing the quantification of the number of CFU and, therefore, absolute intestinal abundance (Fig. 1B). All abundance data are provided in Data Set S1 in the supplemental material.
The five species examined were selected as diverse representatives of genera commonly found in the zebrafish intestine. Full names and species identifiers are given in Materials and Methods. As expected given their association with the zebrafish gut microbiome, each species in monoassociation, i.e., as the sole species inoculated in germfree fish, colonizes robustly to an abundance of 103 to 104 CFU/gut, corresponding to an in vivo density of approximately 109 to 1010 bacteria/ml (Fig. 1C).

Pairwise interactions in diassociations.

We first examined all 10 possible coinoculations of two species, which enables the assessment of pairwise interactions in the absence of higher-order effects. Intestinal CFU data show a wide range of outcomes for different species pairs. As exemplars, the number of CFU per gut for each of two species, Acinetobacter calcoaceticus and Enterobacter sp. strain ZOR0014, in the presence of each of the other four is displayed in Fig. 2A and B, respectively. The abundance of A. calcoaceticus is similar in the presence of any second species to its value in monoassociation. In contrast, the mean Enterobacter sp. strain ZOR0014 abundance is similar to its monoassociation value if coinoculated with Plesiomonas sp. strain ZOR0011 or Pseudomonas mendocina, about 10 times lower if coinoculated with A. calcoaceticus, and over 2 orders of magnitude lower if coinoculated with Aeromonas sp. strain ZOR0001, implying in the latter cases strong negative interactions.
FIG 2 Strong negative pairwise interactions dominate diassociation experiments. Abundances per zebrafish gut of A. calcoaceticus (A) and Enterobacter sp. strain ZOR0014 (B) in monoassociation (gray) and in diassociation with each of the other bacterial species (blue/green). Each circular data point is a CFU value from an individual fish (N = 13, 21, 19, 20, and 27 [A] and N = 15, 19, 22, 18, and 23 [B], from left to right), with the means and standard deviations indicated by the square markers and error bars. (C) Matrix of pairwise interaction coefficients, CijII, characterizing the effect of species j on the abundance of species i. Coefficients that differ from zero by more than three standard deviations (provided in Fig. S2A) are outlined in black. (D) The average bacterial load per zebrafish in each of the diassociation combinations, expressed as log10 total CFU. The standard deviations are between 0.3 and 1.1 and are displayed in Fig. S1. Values on the diagonal are the monoassociation load for each of the five species. AC, Acinetobacter calcoaceticus; AE, Aeromonas sp. strain ZOR0001; EN, Enterobacter sp. strain ZOR0014; PL, Plesiomonas sp. strain ZOR0011; PS, Pseudomonas mendocina.
Parameterizing the strength of interactions between species is necessarily model dependent, contingent on the functional form of the relationship between one species’ abundance and that of the other. We show that the conclusions we reach regarding interaction strengths, especially their shifts when multiple species are present, are qualitatively similar for a wide range of models and, therefore, robust. We first consider a phenomenological interaction coefficient, CijII, that is linear in log abundance, characterizing the effect of species j on species i as
where Pi denotes the abundance of species i and the superscript I or II denotes a mono- or diassociation experiment. This form is motivated by the distribution of gut bacterial abundances being roughly lognormal, with species addition capable of inducing orders-of-magnitude changes (Fig. 2A and B). This CijII value can be derived as the interaction parameter in a competitive Lotka-Volterra model modified to act on log abundances (see Text S1). Qualitatively, a positive Cij implies that the abundance of species i increases in the presence of j. Similarly, a negative Cij indicates that the abundance of species i declines in the presence of species j. Subsampling from the measured sets of bacterial abundances gives the mean and standard deviation of the estimated interaction parameters (Text S1).
We plot in Fig. 2C the CijII defined by Equation 1, calculated from all diassociation data of all species pairs (N = 190 fish in total). For determining CijII, we only use data from fish in which both species were detected so that abundance changes of one species can definitively be ascribed to the presence of the other within the gut. Uncertainties in CijII are estimated from bootstrap subsampling (Text S1). The interactions are predominantly negative. Thirteen out of 20 coefficients differ from zero by over three standard deviations, indicating both a large magnitude and a less than 0.001 probability that the interaction strength is zero or of the opposite sign. The total bacterial load, i.e., the sum of the bacterial abundances, is similar for all the diassociations, suggesting that the interaction effects do not stem from changes in intestinal capacity (Fig. 2D).
Although the physical and chemical environment of the zebrafish gut is likely very dissimilar to test tubes of standard growth media, we examined abundances of each of the pairs of species in in vitro competition experiments, growing overnight cultures in lysogeny broth (LB) medium and plating for CFU (see Materials and Methods). Assessing CijII as described above, we find, as expected, that interaction coefficients calculated from the in vitro experiments are markedly different from those measured in vivo (Fig. 3B and Fig. S3).
FIG 3 Weak pairwise interactions in five-species experiments. (A) Abundance per zebrafish gut of one of the bacterial species, Enterobacter sp. strain ZOR0014, when all five species are coinoculated (gray) and in each four species coinoculation experiments (green), with the omitted species indicated on the axis. Each circular data point is a CFU value from an individual fish (N = 40, 12, 12, 11, and 9, from left to right), with the mean and standard deviation indicated by the square marker and error bars. (B) Matrix of pairwise interaction coefficients, CijV, when 5 bacterial species are present. The coefficients outlined in black differ from zero by over three standard deviations (see Fig. S2B). (C) The pairwise interaction coefficients inferred from 4- to 5-species experiments versus those from 1- to 2-species experiments. The colors label species i for each interaction pair. (D) The minimum interaction coefficient calculated from a power law interaction model for different values of the exponent α for the 1- to 2-species (square filled markers) and the 4- to 5-species (square markers) experiments. AC, Acinetobacter calcoaceticus; AE, Aeromonas sp. strain ZOR0001; EN, Enterobacter sp. strain ZOR0014; PL, Plesiomonas sp. strain ZOR0011; PS, Pseudomonas mendocina.
Our characterizations of interactions within the zebrafish gut are not qualitatively altered by using a more general power law model to compute CijII from absolute abundance data (see “Interactions under more general models,” below) following the presentation of measurements of interactions between more than two species.

Pairwise interactions in multispecies communities.

To assess whether the strong competitive interactions we found in two-species experiments are conserved in multispecies communities, we quantified pairwise interactions in experiments inoculating fish with four or five bacterial species. To assess CijV, we adopted a method similar to the leave-one-out approach often used in macroscopic ecological studies, dating at least to classic experiments in which single species were removed from tide pools and the abundances of the remaining species were measured to evaluate interspecies interactions (53). Here, we performed coinoculation experiments, leaving out one of the five species of bacteria, and compared intestinal abundances for these four-species communities to those measured in five-species coinoculation experiments.
In approximately N = 10 fish each, we performed all five different coinoculations of four bacterial species. The difference in the abundance of species i in fish inoculated with all five species compared to fish inoculated with four, missing species j, gives a measure of the impact of species j on species i in the multispecies environment. As an example, Enterobacter sp. strain ZOR0014 abundance in inoculations lacking A. calcoaceticus, Aeromonas sp. strain ZOR0001, Plesiomonas sp. strain ZOR0011, and P. mendocina, and its abundance in five-species inoculations, are shown in Fig. 3A. In contrast to the diassociation experiments (Fig. 2B), we see that Enterobacter sp. strain ZOR0014 does not show large abundance differences, in either its mean or its distribution, as a result of any fifth species being present. Independent of any model, this suggests that nonpairwise, i.e., higher-order, interactions are present in the multispecies community.
Again, a variety of options are possible for quantifying interaction coefficients in the multispecies system. We first consider interaction coefficients as modifying mean log abundances, analogous to the pairwise model of Equation 1:
The interaction coefficients, CijV, that we obtain, displayed in Fig. 3B, are different and in general are considerably weaker than the CijII found in the two-species case (Fig. 2C). There are only three interactions that differ from zero by over three standard deviations. Strikingly, all three of these interactions are positive. This shift toward weaker and more positive interactions between the two-species and multispecies interactions is further illustrated in Fig. 3C, in which the multispecies interaction coefficients, CijV, are plotted against the 2-species interaction coefficients, CijII.

Interactions under more general models.

As noted, a model that is linearly additive in logarithmic abundances is only one of an infinite number of choices and, moreover, may not adequately capture the complexity of interactions in the gut. Earlier experiments investigating the spatial structure of specific microbial communities in the larval zebrafish intestine have shown that species such as Aeromonas sp. strain ZOR0001, Enterobacter sp. strain ZOR0014, and P. mendocina form dense three-dimensional aggregates (54). The size and location of aggregates and the locations of cells, conspecific or otherwise, within these aggregates may affect their interactions in ways that could be sublinear, linear, or superlinear in population size. Previous work has also established that gut bacteria also influence intestinal mechanics (40), highlighting one of many possible indirect interaction mechanisms whose functional forms are unknown. Furthermore, other studies have shown that different modes of physical and chemical communication could result in long-range interactions between different species (5557). To address these possibilities, we evaluated species interactions with a more general power law model, wherein the interaction effects between species could be nonlinear in the abundance of the effector species. Here, the interaction coefficient Cij depends on a power, α, of the abundance of the effector species j, which we evaluate in the range of α = 0.1 to 2, spanning sublinear and superlinear interactions. Modified versions of Equations 1 and 2 give
from which we can evaluate CijII and CijV, respectively. Note that α = 1 in Equations 3 and 4, i.e., interactions that are linear in abundance, is simply the steady-state behavior of the competitive Lotka-Volterra model commonly used in population modeling and is shown in Fig. S5. We provide the CijII and CijV for several different α values in Fig. S4. Throughout, as in the logarithmic model shown above, pairwise interactions in diassociation are, in many cases, strongly negative, while the multispecies interactions are weaker. This is summarized by studying the trends in the most negative CijII and CijV for different values of α, depicted in Fig. 3D, which shows that for all α values, the strongest CijV is significantly weaker than the strongest CijII, suggesting that our results are robust to choice of model.

Five-species coexistence.

We next consider coinoculation of all five bacterial species. Examination of over 200 fish shows a large variety in abundances, depicted in Fig. 4A as the relative abundance of each species in each larval gut. Multiple species are able to coexist, with the median number of species present being 4 (Fig. 4B). The mean total bacterial load as well as its distribution (Fig. 4C) are similar to the mean and distribution of the mono- and diassociation experiments, as well as four-species coinoculation experiments discussed earlier. We calculated the expected abundance of each bacterial species if the interactions governing the five-species community were simply a linear combination of the pairwise interactions governing diassociations, CijII. Any of the additive models we evaluated can be extended to combinations of species. Considering the model focused on above, with interaction coefficients modifying log abundances, the predicted abundance of species i in the presence of another species j is given by
FIG 4 Communities are more diverse and abundant in five-species experiments than would be predicted solely based on two-species pairwise interactions. (A) Stacked bar plot of the relative abundances of the five bacterial species when all five were coinoculated. Each bar is from a single dissected fish. The bars are ordered by total bacterial load. (B) Histogram of the total number of bacterial species present in the gut when all five species were coinoculated. (C) The total bacterial load as a function of the number of inoculated species. Each circular data point is a CFU value from an individual fish (N = 63, 232, 187, and 202, from left to right), with the mean and standard deviation indicated by the square marker and error bars. (D) The predicted (blue Xs) and measured (brown circles) abundances of each bacterial species in the zebrafish gut when all five species are coinoculated. Predictions are based on an interaction model that is linear in log abundance using the pairwise CijII coefficients, as described in the text. Solid square markers indicate the mean and standard deviation of the distributions excluding null counts. The dotted line indicates the experimental detection limit of 25 cells. The experimental data are from N = 202 fish in total, and the predicted distributions arise from 250 samples of the distribution of interaction coefficients. (E) The observed frequency of occurrence in the gut from the five-species coinoculation experiment versus the predicted frequencies for each of the five species. (F) The Pearson correlation coefficients calculated from the relative abundances of pairs of species when all five species were coinoculated. AC, Acinetobacter calcoaceticus; AE, Aeromonas sp. strain ZOR0001; EN, Enterobacter sp. strain ZOR0014; PL, Plesiomonas sp. strain ZOR0011; PS, Pseudomonas mendocina.
where the superscript V denotes the five-species coinoculation experiment. A model linear in species abundance (α = 1) is also considered in the Text S1 and gives qualitatively similar outputs and conclusions. Sampling from the measured distributions of each of the interaction coefficients and the mean abundance in monoassociation allows the calculation of the distribution of expected PiV values (see Data Set S3).
We plot the measured and predicted distributions of intestinal abundances of each of the five species for the five-species coinoculation experiment in Fig. 4D. The measured distributions of each of the species are very similar. In contrast, the distributions of the predicted abundances vary significantly by species. For two of the species, A. calcoaceticus and P. mendocina, the means of the observed and predicted distributions are similar. For the other three species, in contrast, the observed and predicted populations are in strong disagreement, with the pairwise prediction being at least an order of magnitude lower than the observed abundances. For Enterobacter sp. strain ZOR0014 and Plesiomonas sp. strain ZOR0011 in particular, we would expect extinction in a large fraction of fish due to strong negative pairwise interactions; in actuality, both species are common and abundant.
Similarly, we can extract from the model the predicted frequency of occurrence of each of the species, regardless of abundance, i.e., the fraction of fish with a nonzero population of that species. We find that the predicted frequency is much lower than the experimentally observed frequency for Plesiomonas sp. strain ZOR0011 and Enterobacter sp. strain ZOR0014 (Fig. 4E).
By measuring absolute abundances of bacterial populations in the gut, we provide direct assessments of interspecies interactions. More common sequencing-based methods, applied, for example, to the human gut microbiome, typically provide relative measures of species abundance, i.e., each taxonomic unit’s fraction of the total load. Correlations among relative abundances are often used as measures of interaction strengths (58, 59). Calculating the Pearson correlation coefficients of the relative abundances of each pair of species in fish inoculated with all five bacterial species, we find a strikingly different interaction matrix (Fig. 4F) than that inferred from absolute abundance changes (Fig. 3B), with many strong negative values. There are many likely reasons for the difference between Pearson correlations and our more directly measured interaction coefficients. The Pearson r necessarily attributes correlations between pairs of species as being indicative of the dynamics of that pair independent of other species, is confounded by overall changes in total bacterial load, and, perhaps most importantly, is necessarily symmetric (Cij = Cji). Our Cij values, inferred from absolute abundance data, are notably asymmetric (Fig. 2C and 3B).


Using a model system comprising five commensal bacterial species in the larval zebrafish intestine, we have characterized aspects of gut microbiome assembly. Controlled combinations of inoculated species and measurements of absolute abundance in the gut, both challenging to perform in other vertebrate systems, reveal clear signatures of interactions among species. We find strong, competitive interactions among certain pairs in fish inoculated with two bacterial species. In contrast, pairwise interactions are weak in intestines colonized by four to five species, and all species are present at equal or greater abundance than would be predicted based on two-species data.
Our quantification of interaction strengths relies on a minimal set of assumptions that serve as a general test of additive models. Interaction strengths are necessarily parameters of particular models. We made use of a model in which the log-transformed population of a species is a linear function of the other species’ log-transformed populations and a more general power law model that spans both sublinear and superlinear dependences on population sizes. There are good reasons to be skeptical of such frameworks. First, intestinal populations may not be well described by equilibrium, steady-state values. Second, these models lack spatial structure information. In vivo microscopy of one or two species in the zebrafish gut (39, 40, 60) underscores both of these concerns: populations are very dynamic, with rapid growth and stochastic expulsions; interactions can be mediated by complex intestinal mechanics; and aggregation and localization behaviors are species specific.
Imaging, however, provides justifications for these rough models. Prior microscopy-based studies have shown that growth rates are rapid, with populations reaching carrying capacities within roughly 12 h (39, 60), well below the 48-h assessment time considered here. Because of strong aggregation observed in nearly all bacterial species, most individual bacteria residing in the bulk of clusters will not directly interact with other species, leading to interactions that are sublinear in population size, suggesting a logarithmic or α < 1 power law functional form. Furthermore, stochastic dynamics can be mapped onto robust average properties for populations (39, 61). Therefore, it is reasonable to make use of simple models not as rigorous descriptions of the system but as approximations whose parameters characterize effective behaviors. We note that all these issues also affect more commonly used models, such as standard competitive Lotka-Volterra models that are linear in population sizes. These models are often applied to gut microbiome data and used to infer interaction parameters (26, 62, 63) despite a lack of information about their realism. The power law model of interactions provides the strongest indication of the generality of our conclusions. Over a range of interaction forms extending from highly sublinear (α = 0.1) to superlinear (α = 2.0), strong competitive interactions are dampened when four or five species are present (Fig. 3D), suggestive of higher-order interactions among intestinal bacteria.
The ecological potential for higher-order or indirect interactions, i.e., interactions that cannot be reduced to pairwise additive components but rather result from the activities of three or more species, to be important determinants of community structure has been appreciated for decades (41, 42, 49). Identification of higher-order interactions among constituent species is important for accurate prediction of responses to ecological perturbations such as species invasion or extinction, as well as functions of multispecies communities, as such features will not be adequately forecast by the examination of direct interactions in subsystems (41, 64).
Inferring and quantifying indirect interactions in natural ecosystems is, however, challenging, calling for subtle and model-dependent statistical tests (41, 42, 65). Constructed or manipulated systems enable more straightforward assays in which particular species are introduced or removed amid a backdrop of others. Several such systems involving macroscopic organisms (6670), as well as microorganisms (32, 52), have uncovered significant indirect interactions. However, some studies of microbial communities have found weak or negligible higher-order interactions (49, 50), including one study examining combinations of species introduced to the C. elegans gut (31). The complexity of interactions in a vertebrate gut has remained unclear, and correlation-based methods for inferring interactions from sequencing-based data have assumed that pairwise interactions suffice (58, 71, 72).
Our measurements using gnotobiotic larval zebrafish, a model vertebrate, show strong pairwise interactions when only two bacterial species are present in the intestine and weak pairwise interactions when four to five species are present, indicating that higher-order interactions are important (Fig. 3). In many cases, the effect is evident from the raw data itself. For example, Enterobacter sp. strain ZOR0014 is strongly suppressed by Aeromonas sp. strain ZOR0001 if the two are inoculated together (Fig. 2B). Comparing Enterobacter sp. strain ZOR0014 abundance in fish colonized by all species except Aeromonas sp. strain ZOR0001 with its abundance in fish colonized by all species, in contrast, shows little difference, indicating that the Enterobacter sp. strain ZOR0014-Aeromonas sp. strain ZOR0001 interaction is strongly attenuated by the presence of the other bacterial groups (Fig. 3A).
Two additional observations also imply the presence of strong higher-order interactions in our intestinal ecosystem. Considering fish colonized by all five bacterial species, the mean abundance of each species is at least as high as would be predicted solely from direct interactions (Fig. 4D). Moreover, the diversity of bacterial species is higher than would be predicted, as all of the species occur in more than 50% of fish, contrary to prediction (Fig. 4E). Considering fish colonized by all five bacterial species, the abundance of each species is at least as high and the diversity of bacterial species is higher than the values that would be predicted solely from direct interactions.
Our finding of higher species diversity than expected from pairwise interactions in a system of several gut bacterial species is consistent with recent theoretical studies that suggest, for a variety of reasons, that higher-order interactions are likely to stabilize communities and promote coexistence. The topic of multispecies coexistence has a long history in ecology. Especially since classic work by Robert May showing that a system comprising pairwise interacting constituents will, in general, be less stable as the number of species increases (73), explaining how complex communities can exist has been a theoretical challenge. There are many resolutions to this paradox, such as spatial heterogeneity, interactions across trophic levels, and temporal variation. However, even without such additional structure, incorporating higher-order terms into general random competitive interaction models leads to widespread coexistence (4547). Such large-scale coexistence can also emerge naturally from contemporary resource competition models (48, 74), in which cross-feeding or metabolic tradeoffs necessarily involve multiple interacting species. Intriguingly, the abundance distributions of all five of our gut bacterial species, when inoculated together, are similar to one another. The average Shannon entropy of the five-species community (H = 1.16 ± 0.24) (see Text S1 in the supplemental material) also resembles that of a purely neutrally assembled community (H = 1.61), reminiscent of dynamics mimicking neutral assembly that emerge from multispecies dynamics driven by resource use constraints (48, 75).
Our findings imply that measurements of two-species interactions among microbial residents of the vertebrate gut are insufficient for predictions of community dynamics and composition. Moreover, they imply that inference from microbiome data of interspecies interactions, for example, by fitting Lotka-Volterra-type models with pairwise interaction terms (26, 62, 63, 76), should not be thought of as representing fundamental pairwise interactions that would be manifested if, for example, the constituent species were isolated but rather as effective interactions in a complex milieu.
Our measurements do not shed light on what mechanisms give rise to higher-order interactions in our system. Likely candidates include metabolic interactions among the species, interactions mediated by host activities, such as immune responses, and modulation of spatial structure by coexisting species. Immune responses are sensitive to specific bacterial species (77) and to bacterial behaviors (78). Regarding spatial structure in particular, in vivo imaging of these bacterial species in monoassociation has shown robust aggregation behaviors that correlate with location in the gut (54). Given the physical constraints of the intestinal environment, we think that the modification of spatial organization due to the presence of species with overlapping distributions is a likely mechanism for higher-order interactions. Notably, both immune responses and spatial structure are amenable to live imaging in larval zebrafish (39, 40, 54). Although the parameter space of transgenic hosts, fluorescent labels, and interaction timescales to explore in imaging studies is potentially very large, such future studies are likely to yield valuable insights into the mechanisms orchestrating the strong interactions observed here. Furthermore, examination of the roles of priority effects and other aspects of initial colonization, as well as the stability of diverse communities with respect to invasion, may reveal routes for intentionally manipulating the vertebrate microbiome to engineer desired traits.


Animal care.

All experiments with zebrafish were done in accordance with protocols approved by the University of Oregon Institutional Animal Care and Use Committee and by following standard protocols (79).


Wild-type (ABC×TU strain) larval zebrafish (Danio rerio) were derived germfree as described in reference 36. In brief, embryos were washed at approximately 7 h postfertilization with antibiotic, bleach, and iodine solutions and then moved to tissue culture flasks containing 15 ml sterile embryo medium solution with approximately 1 ml of sterile solution per larva. The flasks were then stored in a temperature-controlled room maintained at 28°C.

Bacterial strains and culture.

The five bacterial strains used in this study, namely, Aeromonas sp. strain ZOR0001, Pseudomonas mendocina (ZWU0006), Acinetobacter calcoaceticus (ZOR0008), Enterobacter sp. strain ZOR0014, and Plesiomonas sp. strain ZOR0011 were originally isolated from the zebrafish intestine and have been fluorescently labeled to express green fluorescent protein (GFP) and dTomato, facilitating their identification in our experimental assays (80, 81). Stocks of bacteria were maintained in 25% glycerol at −80°C.

Inoculation of tissue culture flasks.

One day prior to inoculation of the tissue culture flasks, bacteria from frozen glycerol stocks were shaken overnight in lysogeny broth (LB medium; 10 g/liter NaCl, 5 g/liter yeast extract, 12 g/liter tryptone, 1 g/liter glucose) and grown for 16 h overnight at 30°C. Samples of 1 ml of each of the overnight cultures were washed twice by centrifuging at 7,000 × g for 2 min, removing the supernatant, and adding 1 ml of fresh sterile embryo medium. At 5 dpf, the tissue culture flasks were inoculated with this solution at a concentration of 106 CFU/ml. For each of the competition experiments involving 2, 4, and 5 bacterial species, equal concentrations were inoculated into the flasks. After inoculation, the flasks were maintained at 30°C until dissection at 7 dpf.

Dissection and plating.

To determine the intestinal abundance of bacterial species, dissections of larval zebrafish were performed at 7 dpf. Zebrafish were euthanized by hypothermal shock. Intestines were removed by dissection, placed in 500 μl of sterile embryo medium, and homogenized with zirconium oxide beads using a bullet blender. The homogenized gut solution was diluted to 10−1 and 10−2, and 100-μl volumes of these dilutions were spread onto agar plates. For mono- and diassociated inoculations, tryptic soy agar (TSA) plates were used in which fluorescence could be used to differentiate up to two species. For inoculations of more than two species, Universal HiChrome agar (Sigma-Aldrich) plates were used, allowing for visual differentiation of each species using a colorimetric indicator. The abundances of each of the species in the zebrafish gut were determined by counting the numbers of CFU on the plates. These abundances for different experiments are provided in Data Set S1.

In vitro competition experiments.

To determine the in vitro competition coefficients, all of the different pairwise combinations of the five species were grown in overnight cultures of LB medium as described above. On the following day, cultures were plated at 10−7 or 10−6 dilutions, depending on the ability to detect both species in a given dilution. Abundances were obtained by counting the number of CFU of each species on the plates. These values are provided in Data Set S1.


We thank Rose Sockol and the University of Oregon Zebrafish Facility staff for fish husbandry and Laura-Taggart Murphy for germfree derivations of zebrafish larvae. We also thank Jayson Paulose and Karen Guillemin for insightful comments and suggestions.
Research was supported by an award from the Kavli Microbiome Ideas Challenge, a project led by the American Society for Microbiology in partnership with the American Chemical Society and the American Physical Society and supported by The Kavli Foundation. Work was also supported by the National Science Foundation under awards 1427957 (R.P.) and 0922951 (R.P.), and the National Institutes of Health (, under awards P50GM09891 and P01GM125576-01. The content is solely the responsibility of the authors and does not represent the official views of the NSF, National Institutes of Health, or other funding agencies.

Supplemental Material

File (mbio.01667-20-s0001.pdf)
File (mbio.01667-20-sd001.xlsx)
File (mbio.01667-20-sd002.xlsx)
File (mbio.01667-20-sd003.xlsx)
File (mbio.01667-20-sf001.pdf)
File (mbio.01667-20-sf002.pdf)
File (mbio.01667-20-sf003.pdf)
File (mbio.01667-20-sf004.pdf)
File (mbio.01667-20-sf005.pdf)
File (mbio.01667-20-sf006.pdf)
ASM does not own the copyrights to Supplemental Material that may be linked to, or accessed through, an article. The authors have granted ASM a non-exclusive, world-wide license to publish the Supplemental Material files. Please contact the corresponding author directly for reuse.


Semova I, Carten JD, Stombaugh J, Mackey LC, Knight R, Farber SA, Rawls JF. 2012. Microbiota regulate intestinal absorption and metabolism of fatty acids in the zebrafish. Cell Host Microbe 12:277–288.
Martinez-Guryn K, Hubert N, Frazier K, Urlass S, Musch MW, Ojeda P, Pierre JF, Miyoshi J, Sontag TJ, Cham CM, Reardon CA, Leone V, Chang EB. 2018. Small intestine microbiota regulate host digestive and absorptive adaptive responses to dietary lipids. Cell Host Microbe 23:458–469.
Hill JH, Franzosa EA, Huttenhower C, Guillemin K. 2016. A conserved bacterial protein induces pancreatic beta cell expansion during zebrafish development. Elife 5:e20145.
Yan J, Herzog JW, Tsang K, Brennan CA, Bower MA, Garrett WS, Sartor BR, Aliprantis AO, Charles JF. 2016. Gut microbiota induce IGF-1 and promote bone formation and growth. Proc Natl Acad Sci U S A 113:E7554–E7563.
Mazmanian SK, Liu CH, Tzianabos AO, Kasper DL. 2005. An immunomodulatory molecule of symbiotic bacteria directs maturation of the host immune system. Cell 122:107–118.
Kamada N, Núñez G. 2014. Regulation of the immune system by the resident intestinal bacteria. J Gastroenterol 146:1477–1488.
Arpaia N, Campbell C, Fan X, Dikiy S, van der Veeken J, deRoos P, Liu H, Cross JR, Pfeffer K, Coffer PJ, Rudensky AY. 2013. Metabolites produced by commensal bacteria promote peripheral regulatory T-cell generation. Nature 504:451–455.
Honda K, Littman DR. 2012. The microbiome in infectious disease and inflammation. Annu Rev Immunol 30:759–795.
Scharschmidt TC, Vasquez KS, Truong HA, Gearty SV, Pauli ML, Nosbaum A, Gratz IK, Otto M, Moon JJ, Liese J, Abbas AK, Fischbach MA, Rosenblum MD. 2015. A wave of regulatory T cells into neonatal skin mediates tolerance to commensal microbes. J Immun 43:1011–1021.
Maynard CL, Elson CO, Hatton RD, Weaver CT. 2012. Reciprocal interactions of the intestinal microbiota and immune system. Nature 489:231–241.
Boulangé CL, Neves AL, Chilloux J, Nicholson JK, Dumas ME. 2016. Impact of the gut microbiota on inflammation, obesity, and metabolic disease. Genome Med 8:42.
Pflughoeft KJ, Versalovic J. 2012. Human microbiome in health and disease. Annu Rev Pathol 7:99–122.
Wang B, Yao M, Lv L, Ling Z, Li L. 2017. The human microbiota in health and disease. Engineering 3:71–82.
Moloney RD, Desbonnet L, Clarke G, Dinan TG, Cryan JF. 2014. The microbiome: stress, health and disease. Mamm Genome 25:49–74.
Shreiner AB, Kao JY, Young VB. 2015. The gut microbiome in health and in disease. Curr Opin Gastroenterol 31:69–75.
Rajilić Stojanović M, Biagi E, Heilig HGHJ, Kajander K, Kekkonen RA, Tims S, de Vos WM. 2011. Global and deep molecular analysis of microbiota signatures in fecal samples from patients with irritable bowel syndrome. J Gastroenterol 141:1792–1801.
Durack J, Lynch SV. 2019. The gut microbiome: relationships with disease and opportunities for therapy. J Exp Med 216:20–40.
Sokol H, Pigneur B, Watterlot L, Lakhdari O, Bermúdez-Humarán LG, Gratadoux JJ, Blugeon S, Bridonneau C, Furet JP, Corthier G, Grangette C, Vasquez N, Pochart P, Trugnan G, Thomas G, Blottière HM, Doré J, Marteau P, Seksik P, Langella P. 2008. Faecalibacterium prausnitzii is an anti-inflammatory commensal bacterium identified by gut microbiota analysis of Crohn disease patients. Proc Natl Acad Sci U S A 105:16731–16736.
Zhu W, Winter MG, Byndloss MX, Spiga L, Duerkop BA, Hughes ER, Büttner L, de Lima Romão E, Behrendt CL, Lopez CA, Sifuentes-Dominguez L, Huff-Hardy K, Wilson RP, Gillis CC, Tükel A, Koh AY, Burstein E, Hooper LV, Bäumler AJ, Winter SE. 2018. Precision editing of the gut microbiota ameliorates colitis. Nature 553:208–211.
Everard A, Belzer C, Geurts L, Ouwerkerk JP, Druart C, Bindels LB, Guiot Y, Derrien M, Muccioli GG, Delzenne NM, de Vos WM, Cani PD. 2013. Cross-talk between Akkermansia muciniphila and intestinal epithelium controls diet-induced obesity. Proc Natl Acad Sci U S A 110:9066–9071.
Dethlefsen L, Relman DA. 2011. Incomplete recovery and individualized responses of the human distal gut microbiota to repeated antibiotic perturbation. Proc Natl Acad Sci U S A 108:4554–4561.
Modi SR, Collins JJ, Relman DA. 2014. Antibiotics and the gut microbiota. J Clin Investig 124:4212–4218.
Zhao L, Zhang F, Ding X, Wu G, Lam YY, Wang X, Fu H, Xue X, Lu C, Ma J, Yu L, Xu C, Ren Z, Xu Y, Xu S, Shen H, Zhu X, Shi Y, Shen Q, Dong W, Liu R, Ling Y, Zeng Y, Wang X, Zhang Q, Wang J, Wang L, Wu Y, Zeng B, Wei H, Zhang M, Peng Y, Zhang C. 2018. Gut bacteria selectively promoted by dietary fibers alleviate type 2 diabetes. Science 359:1151–1156.
Ye L, Mueller O, Bagwell J, Bagnat M, Liddle RA, Rawls JF. 2019. High fat diet induces microbiota-dependent silencing of enteroendocrine cells. Elife 8:e48479.
Carmody RN, Gerber GK, Luevano JM, Gatti DM, Somes L, Svenson KL, Turnbaugh PJ. 2015. Diet dominates host genotype in shaping the murine gut microbiota. Cell Host Microbe 17:72–84.
Fisher CK, Mehta P. 2014. Identifying keystone species in the human gut microbiome from metagenomic timeseries using sparse linear regression. PLoS One 9:e102451.
Shang Y, Sikorski J, Bonkowski M, Fiore-Donno AM, Kandeler E, Marhan S, Boeddinghaus RS, Solly EF, Schrumpf M, Schöning I, Wubet T, Buscot F, Overmann J. 2017. Inferring interactions in complex microbial communities from nucleotide sequence data and environmental parameters. PLoS One 12:e0173765.
Carr A, Diener C, Baliga NS, Gibbons SM. 2019. Use and abuse of correlation analyses in microbial ecology. ISME J 13:2647–2655.
Lovell D, Pawlowsky-Glahn V, Egozcue JJ, Marguerat S, Bähler J. 2015. Proportionality: a valid alternative to correlation for relative data. PLoS Comput Biol 11:e1004075.
Zmora N, Zilberman-Schapira G, Suez J, Mor U, Dori-Bachash M, Bashiardes S, Kotler E, Zur M, Regev-Lehavi D, Brik RBZ, Federici S, Cohen Y, Linevsky R, Rothschild D, Moor AE, Ben-Moshe S, Harmelin A, Itzkovitz S, Maharshak N, Shibolet O, Shapiro H, Pevsner-Fischer M, Sharon I, Halpern Z, Segal E, Elinav E. 2018. Personalized gut mucosal colonization resistance to empiric probiotics is associated with unique host and microbiome features. Cell 174:1388–1405.
Lopez AO, Vega NM, Gore J. 2019. Interspecies bacterial competition determines community assembly in the C. elegans intestine. bioRxiv.
Gould AL, Zhang V, Lamberti L, Jones EW, Obadia B, Korasidis N, Gavryushkin A, Carlson JM, Beerenwinkel N, Ludington WB. 2018. Microbiome interactions shape host fitness. Proc Natl Acad Sci U S A 115:E11951–E11960.
Aranda-Díaz A, Obadia B, Dodge R, Thomsen T, Hallberg ZF, Güvener ZT, Ludington WB, Huang KC. 2020. Bacterial interspecies interactions modulate pH-mediated antibiotic tolerance. Elife 9:e51493.
Burns AR, Stephens WZ, Stagaman K, Wong S, Rawls JF, Guillemin K, Bohannan BJ. 2016. Contribution of neutral processes to the assembly of gut microbial communities in the zebrafish over host development. ISME J 10:655–664.
Flores EM, Nguyen AT, Odem MA, Eisenhoffer GT, Krachler AM. 2020. The zebrafish as a model for gastrointestinal tract-microbe interactions. Cell Microbiol 22:e13152.
Melancon E, De La Torre Canny SG, Sichel S, Kelly M, Wiles TJ, Rawls JF, Eisen JS, Guillemin K. 2017. Best practices for germ-free derivation and gnotobiotic zebrafish husbandry. Methods Cell Biol 138:61–100.
Dooley K, Zon LI. 2000. Zebrafish: a model system for the study of human disease. Curr Opin Genet Dev 10:252–256.
Rawls JF, Samuel BS, Gordon JI. 2004. Gnotobiotic zebrafish reveal evolutionarily conserved responses to the gut microbiota. Proc Natl Acad Sci U S A 101:4596–4601.
Wiles TJ, Jemielita M, Baker RP, Schlomann BH, Logan SL, Ganz J, Melancon E, Eisen JS, Guillemin K, Parthasarathy R. 2016. Host gut motility promotes competitive exclusion within a model intestinal microbiota. PLoS Biol 14:e1002517.
Logan SL, Thomas J, Yan J, Baker RP, Shields DS, Xavier JB, Hammer BK, Parthasarathy R. 2018. The Vibrio cholerae type VI secretion system can modulate host intestinal mechanics to displace gut bacterial symbionts. Proc Natl Acad Sci U S A 115:E3779–E3787.
Wootton JT. 1994. The nature and consequences of indirect effects in ecological communities. Annu Rev Ecol Syst 25:443–466.
Billick I, Case TJ. 1994. Higher order interactions in ecological communities: what are they and how can they be detected? Ecology 75:1530–1543.
Sanchez A. 2019. Defining higher-order interactions in synthetic ecology: lessons from physics and quantitative genetics. Cell Syst 9:519–520.
Mittelbach GG, McGill BJ. 2019. Community ecology. Oxford University Press, Oxford, United Kingdom.
Bairey E, Kelsic ED, Kishony R. 2016. High-order species interactions shape ecosystem diversity. Nat Commun 7:1–7.
Levine JM, Bascompte J, Adler PB, Allesina S. 2017. Beyond pairwise mechanisms of species coexistence in complex communities. Nature 546:56–64.
Grilli J, Barabás G, Michalska-Smith MJ, Allesina S. 2017. Higher-order interactions stabilize dynamics in competitive network models. Nature 548:210–213.
Posfai A, Taillefumier T, Wingreen NS. 2017. Metabolic trade-offs promote diversity in a model ecosystem. Phys Rev Lett 118:e028103.
Vandermeer JH. 1969. The Competitive structure of communities: an experimental approach with protozoa. Ecology 50:362–371.
Friedman J, Higgins LM, Gore J. 2017. Community structure follows simple assembly rules in microbial microcosms. Nat Ecol Evol 1:1–7.
Morin M, Pierce EC, Dutton RJ. 2018. Changes in the genetic requirements for microbial interactions with increasing community complexity. Elife 7:e37072.
Mickalide H, Kuehn S. 2019. Higher-order interaction between species inhibits bacterial invasion of a phototroph-predator microbial community. Cell Syst 9:521–533.
Paine RT. 1980. Food webs: linkage, interaction strength and community infrastructure. J Anim Ecol 49:667–685.
Schlomann BH, Wiles TJ, Wall ES, Guillemin K, Parthasarathy R. 2018. Bacterial cohesion predicts spatial distribution in the larval zebrafish intestine. Biophys J 115:2271–2277.
Ozgen VC, Kong W, Blanchard AE, Liu F, Lu T. 2018. Spatial interference scale as a determinant of microbial range expansion. Sci Adv 4:eaau0695.
Stubbendieck RM, Straight PD. 2016. Multifaceted interfaces of bacterial competition. J Bacteriol 198:2145–2155.
Prindle A, Liu J, Asally M, Ly S, Garcia-Ojalvo J, Süel GM. 2015. Ion channels enable electrical communication within bacterial communities. Nature 527:59–63.
Hoffmann C, Dollive S, Grunberg S, Chen J, Li H, Wu GD, Lewis JD, Bushman FD. 2013. Archaea and fungi of the human gut microbiome: correlations with diet and bacterial residents. PLoS One 8:e66019.
Steele JA, Countway PD, Xia L, Vigil PD, Beman JM, Kim DY, Chow CET, Sachdeva R, Jones AC, Schwalbach MS, Rose JM, Hewson I, Patel A, Sun F, Caron DA, Fuhrman JA. 2011. Marine bacterial, archaeal and protistan association networks reveal ecological linkages. ISME J 5:1414–1425.
Schlomann BH, Wiles TJ, Wall ES, Guillemin K, Parthasarathy R. 2019. Sublethal antibiotics collapse gut bacterial populations by enhancing aggregation and expulsion. Proc Natl Acad Sci U S A 116:21392–21400.
Schlomann BH. 2018. Stationary moments, diffusion limits, and extinction times for logistic growth with random catastrophes. J Theor Biol 454:154–163.
Stein RR, Bucci V, Toussaint NC, Buffie CG, Rätsch G, Pamer EG, Sander C, Xavier JB. 2013. Ecological modeling from time-series inference: insight into dynamics and stability of intestinal microbiota. PLoS Comput Biol 9:e1003388.
Marino S, Baxter NT, Huffnagle GB, Petrosino JF, Schloss PD. 2014. Mathematical modeling of primary succession of murine intestinal microbiota. Proc Natl Acad Sci U S A 111:439–444.
Sanchez-Gorostiaga A, Bajić D, Osborne ML, Poyatos JF, Sanchez A. 2019. High-order interactions distort the functional landscape of microbial consortia. PLoS Biol 17:e3000550.
Lawlor LR. 1979. Direct and indirect effects of n-species competition. Oecologia 43:355–364.
Dungan M. 1986. Three-way interactions: barnacles, limpets, and algae in a sonoran desert rocky intertidal zone. Am Nat 127.
Losey JE, Denno RF. 1998. Positive predator-predator interactions: enhanced predation rates and synergistic suppression of aphid populations. Ecology 79:2143–2152.
Fauth JE. 1990. Interactive effects of predators and early larval dynamics of the treefrog hyla chrysocelis. Ecology 71:1609–1616.
Worthen WB, Moore JL. 1991. Higher-order interactions and indirect effects: a resolution using laboratory drosophila communities. Am Nat 138:1092–1104.
Case TJ, Bender EA. 1981. Testing for higher order interactions. Am Nat 118:920–929.
Qin J, Li R, Raes J, Arumugam M, Burgdorf KS, Manichanh C, Nielsen T, Pons N, Levenez F, Yamada T, Mende DR, Li J, Xu J, Li S, Li D, Cao J, Wang B, Liang H, Zheng H, Xie Y, Tap J, Lepage P, Bertalan M, Batto JM, Hansen T, Le Paslier D, Linneberg A, Nielsen HB, Pelletier E, Renault P, Sicheritz-Ponten T, Turner K, Zhu H, Yu C, Li S, Jian M, Zhou Y, Li Y, Zhang X, Li S, Qin N, Yang H, Wang J, Brunak S, Doré J, Guarner F, Kristiansen K, Pedersen O, Parkhill J, Weissenbach J, MetaHIT Consortium, et al. 2010. A human gut microbial gene catalogue established by metagenomic sequencing. Nature 464:59–65.
Zhang Z, Geng J, Tang X, Fan H, Xu J, Wen X, Ma ZS, Shi P. 2014. Spatial heterogeneity and co-occurrence patterns of human mucosal-associated intestinal microbiota. ISME J 8:881–893.
May RM. 1972. Will a large complex system be stable? Nature 238:413–414.
Goldford JE, Lu N, Bajić D, Estrela S, Tikhonov M, Sanchez-Gorostiaga A, Segré D, Mehta P, Sanchez A. 2018. Emergent simplicity in microbial community assembly. Science 361:469–474.
D’Andrea R, Gibbs T, O’Dwyer JP. 2019. Emergent neutrality in consumer-resource dynamics. bioRxiv.
Shaw GTW, Pao YY, Wang D. 2016. MetaMIS: a metagenomic microbial interaction simulator based on microbial community profiles. BMC Bioinformatics 17:488.
Rolig AS, Parthasarathy R, Burns AR, Bohannan BJ, Guillemin K. 2015. Individual members of the microbiota disproportionately modulate host innate immune responses. Cell Host Microbe 18:613–620.
Wiles TJ, Schlomann BH, Wall ES, Betancourt R, Parthasarathy R, Guillemin K. 2020. Swimming motility of a gut bacterial symbiont promotes resistance to intestinal expulsion and enhances inflammation. PLoS Biol 18:e3000661.
Westerfield M, Westernfield M. 2007. The zebrafish book: a guide for the laboratory use of zebrafish (Danio rerio), 5th ed. University of Oregon, Eugene, OR.
Wiles TJ, Wall ES, Schlomann BH, Hay EA, Parthasarathy R, Guillemin K. 2018. Modernized tools for streamlined genetic manipulation and comparative study of wild and diverse proteobacterial lineages. mBio 9:e01877-18.
Stephens WZ, Burns AR, Stagaman K, Wong S, Rawls JF, Guillemin K, Bohannan BJM. 2016. The composition of the zebrafish intestinal microbial community varies across development. ISME J 10:644–654.

Information & Contributors


Published In

cover image mBio
Volume 11Number 527 October 2020
eLocator: 10.1128/mbio.01667-20
Editors: William Ludington, Carnegie Institution for Science, Edward G. Ruby, University of Hawaii at Manoa


Received: 22 June 2020
Accepted: 10 September 2020
Published online: 13 October 2020


  1. gut microbiome
  2. microbial interactions
  3. zebrafish
  4. gnotobiotic



Deepika Sundarraman
Department of Physics, University of Oregon, Eugene, Oregon, USA
Institute of Molecular Biology, University of Oregon, Eugene, Oregon, USA
Materials Science Institute, University of Oregon, Eugene, Oregon, USA
Edouard A. Hay
Department of Physics, University of Oregon, Eugene, Oregon, USA
Institute of Molecular Biology, University of Oregon, Eugene, Oregon, USA
Materials Science Institute, University of Oregon, Eugene, Oregon, USA
Dylan M. Martins
Department of Physics, University of Oregon, Eugene, Oregon, USA
Institute of Molecular Biology, University of Oregon, Eugene, Oregon, USA
Materials Science Institute, University of Oregon, Eugene, Oregon, USA
Drew S. Shields
Department of Physics, University of Oregon, Eugene, Oregon, USA
Institute of Molecular Biology, University of Oregon, Eugene, Oregon, USA
Materials Science Institute, University of Oregon, Eugene, Oregon, USA
Noah L. Pettinari
Department of Physics, University of Oregon, Eugene, Oregon, USA
Institute of Molecular Biology, University of Oregon, Eugene, Oregon, USA
Materials Science Institute, University of Oregon, Eugene, Oregon, USA
Department of Physics, University of Oregon, Eugene, Oregon, USA
Institute of Molecular Biology, University of Oregon, Eugene, Oregon, USA
Materials Science Institute, University of Oregon, Eugene, Oregon, USA


William Ludington
Invited Editor
Carnegie Institution for Science
Edward G. Ruby
University of Hawaii at Manoa


Address correspondence to Raghuveer Parthasarathy, [email protected].
Deepika Sundarraman, Edouard A. Hay, and Dylan M. Martins contributed equally to this work. Author order was determined based on contributions to data analysis and manuscript writing.

Metrics & Citations


Note: There is a 3- to 4-day delay in article usage, so article usage will not appear immediately after publication.

Citation counts come from the Crossref Cited by service.


If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

View Options

Figures and Media






Share the article link

Share with email

Email a colleague

Share on social media

American Society for Microbiology ("ASM") is committed to maintaining your confidence and trust with respect to the information we collect from you on websites owned and operated by ASM ("ASM Web Sites") and other sources. This Privacy Policy sets forth the information we collect about you, how we use this information and the choices you have about how we use such information.
FIND OUT MORE about the privacy policy