Stochastic block models reveal a robust nested pattern in healthy human gut microbiomes — TEI-tagged view

Sergio Cobo-LópezpersName, Vinod K GuptapersName, Jaeyun SungpersName, Roger GuimeràpersName, Marta Sales-PardopersName (2022date). “Stochastic block models reveal a robust nested pattern in healthy human gut microbiomes.” PNAS Nexustitle·j.

Object of study — robust structural patterns in the taxonomic composition of healthy human gut microbiomesterm·key

Research question — what robust structural patterns underlie gut microbiometerm·key taxonomic composition

Contribution — mixed-membershipterm·key enterotypesterm·dom via stochastic block modelsterm·mth, revealing a nested generalist-specialist structure

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Stochastic block modelsterm·mth reveal a robust nested pattern in healthy human gut microbiomesterm·key Abstract A key question in human gut microbiometerm·key research is what are the robust structural patterns underlying its taxonomic composition. Herein, we use whole metagenomic datasets from healthy human guts to show that such robust patterns do exist, albeit not in the conventional enterotypeterm·dom sense. We first introduce the concept of mixed-membershipterm·key enterotypesterm·dom using a network inference approach based on stochastic block modelsterm·mth. We find that gut microbiomesterm·key across a group of people (hosts) display a nested structure, which has been observed in a number of ecological systems. This finding led us to designate distinct ecological roles to both microbes and hosts: generaliststerm·key and specialiststerm·key. Specifically, generalist hoststerm·dom have microbiomes with most microbial species, while specialist hoststerm·dom only have generalistterm·key microbes. Moreover, specialistterm·key microbes are only present in generalist hoststerm·dom. From the nested structure of microbial taxonomies, we show that these ecological roles of microbes are generally conserved across datasets. Our results show that the taxonomic composition of healthy human gut microbiomesterm·key is associated with robustly structured combinations of generalistterm·key and specialistterm·key species.

Is the gut microbiometerm·key a random ecosystem or does it have an internal structure? We show that gut microbial communities in healthy individuals have a well-defined and robust ecological order consisting of gradually increasing degrees of specialization: while some microbes are found in the majority of humans, others are very specific to some people. Similarly, hosts have varying degrees of microbial diversity. Our findings describe the principles of how healthy human gut microbiomesterm·key are associated with robustly structured combinations of generalistterm·key and specialistterm·key species.

Introduction Among the different communities of microbes that live in symbiosis inside and on the human body, the gut microbiometerm·key has received the most attention by far ( 1–3 ). During the last decade, gut microbiometerm·key dysbiosis was found to be linked to a variety of diseases related to the gut ( 1 , 4–6 ), and even to those not overtly gut-related, including fertility disruption ( 7 ), neurological pathologies ( 8 ), or atopy and asthma in infants ( 9–11 ). Recent studies continue to show that the gut microbiometerm·key plays an important role in health and disease, such as influencing immune response in infants ( 11 ) and regulating inflammation ( 12 ).

The need to better understand why the gut microbiometerm·key is such an important component in human health has sent scientists on a quest to discover “healthy” gut microbiometerm·key compositions. Such an understanding is critical for the design of better therapies based on fecal microbiota transplantation (FMT) ( 13 , 14 ) or gut microbiome-targeting drug treatments ( 3 ). Despite recent progress, the scientific community has yet to reach a consensus that pinpoints the essential microbial components for our well-being. By and large, this has been due to the large variability observed in gut microbiometerm·key composition and that only few species are present in the majority of samples ( 15–18 ). This can be partially explained by an ecological perspective, in which the occupation of niches determines the governing functions of the ecological system ( 15 ). From birth, the microbial niches of the gastrointestinal tract are affected by many dietary and environmental factors, resulting in much of the observed variability ( 2 , 19 , 20 ).

Despite their known heterogeneity in taxonomic composition from person to person, gut microbiomesterm·key from healthy subjects are able to perform nearly the same critical biological functions ( 21–23 ); this suggests that there must be recurrent patterns in the structure of the human gut microbiometerm·key that are linked to the healthy phenotype. An approach to characterize and classify such patterns is the identification of so-called enterotypesterm·dom, defined by specific traits of microbiome composition profiles that are common to particular groups of hosts. However, despite some consistencies in the most abundant species, which suggests that there are taxa more preferred than others, enterotypesterm·dom have not been found to be convincingly robust across, for example, geographically distinct populations ( 23 ). As a result, the search for robust, rigid classification criteria of microbiome compositions (such as the original concept of enterotypesterm·dom) has been somewhat abandoned in favor of more nuanced studies aiming to understand the taxonomic and functional landscape of the gut microbiometerm·key and its connection to various pathologies ( 23 ).

The fundamental assumption behind the definition of enterotypesterm·dom is that there are groups of hosts that have similar microbe abundance profiles. On the other hand, there is a class of generative models called bipartite stochastic block modelsterm·mth (SBMsterm·mth) ( 24 ), which (for microbiome research purposes) explicitly assume the existence of groups of hosts and groups of microbes that are associated according to similarities in the microbe abundance profiles. These statistical models can be seen, in practice, as a biclustering techniqueterm·mth that allows simultaneous clustering of both microbes and hosts. Additionally, they are easily interpretable and are amenable to rigorous and unbiased model selection tools to identify optimal groupings of hosts and microbes ( 25 ). In fact, SBMsterm·mth have been successfully used in the complex networks literature to characterize the large-scale structure of complex networks and to make predictions of unobserved data ( 26–28 ).

In our study, we represent gut microbe abundance profiles from a cohort of healthy people as a network and follow a Bayesian SBMterm·mth approach to show: first, that the SBMsterm·mth are able to predict unobserved microbe abundances better than agglomerative hierarchical clustering methods; and second, that there exists a robust organization structure in healthy gut microbiomesterm·key based on taxonomic composition. Our analysis reveals that microbe abundances display the same nested structure that was first observed in faunas of fragmented habitats and archipelagos ( 29 ), and later discovered in a wide range of ecological systems ( 30–33 ).

In this nested structure of gut microbiomesterm·key, we find that a few groups of hosts have significant abundances of a large number of microbial groups (generalist hoststerm·dom), while other groups of hosts have very few microbial groups of significant abundance (specialist hoststerm·dom). Analogously, some microbial groups are present in nearly all groups of hosts (generalistterm·key microbial groups), while other microbial groups are present in just a few groups of hosts (specialistterm·key microbial groups).

