Article21 June 2018Open Access Transparent process Deciphering microbial interactions in synthetic human gut microbiome communities Ophelia S Venturelli Corresponding Author Ophelia S Venturelli [email protected] orcid.org/0000-0003-2200-1963 Department of Biochemistry, University of Wisconsin-Madison, Madison, WI, USA Search for more papers by this author Alex V Carr Alex V Carr Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Garth Fisher Garth Fisher Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Ryan H Hsu Ryan H Hsu orcid.org/0000-0001-8221-224X California Institute for Quantitative Biosciences, University of California Berkeley, Berkeley, CA, USA Search for more papers by this author Rebecca Lau Rebecca Lau Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Benjamin P Bowen Benjamin P Bowen Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Susan Hromada Susan Hromada Department of Biochemistry, University of Wisconsin-Madison, Madison, WI, USA Search for more papers by this author Trent Northen Trent Northen Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Adam P Arkin Adam P Arkin Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA California Institute for Quantitative Biosciences, University of California Berkeley, Berkeley, CA, USA Department of Bioengineering, University of California Berkeley, Berkeley, CA, USA Energy Biosciences Institute, University of California Berkeley, Berkeley, CA, USA Search for more papers by this author Ophelia S Venturelli Corresponding Author Ophelia S Venturelli [email protected] orcid.org/0000-0003-2200-1963 Department of Biochemistry, University of Wisconsin-Madison, Madison, WI, USA Search for more papers by this author Alex V Carr Alex V Carr Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Garth Fisher Garth Fisher Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Ryan H Hsu Ryan H Hsu orcid.org/0000-0001-8221-224X California Institute for Quantitative Biosciences, University of California Berkeley, Berkeley, CA, USA Search for more papers by this author Rebecca Lau Rebecca Lau Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Benjamin P Bowen Benjamin P Bowen Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Susan Hromada Susan Hromada Department of Biochemistry, University of Wisconsin-Madison, Madison, WI, USA Search for more papers by this author Trent Northen Trent Northen Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Search for more papers by this author Adam P Arkin Adam P Arkin Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA California Institute for Quantitative Biosciences, University of California Berkeley, Berkeley, CA, USA Department of Bioengineering, University of California Berkeley, Berkeley, CA, USA Energy Biosciences Institute, University of California Berkeley, Berkeley, CA, USA Search for more papers by this author Author Information Ophelia S Venturelli *,1, Alex V Carr2,‡,‡, Garth Fisher2,‡, Ryan H Hsu3, Rebecca Lau2, Benjamin P Bowen2, Susan Hromada1, Trent Northen2 and Adam P Arkin2,3,4,5 1Department of Biochemistry, University of Wisconsin-Madison, Madison, WI, USA 2Environmental Genomics and Systems Biology, Lawrence Berkeley National Laboratory, Berkeley, CA, USA 3California Institute for Quantitative Biosciences, University of California Berkeley, Berkeley, CA, USA 4Department of Bioengineering, University of California Berkeley, Berkeley, CA, USA 5Energy Biosciences Institute, University of California Berkeley, Berkeley, CA, USA ‡These authors contributed equally to this work ‡Correction added on 27 June 2018 after first online publication: Alex C Carr was corrected to Alex V Carr *Corresponding author. Tel: +1 608 263 7017; E-mail: [email protected] Molecular Systems Biology (2018)14:e8157https://doi.org/10.15252/msb.20178157 See also: C Abreu et al (June 2018) PDFDownload PDF of article text and main figures. Peer ReviewDownload a summary of the editorial decision process including editorial decision letters, reviewer comments and author responses to feedback. ToolsAdd to favoritesDownload CitationsTrack CitationsPermissions ShareFacebookTwitterLinked InMendeleyWechatReddit Figures & Info Abstract The ecological forces that govern the assembly and stability of the human gut microbiota remain unresolved. We developed a generalizable model-guided framework to predict higher-dimensional consortia from time-resolved measurements of lower-order assemblages. This method was employed to decipher microbial interactions in a diverse human gut microbiome synthetic community. We show that pairwise interactions are major drivers of multi-species community dynamics, as opposed to higher-order interactions. The inferred ecological network exhibits a high proportion of negative and frequent positive interactions. Ecological drivers and responsive recipient species were discovered in the network. Our model demonstrated that a prevalent positive and negative interaction topology enables robust coexistence by implementing a negative feedback loop that balances disparities in monospecies fitness levels. We show that negative interactions could generate history-dependent responses of initial species proportions that frequently do not originate from bistability. Measurements of extracellular metabolites illuminated the metabolic capabilities of monospecies and potential molecular basis of microbial interactions. In sum, these methods defined the ecological roles of major human-associated intestinal species and illuminated design principles of microbial communities. Synopsis Analysis of microbial interactions in a synthetic human gut microbiome community shows that pairwise microbial interactions are major drivers of multi-species community dynamics. The study reveals ecological drivers, metabolite hub species and ecologically sensitive organisms in the network. A data-driven pipeline is used to construct a predictive dynamic model of a diverse anaerobic human gut microbiome community. Design principles of stable coexistence and history-dependence are elucidated. Ecological roles and metabolite profiles are analyzed for each organism. The study highlights challenges in using phylogenetic and exo-metabolomic “signals” to predict microbial interactions and community functions. Introduction Microbes have evolved in diverse microbial communities that occupy nearly every environment on Earth, spanning extreme environments such as acid mine drains and hot springs to multicellular organisms. The gut microbiome is a dense collection of microorganisms that inhabits the human gastrointestinal tract (Lozupone et al, 2012; Earle et al, 2015; Tropini et al, 2017) and performs numerous functions to impact human physiology, nutrition, behavior, and development (Ley et al, 2005; Fischbach & Sonnenburg, 2011; Foster & McVey Neufeld, 2013; Louis et al, 2014; Sharon et al, 2014; Rooks & Garrett, 2016). Functions of the gut microbiota are partitioned among genetically distinct populations that interact to perform complex chemical transformations and exhibit emergent properties such as colonization resistance at the community level. Such collective functions are realized by the combined interactions of diverse microbial species operating on multiple time and spatial scales and could not be achieved by a single monospecies population. The degree of spatial structuring in the gut microbiota varies across length scales: At a macroscale of hundreds of micrometers, bacteria cluster into distinct habitats, whereas at a scale of micrometers, intermixing of community members has been observed (Donaldson et al, 2015; Earle et al, 2015; Mark Welch et al, 2017). The gut microbiota is composed of hundreds of bacterial species, the majority of which span the Firmicutes, Bacteroidetes, and Actinobacteria phyla (Ley et al, 2006). Constituent strains of the gut microbiota have been shown to persist in an individual over long periods of time, demonstrating that the gut microbiota exhibits stability over time (Faith et al, 2013). Perturbations to the system such as dietary shifts or antibiotic administration can shift the operating point of the gut microbiota to an alternative state (Relman, 2012). While the identities of the organisms and microbial co-occurrence relationships across individuals have been elucidated (Faust et al, 2012), we lack a quantitative understanding of how microbial interactions shape community assembly, stability, and response to perturbations. For example, the ecological and molecular forces that enable stable coexistence of the dominant phyla Firmicutes, Bacteroidetes and Actinobacteria are not well understood (Fischbach & Sonnenburg, 2011). The resilience of microbiomes, defined as the capacity to recover from perturbations, is strongly linked to microbial diversity. Indeed, a reduction in microbial diversity of the human gut microbiome is associated with multiple diseases, suggesting that a high-dimensional and functionally heterogeneous ecosystem promotes human health (Sommer et al, 2017). Understanding the molecular and ecological factors influencing the stability and resilience of the gut microbiota has implications for the development of targeted interventions to modulate microbiome states. Central to this problem is inferring unknown microbial interactions and developing tools to predict temporal changes in community behaviors in response to environmental stimuli. Cooperation and competition generate positive and negative feedbacks in microbial communities and influence functional activities and stability. Negative interactions have been shown to dominate microbial inter-relationships in synthetic aquatic microcosms (Foster & Bell, 2012). However, the prevalence of competition and cooperation in microbial communities occupying other diverse environments such as the human gut microbiota remains elusive. Direct negative interactions in microbial consortia can originate from competition for resources or space, biomolecular warfare, or production of toxic waste products (Hibbing et al, 2010). Positive interactions can stem from secreted metabolites that are utilized by a community member or detoxification of the environment. Pairwise microbial interactions can be modified by a third organism, leading to higher-order effects that influence community behaviors (Bairey et al, 2016). Ecological driver species, which exhibit a large impact on community structure and function, represent key nodes in the network that could be manipulated to control community states (Gibson et al, 2016). Predicting community dynamics is a key step toward understanding the organizational principles of microbial communities. Computational models at different resolutions can be used to analyze and predict the behaviors of microbial communities (Faust & Raes, 2012). Dynamic computational models can be used to investigate temporal changes in community structure, and tools from dynamical systems theory can be used to analyze system properties including stability and parameter sensitivity (Astrom & Murray, 2010). Generalized Lotka–Volterra (gLV) is an ordinary differential equation model that represents microbial communities with a limited number of parameters that can be deduced from time-series data. Here, we develop a systematic modeling and experimental pipeline to construct a predictive computational model of microbial community dynamics and interrogate microbial interactions mediating community assembly. Time-resolved measurements of monospecies and pairwise assemblages were used to train a dynamic computational model of a diverse synthetic human gut microbiome community. Our model revealed a high proportion of negative and frequent positive interactions. Specific ecological driver species and responsive organisms to community context were identified in the network. Bacteroidetes exhibited an overall negative impact on the community, whereas specific members of Actinobacteria and Firmicutes displayed numerous positive outgoing interactions. A prevalent pairwise sub-network composed of positive and negative interactions exhibited robust species coexistence to variations in model parameters. Our model showed that the majority of history-dependent responses in pairwise consortia were due to slow convergence to a steady state composition and these networks were enriched for negative interactions. The metabolic capabilities of monospecies were elucidated using exo-metabolomics profiling, and these data pinpointed a set of metabolites predicted to mediate positive and negative interactions. However, the metabolite profiles failed to forecast specific influential organisms modulating community assembly. Together, these results show that combinations of pairwise interactions can represent the assembly of multi-species communities and such pairwise couplings can realize a diverse repertoire of dynamic behaviors. Results Probing the temporal behaviors of monospecies and pairwise assemblages We aimed to dissect the microbial interactions influencing community assembly in a reduced complexity model gut community spanning the major phyla Bacteroidetes, Firmicutes, Actinobacteria, and Proteobacteria. To this end, a synthetic ecology encompassing prevalent human-associated intestinal species Bacteroides thetaiotaomicron (BT), Bacteroides ovatus (BO), Bacteroides uniformis (BU), Bacteroides vulgatus (BV), Blautia hydrogenotrophica (BH), Collinsella aerofaciens (CA), Clostridium hiranonis (CH), Desulfovibrio piger (DP), Eggerthella lenta (EL), Eubacterium rectale (ER), Faecalibacterium prausnitzii (FP), and Prevotella copri (PC) was designed to mirror the functional and phylogenetic diversity of the natural system (Fig 1A; Qin et al, 2010). These species have been shown to contribute significantly to human health and are implicated in multiple human diseases (Watterlot et al, 2008; Larsen et al, 2010; Thota et al, 2011; Fujimoto et al, 2013; Haiser et al, 2013; Scher et al, 2013; Table 1). Figure 1. Experimental design for high-throughput characterization of synthetic human gut microbiome consortia Phylogenetic tree of the 12-member synthetic ecology spanning the major phyla in the gut microbiome including Actinobacteria, Bacteroidetes, Firmicutes, and Proteobacteria. Phylogenetic analysis was performed using a concatenated alignment of single-copy marker genes obtained via PhyloSift (preprint: Darling et al, 2014). Maximum likelihood trees were generated using default options. The scale bar represents the number of substitutions per site in the alignment. Schematic of the experimental design for this study. Species were combined using an approximately 1:1 or 19:1 initial proportion based on absorbance measurements at 600 nm (OD600) into microtiter plates using liquid-handling robotic manipulation. Approximately every 12 h, samples were collected for multiplexed 16S rRNA gene sequencing (black circles). Relative abundance was measured using multiplexed 16S rRNA gene sequencing of the V3–V4 region using dual-indexed primers compatible with an Illumina platform (stacked bar plot, bottom right). Serial transfers were performed in approximately 24 intervals (purple bars, top) by transferring an aliquot of the communities into fresh media using a 1:20 dilution. In parallel, time-resolved OD600 measurements of monospecies and consortia were performed. Download figure Download PowerPoint Table 1. Table of species used in study and associations with human diseases based on previous literature. Arrows pointing up or down denote positive or negative associations, respectively Species Association(s) Prevotella copri (PC) Inflammatory and autoimmune disease (↑) (Scher et al, 2013), autism (↓) (Kang et al, 2013) Bacteroides vulgatus (BV) Ulcerative colitis (↑) (Bamba et al, 1995) Bacteroides uniformis (BU) Metabolic/immunological dysfunction (↓) (Gauffin Cano et al, 2012) Bacteroides ovatus (BO) Type I diabetes (↑) (Giongo et al, 2011) Bacteroides thetaiotaomicron (BT) Ulcerative colitis (↑)(Bloom et al, 2011) Faecalibacterium prausnitzii (FP) Crohn's disease (↓) (Watterlot et al, 2008), inflammatory bowel disease (↓) (Segain et al, 2000), Celiac disease (↓) (De Palma et al, 2010) Blautia hydrogenotrophica (BH) Healthy human colon (↑) (Nava et al, 2012) Eubacterium rectale (ER) Type II diabetes (↑) (Larsen et al, 2010) Collinsella aerofaciens (CA) Colon cancer (↓) (Moore & Moore, 1995), rheumatoid arthritis (↑) (Chen et al, 2016) Eggerthella lenta (EL) Cardiac drug transformations (↑) (Haiser et al, 2013), Crohn's disease (↑) (Thota et al, 2011), rheumatoid arthritis (↑) (Chen et al, 2016) Desulfovibrio piger (DP) Regressive autism (↑) (Finegold et al, 2012) Clostridium hiranonis (CH) None reported Synthetic assemblages were arrayed in microtiter plates in an anaerobic chamber using an automated liquid-handling procedure (see Materials and Methods). A rich media (see Materials and Methods) was selected to support the growth of all monospecies. The communities were serially transferred at 24-h intervals to prevent strains with long lag phases from being eliminated and allow communities to approach a steady state composition by monitoring assembly over many cell generations. Further, serial transfers can also reflect recurrent temporal perturbations to the gut microbiota such as diet and colonic transit time (Fig 1B). Multiplexed 16S rRNA gene sequencing was performed in approximately 12-h intervals to elucidate the temporal variations in community structure at different community growth stages. The relative abundance for each species was computed as the sum of the read counts for each organism divided by the total number of reads per condition (see Materials and Methods). Since model construction is aided by absolute abundance information (Bucci et al, 2016; Widder et al, 2016), the total biomass of the communities was monitored approximately every 30 min using absorbance at 600 nm (OD600). Cellular traits such as cell adhesion, size, and shape can influence OD600 measurements (Stevenson et al, 2016). In addition, counting of colony-forming units (CFU) is biased by cell adhesion, dormant sub-populations, growth selection on solid vs. liquid media, and growth stage (Jansson & Prosser, 1997; Volkmer & Heinemann, 2011; Ou et al, 2017). Our model was trained on absolute abundance estimated from OD600 measurements and used to predict absolute abundance based on OD600 and thus automatically accounts for any potential biases. To infer microbial interactions, time-resolved measurements of all monospecies and pairwise communities (66 combinations) were performed using an approximately 1:1 initial abundance ratio based on OD600 values (PW1 dataset, Appendix Fig S1, Dataset EV1). Monospecies growth and community composition were measured using OD600 measurements of total biomass and multiplexed 16S rRNA gene sequencing, respectively (see Materials and Methods). The monospecies displayed a broad range of growth rates, carrying capacities and lag phases (M dataset, Appendix Fig S2). Pairwise consortia exhibited diverse growth responses and dynamic behaviors including coexistence and single-species dominance (Appendix Fig S1). The distribution of absolute abundance of each species across communities in PW1 provided insight into variability in growth in the presence of a second organism (Appendix Fig S3A). Absolute species abundance was normalized to the monospecies maximum OD600 value to evaluate relative changes in the baseline fitness of each organism in the presence of second species. Bacteroides and CH displayed the lowest coefficient of variation (CV), indicating that the fitness levels of these organisms were not significantly modified by a second species (Appendix Fig S3B). The remaining species displayed a bimodal (FP, DP, PC, BH, and CA), long-tail distribution (ER and EL), and/or high CV (ER, BH, CA, and PC), demonstrating that growth was significantly altered in the presence of specific organisms. To further probe the dynamic responses of pairwise consortia, a set of 15 consortia (Fig 2B, Dataset EV2) inoculated at different initial species proportions based on OD600 values (95% species A, 5% species B, and the second wherein these percentages were reversed) were characterized using our experimental workflow (PW2 dataset, Appendix Fig S4). The community behaviors were classified into the following categories based on a quantitative threshold in species proportions at 72 h: (i) single-species dominance; (ii) stable coexistence wherein both species persisted above an abundance threshold; (iii) history dependence whereby communities inoculated using distinct initial species proportions mapped to different community structures; or (iv) other for communities that did not quantitatively satisfy the relative abundance thresholds for cases 1–3. A subset of the communities classified in the other category displayed weak history-dependent responses potentially attributed to variations in biological replicates. The qualitative behaviors of the remaining 51 pairwise communities were classified based on community structure at an initial (t = 0) and final (t = 72 h) time point using the PW2 experimental design wherein the organisms were inoculated using different initial species proportions (95% species A, 5% species B, and the reciprocal percentages, Appendix Fig S5A). Together, these results demonstrated that approximately 50, 24, and 12% of pairwise communities displayed dominance, stable coexistence, and history dependence, respectively (Appendix Fig S5B). Figure 2. Model training of generalized Lotka–Volterra (gLV) to time-resolved measurements of monospecies and pairwise assemblages Species relative abundance as a function of time for all pairwise communities. Experimental measurements and model fits based on T3 are represented as data points and lines, respectively. In each subplot, time and species relative abundance are displayed on the x- and y-axis, respectively. Stars denote datasets with a sum of mean squared errors greater than 0.15. Error bars represent 1 s.d. from the mean of at least three biological replicates. Temporal changes in species relative abundance of a selected set of pairwise assemblages inoculated at 5% species A, 95% species B or 95% species A, 5% species B based on OD600 values. Time and relative abundance are represented on the x- and y-axis, respectively. Data points and lines represent experimental measurements and model fits to T3, respectively. Error bars represent 1 s.d. from the mean of at least three biological replicates. Stars denote datasets with a sum of mean squared errors greater than 0.15. Inferred inter-species interaction coefficients for the gLV model trained on T3. Gray and green edges denote negative (αij < 0) and positive (αij > 0) interaction coefficients. The edge width and node size represent the magnitude of the inter-species interaction coefficient and steady state monospecies abundance (xe = −μiαii−1), respectively. To highlight significant interactions, inter-species interaction coefficients with a magnitude less than 1e-5 were not displayed. Download figure Download PowerPoint Construction of a dynamic computational model of the community A generalizable modeling framework was developed to infer parameters from time-series measurements of relative abundance and total biomass (OD600). The generalized Lotka–Volterra (gLV) model represents microbial growth, intra-species interactions, and pairwise inter-species interactions and can be used to predict the dynamic behaviors of the community and analyze system properties such as stability and parameter sensitivity. The model equations are given by: where n, μ, αii, and αij represent the number of species, growth rates, intra-species, and inter-species interaction coefficients, respectively. To minimize overfitting of the data, a regularized parameter estimation method was implemented that penalized the magnitude of the parameter values (see Materials and Methods). Three training sets were evaluated based on predictive capability: (T1) M; (T2) M, PW1; and (T3) M, PW1, PW2. A range of regularization coefficient values (λ) was scanned to balance the goodness of fit to the training sets and degree of sparsity of the model (Appendix Fig S6). The parameterized gLV model trained on T3 captured the majority of pairwise community temporal responses (Fig 2A and B). However, the model did not accurately represent the dynamic behaviors of a set of communities including BH, EL; PC, CA; BO, CH; ER, BH; and PC, BH based on a threshold in the mean squared error between the model and data. Thresholding the magnitude of the inter-species interaction coefficients using a value of 1e-5 yielded a densely connected network whereby 77% of species pairs exhibited an interaction. The network connectivity varied between 75 and 79% for interaction coefficient thresholds ranging from 1e-6 to 1e-3. Interaction coefficients with a magnitude less than 1e-3 are not expected to change the steady state species abundance based on the inferred growth rate and interaction coefficient values. Of these interactions, 56 and 21% were negative and positive, respectively (Fig 2C). Negative interactions can arise from resource competition, biomolecular warfare, or production of toxic waste by-products. Positive interactions can originate from metabolite secretion or detoxification of the environment. Bacteroides (BO, BV, BU, and BT) displayed a net negative impact on the network, whereas EL, BH, and CH positively stimulated a large number of species (Appendix Fig S7A). Pairwise networks were enriched for unidirectional negative (−/0, 36%), bidirectional negative (−/−, 32%), and positive and negative (+/−, 26%) species couplings (Appendix Fig S7B). The contribution of each species to full community assembly in the model is dictated by a set of coupled ordinary differential equations that are a function of the monospecies growth rates, intra-species interactions, and outgoing and incoming inter-species interactions (see Materials and Methods). Therefore, a prediction about the role of each organism in community assembly requires simulation and analysis of the gLV model. FP was the recipient of five positive interactions, suggesting that the fitness of FP is coupled to the composition of the community (Appendix Fig S7A). To determine the contribution of each incoming positive interaction on FP abundance, we examined a 6-member gLV model composed of FP, BH, BU, BV, CH, and DP. The combined set of five positive inter-species interactions was required to alter FP abundance by more than twofold, and single and dual inter-species interactions moderately increased FP abundance at 72 h (Appendix Fig S7C). Therefore, FP represents an ecologically responsive organism that is significantly enhanced by the presence of multiple organisms in the community in these conditions. Corroborating this notion, FP exhibited significant variability in absolute abundance across PW1 communities and frequent coexistence with other organisms (Appendix Figs S3 and S5B). Strong positive or negative interactions can be deciphered by an enhancement or reduction in community productivity compared to a null model representing the sum of monospecies productivities (Foster & Bell, 2012). We computed the integral of biomass (OD600) over time to evaluate monospecies and pairwise community productivities (Appendix Fig S8). Our results showed that the productivities of 38% of pairwise consortia were less than twofold compared to the predictions based on the null models, consistent with the prevalence of negative interactions. Bacteroides pairwise communities BV, BU; BV, BO; BT, BU; BU, BO; BT, BV; and BT, BO exhibited significantly lower pairwise productivities compared to the null model, which was consistent with the inferred mutual inhibitory network topologies (Fig 2C). The productivities of EL, BH; EL, BU; EL, BT; EL, BO; EL, BV; and CH, ER were significantly enhanced in comparison with the predictions based on the null model. In the inferred network, five of these consortia (EL, BU; EL, BT; EL, BO; EL, BV; and CH, ER) displayed coupled negative and positive interaction topologies and EL, BH exhibited mutualism, demonstrating that both unidirectional and bidirectional positive interactions can augment community productivity. Hierarchical clustering of the gLV interaction coefficients showed that members of Bacteroides or Actinobacteria exhibited similar patterns in outgoing and incoming microbial interactions (Appendix Fig S9A and B). However, the Firmicutes clustering pattern did not reflect the phylogenetic relationships. Together, these data show that distantly related species can display similar microbial interactions (e.g., BH and DP or PC and CA) and closely related species can exhi