Our results suggest that the taxonomic composition of healthy human gut microbiomesterm·key follow a nested structure, which are not only similar to those found in other ecological networks throughout nature, but also predictive of unobserved abundances. By showing that there are generalist hoststerm·dom as well as specialist hoststerm·dom, and that there is a varying gradient of such ecological roles across different hosts (in contrast to discrete enterotypesterm·dom), we mathematically formulate a new description of how gut microbiomesterm·key organize within a small population.

Materials and Methods Downloading and quality control of sequenced reads The current investigation uses publicly available human gut microbiometerm·key datasets from five published studies (Table 1 ). Of note, for studies whose cohort includes microbiome samples from both cases (disease) and controls (healthy), only samples from the controls were considered. Raw sequence files (.fastq) were downloaded from the NCBI Sequence Read Archive (SRA) and European Nucleotide Archive (ENA) databases. Sequence reads were processed with the KneadData quality-control pipeline, which uses Trimmomatic v0.36 and Bowtie2 v0.1 for removal of low-quality read bases and human reads, respectively. Trimmomatic v0.36 was run with parameters SLIDINGWINDOW:4:30, and the Phred quality score threshold was set at “<30”. Illumina adapter sequences were removed, and trimmed nonhuman reads shorter than 60 bp in nucleotide length were discarded. Potential human contamination was filtered by removing reads that aligned to the human genome (reference genome hg19). Furthermore, stool metagenome samples of low read count after quality filtration (< 1 M reads) were excluded from our analysis.

Description of healthy human gut microbiometerm·key datasets used in this study.

Species-level taxonomic profiling Taxonomic profiling was done using the MetaPhlAn2 v2.7.0 phylogenetic clade identification pipeline ( 34 ) using default parameters. Briefly, MetaPhlAn2 classifies metagenomic reads to taxonomies based on a database of clade-specific marker genes derived from ≈17,000 microbial genomes (corresponding to ≈13,500 bacterial and archaeal, ≈3,500 viral, and ≈110 eukaryotic species). After taxonomic profiling, the following stool metagenome samples were discarded from our analysis: (i) samples composed of more than 5% unclassified taxonomies; and (ii) phenotypic outliers according to a dissimilarity measure. More specifically, Bray–Curtis distances were calculated between each sample of a particular phenotype and a hypothetical sample in which the species’ abundances were taken from the medians across samples. A sample was considered as an outlier, and thereby removed from further analysis, when its dissimilarity exceeded the upper and inner fence (i.e., > 1.5 times outside of the interquartile range above the upper quartile and below the lower quartile) of all dissimilarities. Finally, species of viral origin, those of either unclassified or unknown clades, and those of low prevalence (i.e., observed in < 1% of the samples in each study), were excluded from our study.

Discretization of microbiome datasets Each gut microbiometerm·key dataset of a study can be represented as a species matrix M , wherein columns represent gut microbiometerm·key samples of the hosts and rows correspond to microbial species. Each element M hm corresponds to the relative abundance (i.e., proportion) of microbial species m in the sample of host h . Importantly, because our interest was to find the leading patterns in microbiome composition, we only considered as nonzero the relative abundances that were larger than 10 −4 . In addition to matrices of species-level abundances, we considered abundance matrices of phylum, class, order, family, and genus taxonomic ranks.

We model species abundance matrices as bipartite multilink networks, i.e. networks with two classes of nodes (hosts and microbes) in which edges run between nodes of different classes and edges are of different types. Here, edge weights correspond to relative abundances of microbes in each host. Since relative abundances are continuous and our models require a discrete number of edge types, we need to discretize the relative abundances and then define a new species matrix \(\hat{\bf M}\) in which entries \(\widehat{M}_{hm}\) are categorical. We choose to discretize the microbe abundances according to their order of magnitude. In all five datasets, abundances ranged from 10 −4 to 10 −1 , and so we considered three types of edges (categorical abundances): null or negligible, \(\widehat{M}_{hm}=\rm I\) if M mh < 10 −4 ; low, \(\widehat{M}_{hm}=\rm II\) if 10 −4 ≤ M mh ≤ 10 −3 ); and high, \(\widehat{M}_{hm}= \rm III\) if M mh > 10 −3 . Despite the known limitations associated with discretizing data, representing the abundance data by a small set of classes provides a simpler and more straightforward interpretation of “microbial abundance” without significantly altering the inherent distribution of the compositional data.

SBMsterm·mth for microbe abundance matrix modeling and enterotypeterm·dom identification We define a bipartite network with categorical edges whose weights are defined by the microbial relative abundances. This bipartite network is the starting point toward our overall aim, which is to find underlying groups of both microbes and hosts (essentially two classes of nodes) that together best explain (data-grounded) host–microbiome relationships.

Our bipartite SBMsterm·mth assume that hosts can be classified into groups that connect similarly to the same microbes, that is, that have the same relative abundances of the same microbes. At the same time, we assume that microbes can also be classified into groups that have similar relative abundances in similar groups of hosts. Furthermore, to each host and microbe group pairs ( h, m ), we assign a set of probabilities \(p_{hm}=\lbrace p_{hm}^{\rm I}, p_{hm}^{\rm II}, p_{hm}^{\rm III}\rbrace\) , which correspond to the probability that a host in group h has negligible, low, or high relative abundances of microbes in group m . Note that in conventional enterotypeterm·dom identification approaches, each group of hosts would define a so-called enterotypeterm·dom.

If there are L groups of microbes, then each host group h is assigned a vector of sets of probabilities for all the L microbe groups, p h = ( p h 1 , p h 2 ,..., p hM ). We call p h enterotypeterm·dom, since the vector of probability sets determines the most likely set of microbe abundances for those hosts in group h in the same way that conventional enterotypesterm·dom do (Fig. 1 ).

SBMsterm·mth provide a theoretical framework to model microbe abundance patterns. SBMsterm·mth assume the existence of groups of hosts and groups of microbes. Here, we illustrate the main concept of SBMsterm·mth in the context of microbe abundance matrix modeling and mixed-membershipterm·key enterotypeterm·dom identification. (a) Probabilistic matrix of abundances. Each row corresponds to a group of hosts and each column to a group of microbes. Hosts in the same group display a similar microbe abundance pattern (row); microbes in the same group display similar relative abundance patterns across groups of hosts (column). For simplicity, we illustrate the case in which only two possible relative abundances are allowed: non-negligible (relative abundance ≥ 10 −4 ) and negligible (relative abundance < 10 −4 ). Colors represent the probability that a group of microbes is present in a group of hosts with non-negligible relative abundances (dark blue, p (non-negligible) → 1; light blue, p (non-negligible) → 0). (b) Relationship between SBMsterm·mth and enterotypesterm·dom. In single-membership SBMsterm·mth, each host belongs to a host group exclusively (hot pink in the example), so that the corresponding (hot pink) row of the probability abundance matrix (in (a)) indicates the most likely abundances of each microbe group in hosts of that host group. Therefore, single-membership SBMsterm·mth are equivalent to conventional enterotypesterm·dom. In mixed-membershipterm·key SBMsterm·mth, each host belongs to a mixture of host groups with a certain probability given by her membership vector. The enterotypeterm·dom of each host is then a weighted combination of underlying enterotypesterm·dom in the probabilistic matrix of abundances in (a) as depicted in the weighted sum of matrix rows. The latter model is more expressive and is better suited to model complex data ( 28 , 35 ).

Mixed-membershipterm·key stochastic block modelsterm·mth for enterotypeterm·dom modeling and identification To increase the expressibility of our model, we consider a bipartite mixed-membershipterm·key stochastic block modelsterm·mth (MMSBMs) ( 28 , 35 ) in which we allow for real hosts and microbes to belong to all of the groups with a finite probability (Fig. 1 ). In this case, the groups are no longer associated to a group of hosts or microbes and, therefore, they become latent groupsterm·mth, that is underlying groups that we cannot directly measure. Each latent groupterm·mth of hosts is then characterized by its corresponding latent enterotypeterm·dom. Therefore, a possible interpretation is that the observed relative microbe abundances for hosts are the result of mixing underlying enterotypesterm·dom.

Formally, we model the multiedge bipartite network defined by categorical species matrices \(\widehat{{\bf M}}\) as a bipartite MMSBM. MMSBMs assume there exist K latent groupsterm·mth of hosts and L latent groupsterm·mth of microbes. Each host h belongs to host group k with probability θ hk and each microbe m belongs to group l with probability η ml . These probabilities are encoded in \(\boldsymbol{\theta }_h\) and \(\boldsymbol{\eta }_m\) , the so-called mixed-membershipterm·key vectors for host h and microbe m : (1) \[\begin{eqnarray*} \boldsymbol{\theta }_h=\left(\theta _{h1}, \theta _{h2}, \dots , \theta _{hK} \right) \qquad \boldsymbol{\eta }_m=\left(\eta _{m1}, \eta _{m2}, \dots , \eta _{mL} \right)~~. \end{eqnarray*}\] Because θ h and η m are probability vectors, they are subject to the normalization condition that all probabilities have to add to one, (2) \[\begin{eqnarray*} \sum ^{K}_{k=1}\theta _{hk}=1 \qquad \sum ^{L}_{l=1}\eta _{ml}=1 ~. \end{eqnarray*}\] Note that if single membership SBMsterm·mth are the best description for observed abundance matrices, then we will find that hosts have a large probability of belonging to only one group.

If \({\cal A}:=\lbrace \rm I,\rm II,\rm III\rbrace\) is the set of possible values each edge can take, MMSBMs further assume that the probability that microbes in group m have relative abundance \(a\in {\cal A}\) in hosts of group h is p ( a ) k ℓ . Since for each pair of groups ( k , ℓ) each abundance can only take one value, these probabilities are also subject to the normalization condition (3) \[\begin{eqnarray*} \sum _{a\in {\cal A}} p(a)_{k\ell } =1~~~\forall (k,\ell ). \end{eqnarray*}\] Note that p is a tensor of dimensions K × L × A , with \(A\equiv |{\cal A}|\) . Each ’row’ k in this tensor defines the latent enterotypeterm·dom for group k . If there are only two abundance types A = 2, and because of the normalization condition we only need one matrix p of dimensions K × L to define the model as we illustrate in Fig. 1 .

Given the mixed-membershipterm·key vectors \(\boldsymbol{\theta }\) and \(\boldsymbol{\eta }\) and the probability tensor p , the MMSBM defines the probability that microbe m has abundance of type a in host h as (4) \[\begin{eqnarray*} {\rm Pr}[\widehat{M}_{mh}=a] = \sum _{k, \ell } \theta _{hk}\, p_{k \ell } (a) \, \eta _{m\ell } . \end{eqnarray*}\] Inference for model parameter estimation Given a categorical species matrix \(\widehat{\bf M}^O\) of observed microbial abundances of M microbial species in P hosts, the MMSBMs assumes that each relative abundance observation is independent from the others. Therefore, the probability of observing \(\widehat{\bf M}^O\) according to the model is given by the following likelihood function, (5) \[\begin{eqnarray*} p(\widehat{\bf M}^O |{\rm MMSBM})=p(\widehat{\bf M}^O |\boldsymbol{\theta },\boldsymbol{\eta }, \mathbf {p}) =\prod _{(h,m)\in \widehat{\bf M}^O } {\rm Pr}[\widehat{M}^O_{hm}]. \end{eqnarray*}\] Note that the product runs over all observed abundances. If some abundances are missing (i.e. not observed), then these terms do not appear in the likelihood.

If we assume that we have flat priors over model parameters, then we can obtain the set of parameters \((\boldsymbol{\theta }^{*},\boldsymbol{\eta }^{*}, \mathbf {p}^{*})\) that best describe the observed data by maximizing the likelihood (6) \[\begin{eqnarray*} \{\boldsymbol{\theta }^{*},\boldsymbol{\eta }^{*}, \mathbf{p}^{*}\} = \displaystyle\mathop{{\rm arg\,\, max}}_{\{{\boldsymbol\theta },\,{\boldsymbol\eta },\,{\bf {p}}\}}\,\, \{p\,(\widehat{\bf M}^{o}|{\boldsymbol\theta },\,{\boldsymbol\eta },\,{\bf{p}})\} . \end{eqnarray*}\] Equations for the Expectation Maximization algorithm We follow a variational approach to obtain iterative equations for the model parameters and use an Expectation Maximization algorithm to obtain estimates for \(\lbrace \boldsymbol{\theta }^{*},\boldsymbol{\eta }^{*}, \mathbf {p}^{*}\rbrace\) ( 28 ). Specifically, we find for the model parameters \(\boldsymbol{\theta }, \boldsymbol{\eta }, \mathbf {p}\) iterative equations that depend on auxiliary probability distributions as it is usual in Expectation Maximization approaches.

In particular, we find the following equations ( 28 ): (7) \[\begin{eqnarray*} \theta _{hk}& = & \frac{1}{d_h} \bigg [ \sum _{a \in \lbrace \rm I,\rm II,\rm III\rbrace } \sum _{m \in \partial _h^{a}} \sum _{\ell } \omega _{hm}(k,\ell ) \bigg ] \, , \end{eqnarray*}\] (8) \[\begin{eqnarray*} \eta _{m\ell }& =& \frac{1}{d_{m}} \bigg [\sum _{a \in \lbrace \rm I,\rm II,\rm III\rbrace } \sum _{h \in \partial _m^{a}} \sum _{k} \omega _{hm}(k,\ell ) \bigg ] \, , \end{eqnarray*}\] (9) \[\begin{eqnarray*} p_{k\ell }(a)& =& \frac{\sum _{ [(h,m) \in R^o|a_{hm}=a]} \omega _{hm}(k,\ell )}{\sum _{(h^{\prime },m^{\prime }) \in R^o} \omega _{h^{\prime }m^{\prime }}(k,\ell )} \, , \end{eqnarray*}\] where \(\partial _m^{a}=\lbrace h |\widehat{M}_{mh}=a\rbrace , \, a \in \lbrace \rm I,\rm II,\rm III\rbrace\) , and \(\partial _h^{a}=\lbrace m|\widehat{M}_{mh}=a\rbrace , \, a \in \lbrace \rm I,\rm II,\rm III\rbrace\) . d h and d m are the degrees of host h and microbe m in the bipartite representation of \(\widehat{M}\) , that is, the number of microbes with measured relative abundance reported for host h , and the total number of hosts for which the relative abundance of microbe m has been reported, respectively.

In the expressions above, ω pm ( k , ℓ) is an auxiliary distribution that appears in the variational approach used to maximize likelihood functions using Expectation Maximization algorithms ( 28 ). This distribution represents the probability that host h having the microbe m with relative abundance a hm is due to h and m belonging to latent groupsterm·mth k and ℓ, that is (10) \[\begin{eqnarray*} \omega _{hm}(k,l)=\frac{\theta _{hk} \eta _{m\ell } p_{k,\ell }(a)}{\sum _{k^{\prime }\ell ^{\prime }}\theta _{hk^{\prime }} \eta _{m\ell ^{\prime }} p_{k^{\prime }\ell ^{\prime }}(a)} . \end{eqnarray*}\] Selection of K and L We assume that the best model is the most predictive one. We, therefore, look for the combination of model parameters K and L that maximizes the predictive accuracy (see section on metrics and Figure S1 in Supplementary Material ). In order to be consistent, we test our results on the combination of all five datasets. We find that the most predictive combination of parameters corresponds to K = 10, L = 20, that is 10 latent enterotypesterm·dom and 20 latent groupsterm·mth of microbes (see Figure S2, Supplementary Material ) Host-specific cross-validation experiments For each microbe abundance profile of host h , we randomized the order of abundances and partitioned the profile into five equally sized groups (folds) for cross-validation. In each of the cross-validation experiments for host h , we used as the training data for the model (cluster identification or model parameter estimation) all other microbe abundance profiles (columns of the species matrix M ) that do not correspond to host h along with four of the five-folds from host h . After model training, we used the (parameter) fitted model to make predictions on the abundances of the fold held out during training (see Supplementary Methods for further details). We performed this process for all five folds in each and every host.

Predictions in cross-validation experiments SBMterm·mth approach For each missing microbial abundance, we have three different scores associated to the probability that host h has microbe m with null, low or high abundance: \({\rm Pr}[a_{hm}=\rm I]\) , \({\rm Pr}[a_{hm}=\rm II]\) , and \({\rm Pr}[a_{hm}=\rm III]\) . Our estimation in each case is given by the maximum of these probabilities, i.e. (11) \[\begin{eqnarray*} \rm Prediction = \rm \max \limits _x \lbrace \rm Pr[a_{hm} = x] \rbrace \,. \end{eqnarray*}\] We, thus look for the membership vectors \(\boldsymbol{\theta }\) and \(\boldsymbol{\eta }\) that yield the highest predictive power of unobserved abundances.

Conventional approach The agglomerative clustering pipeline we consider rigorously selects the best out of 90 possible clusterings. For each host h and each fold, we select the best clustering to estimate unobserved abundances. We first consider all the hosts that belong to the same cluster as host h . Then, for each missing microbe m i , we estimate its relative abundance as the most likely abundance of that microbe in the other hosts within the cluster.

Similarity distance between membership vectors of microbes according to taxonomic classification Suppose we have a taxonomic level T with t 1 , …, t q taxa. For each taxa t k , we take all possible pairs of microbes’ membership vectors \((\eta ^i_k, \eta ^j_k)\) , such that i, j ∈ t k , and i ≠ j and compute the Euclidean distance d ij ( t k ) between them. Note that the smaller the distance, the more similar the membership vectors of the two microbes. We repeat this process for all n taxa and then compute the mean intrataxa distance as (12) \[\begin{eqnarray*} \langle d_{{\it intra}} \rangle =\frac{1}{N^P_{{\it intra}}}\sum ^q_{k=1} \sum _{i \ne j} d_{ij}(t_k) \ , \end{eqnarray*}\] where \(N^P_{{\it intra}}\) is the total number of pairs of microbes from the same taxa. Next, we measure the mean intertaxa distance. To that end, we take all possible pairs of membership vectors from different taxa \((\eta ^i_k, \eta ^j_{k^{\prime }})\) , such that k ≠ k ′ and measure the Euclidean distances between them. We repeat this process for all pairs of different taxa and calculate the mean: (13) \[\begin{eqnarray*} \langle d_{{\it inter}} \rangle =\frac{1}{N^P_{{\it inter}}}\sum ^q_{k \ne k^{\prime }} \sum _{i,j} d_{ij}(t_k, t_{k^{\prime }}) \ . \end{eqnarray*}\] Finally, we take the log of the ratio intra/inter distances log (〈 d intra 〉/〈 d inter 〉) for all five taxonomic levels.

Assessment of nestednessterm·key In order to measure the nestednessterm·key of a binary matrix in which elements are either 1s or 0s, we use the metric presented in ( 36 ): (14) \[\begin{eqnarray*} n=\frac{ \sum _{i\lt j}q^{(h)}_{ij} + \sum _{i\lt j}q^{(m)}_{ij}}{\sum _{i\lt j} {\rm min}( q^{(h)}_{i}, q^{(h)}_{j}) + \sum _{i\lt j} {\rm min}( q^{(m)}_{i}, q^{(m)}_{j} )} \ , \end{eqnarray*}\] where q ij represents the number of shared interactions between rows/columns i and j and q i is the number of interactions of row/colum i. h and m represent rows (hosts) and microbes(columns), respectively. The numerical value of the nestednessterm·key n ranges from 0 to 1.

Genome sizes and protein-coding gene counts for microbial species Genome sizes and protein coding genes for all strains available in NCBI were retrieved for the species found to be present in the five datasets mentioned above. For each species, median genome size and median count of protein-coding genes were calculated from all of its strains (see Tables S1-3, Supplementary Material ).

Results SBMsterm·mth are predictive of microbe abundances and can identify biologically relevant groups The most common approaches for enterotypeterm·dom identification rely on clustering heuristics to identify groups (clusters) of hosts who share commonalities in gut microbiomesterm·key ( 23 , 42 , 43 ) (the exception being approaches based on multinomial mixture models ( 44 , 45 )). In contrast, we use SBMsterm·mth as generative models to identify group patterns in microbe abundance profiles. SBMsterm·mth first assume that there are distinct groups of hosts and groups of microbes, and that the probability that a host has a significant abundance of a certain microbe depends exclusively on the groups to which the host and the microbe belong ( 24 ). These generative models have only a few underlying assumptions (i.e., that there are groups of microbes and groups of hosts, and that links connecting a microbe group to a host group are independent of other links), and can be used to find the most plausible group memberships of hosts and of microbes that best describe/compress a given dataset. The most plausible group memberships derived from observable data are, by definition, also the most predictive of unobserved data.

SBMsterm·mth are fully defined by the group memberships of hosts and microbes; and probability matrices , whose elements p k ℓ ( a ) determine the probability that microbes of group ℓ are present with an abundance of a in hosts of group k . In this model, each row (for k ) in the matrices \(\mathbf {p}(a)\) describes the most likely abundances of each group of microbes in host group k , and therefore, is akin to what has been typically referred to as enterotypesterm·dom (Fig. 1 ).

The same concept can be further generalized to allow hosts and microbes to simultaneously belong to multiple groups at varying probabilities (which increases the expressiveness of the model ( 28 , 35 )); and this is the type of SBMterm·mth we consider in our analysis (see Fig. 1 and Methods for a formal definition of the model). SBMsterm·mth have the additional advantage that they are amenable to efficient and well-grounded inference methods that can be used to optimally divide hosts and microbes into groups ( 27 , 28 ) (see Methods).

To show how well SBMsterm·mth can identify group patterns in microbe abundance matrices compared to previously demonstrated clustering approaches (e.g. for enterotypeterm·dom detection), we assess each model’s ability to predict relative abundances of unobserved microbes using five published gut microbiometerm·key datasets (Table 1 ). To this end, we first discretize the relative abundances r (which span a simplex distribution) into three discrete categorical bins: negligible ( r < 10 −4 ), low (10 −4 ≤ r < 10 −3 ), and high ( r ≥ 10 −3 ). The merit of this discretization technique is that it simplifies our prediction task into a classification problem. Then, we perform a host-specific five-fold cross-validation for all hosts within each dataset. Here, our main goal is to predict abundance categories of unobserved microbes in individual hosts as a way to test whether it would be possible to reconstruct the entirety of partially observed microbial abundance profiles. More specifically, we devised a cross-validation strategy that splits a profile of a given host into five equally sized partitions. For each split, we train a model with 80% (four of five partitions) of the microbe abundance profile of that host and with the entire abundance profiles of the remaining hosts; the held out 20% of the profile (from the same host) was used as the test set to evaluate predictive performance. The average prediction accuracy found in this manner provides us a means to objectively and quantitatively compare the robustness of group patterns identified by the SBMterm·mth and other clustering approaches ( 25 , 26 , 46 ).

To compare against the predictive performance of the SBMterm·mth approach, we use an agglomerative hierarchical clustering method described in ( 42 ), as such methods are generally regarded to be more reliable than divisive clustering (e.g., k-means) in providing accurate prediction models. Of note, the clustering method generates 90 different clusterings and then selects the best one based on multiple metrics. In our prediction experiments, we clearly see that the predictive performance of the SBMterm·mth approach is systematically higher than that of the conventional clustering approachterm·mth (Fig. 2a ). Hence, we can conclude that SBMsterm·mth are a superior predictive model.

SBMsterm·mth accurately predict and describe microbe abundances. (a) Predictive performance in host-specific cross-validation experiments. We compare the predictive power of SBMsterm·mth and conventional approaches for enterotypeterm·dom identification using a cross-validation strategy centered on individual hosts (Methods). For each of the five different datasets indicated by the last name of the first author (number of samples), we compute the average log-ratio (over multiple folds) of the performance in predicting unobserved microbial abundances by SBMsterm·mth to that by the conventional clustering approachterm·mth (X: accuracy, precision, and recall—see Supplementary Material for definitions). Note that, in our prediction experiments of unobserved microbe abundances, the performance across different hosts is highly variable (see Figure S3, Supplementary Material ): accuracy values are in the range of [0.5, 1], while recall and precision values have a wider range of [0.25, 1]. We compute average log-ratios, which allows us to measure performance differences between methods for the same prediction experiment. Note that a log-ratio equal to zero means that both methods have the same performance; a log-ratio above zero means that the SBMterm·mth performs better than conventional heuristic approaches, and vice versa . Bars represent the standard error of the mean. (b) Similarity in microbe abundances across taxonomic ranks. For each of the five datasets, we compute the average distance between group membership vectors of microbes in the same taxa ('intra' taxonomic ranks), and compare it to the average distance between membership vectors of microbes in different taxa ('inter' taxonomic ranks; Methods and Table S1, Supplementary Material ). Log-ratios smaller than zero indicate that the average intrataxa distance is smaller than the average intertaxa distance, that is, that membership vectors of microbes within a taxon are more similar than membership vectors of microbes across different taxa. Bars represent the standard error of the mean.

Our analysis also shows that the majority of hosts and microbes are predictable; that is, that we are able to correctly predict unobserved abundances of microbes across many hosts, and also unobserved abundances in specific hosts for many microbes. Nonetheless, we also observe some variability in host and microbe predictability—some hosts and microbes have very strong patterns and are, therefore, more predictable than others ( Figures S4–S8, Supplementary Material ). Further analysis of the SBMterm·mth groups reveals that the similarity of group memberships increases as we go down the taxonomic rank hierarchy (Fig. 2b ). In particular, we find that the genus taxonomy has the largest intrataxon to intertaxon similarity ratio (Fig. 2b ), suggesting that SBMsterm·mth can infer biologically relevant groups in gut microbiomesterm·key.

SBMsterm·mth reveal the nested organizationterm·dom of microbe abundances in the human gut microbiometerm·key Our results show that SBMsterm·mth enable us to reasonably predict unobserved microbe abundances, indicating that our approach correctly models relationships between hosts and microbes in the five gut microbiometerm·key datasets. The questions we explore next are: What are these patterns? Are they consistent across datasets? Do these patterns reveal universal traits in the organization of gut microbiometerm·key ecologies? To these ends, we analyze the \(\mathbf {p}\) matrices, which describe the abundances of groups of microbes in groups of hosts, for each of the five datasets.

If we rank-order host groups according to their number of microbe groups of non-negligible relative abundances r ( r ≥ 10 −4 ), and rank-order microbe groups in an analogous manner, we find that biclusteringterm·key \(\mathbf {p}\) matrices display, what is called in ecology, a nested organizationterm·dom ( 29 , 30 , 32 , 48 ) (Fig. 3 ; see Methods and Figure S9 (Supplementary Material) for an illustration). Such a nested organizationterm·dom entails that we can observe generalistterm·key and specialistterm·key groups of hosts and microbes. Generalist hostterm·dom groups have non-negligible abundances of the majority of microbes groups, while specialist hostterm·dom groups have non-negligible abundances of very few microbe groups. Importantly, microbe groups present in specialist hoststerm·dom are subsets of those observed in generalist hoststerm·dom; and thereby the nested structure. Moreover, microbe groups display the same behavior: generalistterm·key microbe groups have non-negligible abundances in the majority of host groups, while specialistterm·key microbe groups have non-negligible abundances in a few host groups. The nested structure entails that specialist hoststerm·dom have non-negligible abundances of mainly generalistterm·key microbe groups, while specialistterm·key microbe groups are only present in generalist hoststerm·dom (see Figure S9, Supplementary Material ). Note that in our analysis, generaliststerm·key and specialiststerm·key are defined purely based on what we observe in the data, i.e. their positions within the nested hierarchy of the SBMterm·mth \(\mathbf {p}\) matrices. Therefore, these terms are not defined using the definition of generalistterm·key in other contexts.

SBMsterm·mth reveal underlying nested organizationterm·dom of microbe abundance matrices. ((a)–(e)) Nestednessterm·key of SBMterm·mth \(\mathbf {p}\) matrices for all datasets: (a) W. Liu et al ., (b) N. Qin et al ., (c) M. Schirmer et al ., (d) C. Huttenhower et al . and J. Lloyd-Price et al ., and (e) D. Zeevi et al . Top row: heatmaps represent SBMterm·mth \(\mathbf {p}\) matrices. Rows correspond to host groups, and columns correspond to microbe groups. As in ecological studies (see for instance ( 47 )), we consider two categories of microbe relative abundance r : negligible ( \(\rm 0\) , for r < 10 −4 ) and non-negligible ( \(\rm 1\) , for r ≥ 10 −4 ), and analyze the corresponding SBMterm·mth \(\mathbf {p}\) matrices. Furthermore, because in this case the values \(p_{k\ell }(a=\rm 1)\) matrices are very close to either 0 or 1, we round their values. Note that because our optimal choice of parameters is K = 10 groups of hosts and L = 20 groups of microbes, \(\mathbf {p}\) matrices are 10 × 20 (see Methods). Dark blue matrix elements represent the presence of a microbe group in a host group with non-negligible abundance. Rows and columns are ordered according to their count of non-negligible abundances to show the nested organizationterm·dom (see Figure S9, Supplementary Material ). The inset in each plot shows the numerical value of the nestednessterm·key (see Methods). Bottom row: significance of nestednessterm·key for the corresponding heatmap shown in the top row. For each dataset, we show the distribution of nestednessterm·key values obtained for 1,000 randomizations of the \(\mathbf {p}\) matrix. For this, we randomize each \(\mathbf {p}\) matrix while preserving the average number of connections between groups of hosts and microbes. The black solid line represents the nestednessterm·key value of the observed (actual) \(\mathbf {p}\) matrix.

Strikingly, we find that the nested organizationterm·dom in the SBMterm·mth \(\mathbf {p}\) matrices is also present in the original host–microbe abundance matrices (Fig. 4 ). That is, if we order hosts and microbial species by their degree (the number of non-negligible abundances in their corresponding rows and columns), the resulting ordered relative abundance matrices also display a nested structure (Fig. 4 ), allowing us to again identify generaliststerm·key and specialiststerm·key. Importantly, our results do not depend on the threshold we use to define non-negligible and negligible relative abundances (see Figure S10, Supplementary Material ); even after changing the relative abundance threshold to below what was originally considered as negligible/non-negligible, nested hierarchies in the microbe abundance matrices are still observed.

Microbe abundance matrices display a nested organizationterm·dom of hosts and microbes ((a)–(e)) Nestednessterm·key of microbe abundance matrices from all of the datasets considered: (a) W. Liu et al ., (b) N. Qin et al ., (c) M. Schirmer et al ., (d) C. Huttenhower et al . and J. Lloyd-Price et al ., and (e) D. Zeevi et al . Top row: categorical matrices with two categories of microbial species relative abundance r : negligible ( \(\rm 0\) , for r < 10 −4 , white) and non-negligible ( \(\rm 1\) , for r ≥ 10 −4 , blue). Rows represent individual hosts, and columns represent individual microbial species. Microbial species and hosts have been ordered according to the number of non-negligible abundances in their corresponding columns and rows, revealing a nested pattern (see Figure S9, Supplementary Material ) as previously seen for the host and microbe groups. The inset in each plot shows the numerical value of the nestednessterm·key (see Methods). The bottom plots show the nestednessterm·key of 1,000 randomizations of the datasets. Black solid lines represent the observed (actual) nestednessterm·key value in the host-microbe abundance matrices above. Randomizations preserve the average number of connections between hosts and microbes.

Ecological roles of microbial species are robust In ecology, generalistterm·key species are considered to have a wide dietary breadth (i.e. species that can survive from a variety of food resources) and are less affected by environmental perturbations ( 49 ), whereas specialiststerm·key have a narrower dietary breadth and can exhibit a competitive fitness advantage over generalistterm·key species within the confines of their habitat ( 50 , 51 ). We now investigate whether microbial species in the gut can play the same roles. Indeed, we find that extreme roles (i.e. extreme generaliststerm·key and specialiststerm·key) are well conserved—generalistterm·key and specialistterm·key species play similar ecological roles in all datasets (Fig. 5a ). Specifically, we find that generalistterm·key species concentrate in the Clostridiales and Bacteroidales orders ( Firmicutes and Bacteroidetes phylum, respectively; see Tables S1 and S2, Supplementary Material ), which is coherent with what enterotypeterm·dom identification studies have found to drive enterotypeterm·dom formation ( 23 ). Many specialistterm·key species, such as Anaerotruncus colihominis or Cenarchaeum symbiosum , are also well-conserved across datasets. A large part of these conserved specialistterm·key species belong to the Firmicutes phylum, including many in the Clostridiales order ( Tables S1 and S2, Supplementary Material ).

Ecological roles of microbial species are robust but not related to genome size or number of protein-coding genes. (a) Generalistterm·key and specialistterm·key roles are conserved across datasets. To every microbial species in a nested relative abundance matrix (Fig. 4 ), we assign a rank equal to the position of the column occupied by that microbial species. Low rankings correspond to generalistterm·key species, while high rankings correspond to specialistterm·key species (see Table S2 (Supplementary Material) for normalized ranks of each species in each dataset). For each microbe present in at least four of our datasets, we compute the mean normalized rank in the ordered abundance matrices in Fig. 4 and the standard error of the mean normalized rank. (We normalize the rankings to control for the different number of microbial species in each dataset.) The standard error tends to be smaller for species with low (generaliststerm·key) and high (specialiststerm·key) normalized rank, thus showing that generalistterm·key and specialized microbial species are conserved across datasets. In fact, the variance in normalized ranks across datasets is small for the vast majority of species, which implies that microbes tend to occupy the same position within the nested hierarchy. The line is a quadratic fit to the data and serves to guide the eye. (b)–(f) Relationship between genome size (top) and protein-coding gene count (bottom) with respect to the rank of a microbial species (see Methods and Tables S2 and S3, Supplementary Material ). At the top of each panel, we show Spearman’s ρ and its corresponding P -value. There is no general trend in the datasets (except for the smallest dataset, panel (c)). This indicates that generalistterm·key and specialistterm·key species in nested matrices do not have obvious differences in genome size or the number of protein-coding genes. Study datasets: (b) Liu et al .; (c) Qin et al .; (d) Schirmer et al .; (e) Huttenhower et al . and J. Lloyd-Price et al .; and (f) Zeevi et al .. Points in the scatter-plots indicate individual microbiome samples.

To investigate whether the ecological role of microbial species in the nested hierarchy (i.e. generaliststerm·key and specialiststerm·key) could be connected to their potential functional capacity, we looked at the relationship between genome sizes (and also number of protein-coding genes) of microbes and the ranks of microbes in the rank-ordered abundance matrices in Fig. 4 . Here, low and high rank correspond to generalistterm·key and specialistterm·key microbes, respectively. We find no robust relationship between genome size (Fig. 5 b–f top row) and ecological role, or between the number of protein-coding genes and ecological role (Fig. 5 b–f bottom row; see Figure S11, Supplementary Material ). This suggests that the ecological role microbes play in a nested organizationterm·dom may be independent of the range of functions a given species could potentially perform, but rather related to the occupation of the specific habitat.

Generalistterm·key and specialist hoststerm·dom show varying gut microbiometerm·key diversity Island biogeography studies show that a consequence of the nested structure in fragmented habitats and archipelagos is wide disparities in species richness across habitats ( 29 , 48 ). For individual hosts, microbiome composition is often quantified in terms of a single alpha-diversity index ( 52 ). These diversity indicators have often been regarded as a possible indicator of gut health—the higher the diversity index, the healthier the host ( 3 ). Our results show that, within the nested organizationterm·dom of gut microbiometerm·key datasets, there is a significant correlation between the host’s rank (occupied row number) in the rank-ordered nested abundance matrix in Fig. 4 and that host’s Shannon diversityterm·key index. More specifically, generalist hoststerm·dom have a higher diversity than specialist hoststerm·dom (Fig. 6 ), confirming our earlier finding that generaliststerm·key hosts provide a habitat for a wider range of gut microbes. In all, we conclude that ecological nestednessterm·key, which has been reported in past works in macroecology, extends to the human gut microbiota. Further examinations into this property could lead to new discoveries regarding the ecological assembly of the human gut microbiometerm·key.

Relationship between Shannon diversityterm·key index and host rank. For each dataset, we show the Shannon diversityterm·key index vs. the rank of each host in the ordered relative abundance matrix (i.e. the position of the row which corresponds to a host in the ordered abundance matrix). Low ranks correspond to generalist hoststerm·dom and high ranks correspond to specialist hoststerm·dom. We also show the Spearman’s ρ and its corresponding P -value. In all datasets, we observe a significant inverse correlation between the Shannon diversityterm·key index and host rank. Study datasets: (a) Liu et. al , (b) Qin et al .l, (c) Schirmer et al ., (d) Huttenhower et al . and Lloyd-Price et al ., and (e) Zeevi et al . Points in the scatter-plots indicate individual gut microbiometerm·key samples.

Discussion Ever since links between gut microbiometerm·key dysbiosis and complex, chronic diseases were established, there has been a keen interest in understanding what a healthy gut microbiometerm·key would look like, and how therapies can be designed to restore gut dysbiosis to a healthy state. Fuelled by these interests, scientists have aimed to uncover recurrent patterns in gut microbiometerm·key composition by categorizing gut microbiomesterm·key into enterotypesterm·dom. However, despite some broadly shared properties ( 18 , 23 ), no universally recognized taxonomic property in the gut microbiomesterm·key of healthy subjects has been found. In this study, we demonstrate that gut microbiomesterm·key can indeed show patterns of organization that are common across different studies. Notably, these patterns are defined in terms of the ecological organization when considering both hosts and microbes simultaneously rather than hosts or microbes separately.

We used SBMsterm·mth to formally identify robust organization principles of the healthy human gut microbiometerm·key. By focusing on both groups of microbes and hosts simultaneously, we found that microbe abundance matrices display a nested structure that elucidates well-defined ecological roles (i.e. generaliststerm·key vs. specialiststerm·key). Some hosts and microbial species are generaliststerm·key: generalistterm·key microbes are present in non-negligible abundances in most hosts, while generalist hoststerm·dom contain the majority of microbial species in non-negligible abundances. In contrast, some other hosts and microbes are specialiststerm·key: specialistterm·key microbes are only present in generalist hoststerm·dom, and specialist hoststerm·dom only have non-negligible abundances of generalistterm·key microbes.

There has been much interest in the origin of nested ecologies, as such nested organizationterm·dom has been observed in a number of ecological systems ( 30 ). It was first reported in faunas of insular fragmented habitats and in archipelagos ( 29 , 48 ). Other cases of nested hierarchies were found in endoparasite communities in fish ( 31 ) and parasite communities in geographically dispersed species of bats ( 33 ). Furthermore, nestednessterm·key has been found in ecosystems with mutualistic interactionsterm·dom ( 30 , 32 , 53 , 54 ), i.e. ecological relationships benefiting both parties. Pollinator–plant networks are a common example of mutualistic networks; pollinators obtain their food while pollinizing the plants. Such nested ecological structures have even been suggested to play a role in optimizing and balancing species abundances in the presence of mutualistic interactionsterm·dom ( 55 ). Despite much progress in its study, however, the relationship between mutualism, nestednessterm·key, stability, and species diversity is still relatively poorly understood and remains a subject of close and careful investigation ( 36 , 56–58 ).

We can establish several connections between human gut microbiomesterm·key and other nested ecosystems. First, the relationship between hosts and their microbiomes is well-known to be mutually beneficial. A prominent example is in food digestion, in which microbes in the gut, such as those in the Bacteroidetes phylum, aid in the degradation of dietary complex polysaccharides ( 15 , 59–61 ). Second, we can draw straightforward analogies with findings in the macroecology literature. For the case of island biogeography, hosts can be viewed as islands and microbiota as the island fauna. Several mechanisms have been put forward to explain the nestednessterm·key observed in island faunas, such as selective extinction, colonization ability, habitat nestednessterm·key, and human intervention ( 30 , 47 , 62 , 63 ). In regard to the gut microbiometerm·key, all of these explanations can offer promising clues and insights into how the microbes assemble and respond to perturbations. For instance, Cook and Quinn ( 47 ) found that, in collections of islands, fauna with large dispersal ability (such as birds) show more nested abundance patterns across islands than reptiles and mammals. Their finding could possibly suggest that, in a gut microbiometerm·key context, a gradient in colonization ability (in which dispersal ability is a strong factor) may contribute to the nested organizationterm·dom in microbial species abundances observed in our study.

Other studies have shown that when colonizing new patches, animals with larger trophic breadths (food generaliststerm·key) build the first successful colonies; and thereafter, food specialiststerm·key (typically predators that feed on a narrow range of prey) arrive at a later stage ( 64 , 65 ). An explanation of these observations is that species with large trophic breadths are able to survive comparatively well in ecosystems subject to fluctuations, whereas food specialiststerm·key need more mature, developed, and diverse ecosystems to survive. By drawing a direct analogy with these observations, our results suggest that specialist hoststerm·dom have gut environments that offer relatively restricted living conditions for microbes. In summary, our study offers the opportunity to bridge the gap between macroecology and gut microbiometerm·key studies, and thereby shines a new perspective on human-associated microbial ecology research.

As a methodological point, we note that inferring clusters (or microbial communities) based on Dirichlet multinomial processes or other probabilistic mixture models is indeed a common approach for microbial community detection from metagenomic data (whether the inferred distributions are biologically meaningful or accurate is a different issue) ( 44 , 45 ). To clarify, the Dirichlet multinomial mixture modelterm·mth pipeline developed by Holmes et al . ( 45 ) was designed to find clusters of microbial communities, whereas our SBMterm·mth framework is a biclustering techniqueterm·mth that allows simultaneous clustering of both microbes and hosts; in addition, the Dirichlet model was not designed for the predictive tasks we used to validate our approach, namely, prediction of held-out microbial relative abundances. Specifically, as part of its input, the Dirichlet model needs all the metagenomic reads in each host’s microbiome sample to give the inferred multinomial model parameters. However, its current implementation does not work when a subset of the microbial species’ reads are held-out in cross-validation, as we had done to compare our SBMterm·mth approach against the agglomerative clustering pipeline. While it would be possible to program a different implementation of the Dirichlet model specifically intended for our prediction experiments, this task falls outside of the scope of our current study.

To conclude, we demonstrate that the human gut microbiometerm·key does indeed have robust structural patterns underlying its taxonomic composition. To this point, our study may contribute toward defining the design principles of a healthy gut microbiometerm·key, developing better models to understand gut microbiometerm·key evolution, and establishing strategies to restore or maintain wellness in the form of probiotic communities. An open question is whether diseases can disrupt the nested structure of healthy microbiomes. If that were the case, nestednessterm·key would be a straightforward metric of the health state of gut microbiomesterm·key, and interventions could be designed to restore the original ecological pattern. Explorations into these topics will be the subject of our future studies.

Supplementary Material Click here for additional data file.