Microbial Co-occurrence Relationships in the Human Microbiome

1In nature, organisms rarely live in isolation, but instead coexist in complex ecologies with various symbiotic relationships. As defined in macroecology, observed relationships between organisms span a wide range including win-win (mutualism), win-zero (commensalism), win-lose (parasitism, predation), zero-lose (amensalism), and lose-lose (competition) situations,,. These interactions are also widespread in microbial communities, where microbes can exchange or compete for nutrients, signaling molecules, or immune evasion mechanisms,,. While such ecological interactions have been recently studied in environmental microbial communities,,,, it is not yet clear what the range of normal interactions among human-associated microbes might be, nor how their occurrence throughout a microbial population may influence host health or disease.

2Several previous studies have identified individual microbial interactions that are essential for community stability in the healthy commensal microbiota,,,, and many are further implicated in dysbioses and overgrowth of pathogens linked to disease. Each human body site represents a unique microbial landscape or niche,, and relationships analogous to macroecological “checkerboard patterns” of organismal co-occurrence have been observed due to competition and cooperation,,,. For example, dental biofilm development is known to involve complex bacterial interactions with specific colonization patterns,,. Likewise, disruption of relationships among the normal intestinal microbiota by overgrowth of competitive pathogenic species can lead to diseases, e.g. colonization of Clostridium difficile in the gut. However, no complete catalog of normally occurring interactions in the human microbiome exists, and characterizing these co-occurrence and co-exclusion patterns across body sites would elucidate both their contributions to health and the basic biology of their ecological relationships. Thus, characterizing key microbial interactions of any ecological type within the human body would serve as an important first step for studying and understanding transitions among various healthy microbial states or into disease-linked imbalances.

3As has been also been pointed out in macroecology, however, the analytical methodology needed to comprehensively detect such co-occurrence relationships is surprisingly complex. Most existing studies employ simple measures such as Pearson's or Spearman's correlation to identify significant abundance relationships,,. These methods are suboptimal when applied without modification to organismal relative abundances. Since absolute microbial counts are not known and measurements depend on sampling and sequencing depth, an increase in one relative abundance must be accompanied by a compositional decrease in another, leading to spurious correlations among non-independent measurements. In addition, sparse sequence counts can cause artefactual associations for low-abundance organisms with very few non-zero observations. Conversely, association methods such as log-ratio based distances that have been developed specifically for such compositional data are difficult to assign statistical significance, a vital consideration in high-dimensional microbial communities containing hundreds or thousands of taxa.

4Here, we have addressed these issues to catalog a baseline of normal microbial interactions in the healthy human microbiome. The Human Microbiome Project (HMP) sampled a disease-free adult population of 239 individuals, including 18 body habitats in five areas (oral, nasal, skin, gut, and urogenital), providing 5,026 microbial community compositions assessed using 16S rRNA gene taxonomic marker sequencing. We have developed a suite of methods to characterize microbial co-occurrence and co-exclusion patterns throughout the healthy human microbiome while suppressing spurious correlations. Specifically, these were 1) an ensemble approach including multiple similarity and dissimilarity measures, and 2) a compendium of generalized boosted linear models (GBLMs) describing predictive relationships, both assessed nonparametrically for statistical significance while mitigating the effects of compositionality. Together, these methods provide a microbiome-wide network of associations both among individual microbes and between entire microbial clades.

5Among the 726 taxa and 884 clades in the HMP data, we examined both intra-body site and inter-body site relationships as a single integrated microbial co-occurrence network. Each relationship represents co-occurrence/co-exclusion pattern between a pair of microbes within or between body sites among all subjects in the HMP (in contrast to studies within single subjects of microbial co-occurrences across biogeography, e.g., ). This ecological network proved to contain few highly connected (hub) organisms and was, like most biological networks, scale-free. Co-occurrence patterns of the human microbiome were for the most part highly localized, with most relationships occurring within a body site or area, and there were proportionally few strong correspondences spanning even closely related body sites. Each pair of organisms was assessed for positive (e.g. cooperative) or negative (e.g. competitive) associations, and in many cases these patterns could be explained by comparing the organisms' phylogenetic versus functional similarities. In particular, taxa with close evolutionary relationships tended to positively associate at a few proximal body sites, while distantly related taxa with functional similarities tended to compete. The resulting network of microbial associations thus provides a starting point for further investigations of the ecological mechanisms underlying the establishment and maintenance of human microbiome structure.

6We inferred a microbiome-wide microbial interaction network by analyzing 5,026 samples from the Human Microbiome Project (HMP) comprising 18 body sites, 239 individuals recruited at two clinical centers, and 726 bacterial phylotypes detected by 16S rRNA gene sequencing (Table 1). Our study aimed to determine co-occurrence and co-exclusion relationships among the relative abundances of microbial taxa across all individuals, potentially indicative of their ecological relationships. We thus combined two complementary approaches, namely an ensemble of multiple similarity and dissimilarity measures (henceforth “ensemble approach”) and a compendium of generalized boosted linear models (GBLMs, henceforth “GBLM approach”). Both methods were applied to the HMP data to produce microbial interaction networks in which each node represented a microbial clade (taxon or group of taxa) connected by edges that were weighted by the significance of their association (positive or negative). Spurious correlations due to compositional structure of relative abundance data were prevented by a novel bootstrap and re-normalization approach assessing the degree of association present beyond that expected by compositionality alone. We used Simes method followed by Benjamini-Hochberg-Yekutieli false discovery rate (FDR) correction to combine the resulting networks (Figure 1). A detailed final network is provided in Figure S1, with a comparison of all networks in Figure S7 and additional information in Methods. This provided a single global microbial interaction network capturing 3,005 associations among 197 phylotypes, spanning all available body sites from the human microbiome (Figure 2; Table S1).

7Global properties of the microbiome-wide network of microbial associations are summarized in Figures 2 and 3. A dominant characteristic of the network was its habitat-specific modularity. After grouping the 18 body sites into five broad areas (oral, skin, nasal, urogenital, and gut), the large majority of edges were found clustered within body areas (98.54%), and these clusters were sparsely connected through a minority of edges (1.46%). This is confirmed by the network's high modularity coefficient of 0.28 (as defined by ) and Markov clustering of the network (see Methods and Figure S2). It has long been observed that sites within the human microbiome are distinct in terms of microbial composition, and this proved to be true of microbial interactions as well: microbial relationships within each body area's community were largely unique (Table 2). The microstructure of interaction patterns - and thus in the underlying ecology - was different for different areas, however. For example, all vaginal sites within the urogenital area were interrelated in a single homogeneous community, whereas interactions within the oral cavity suggested microbial cross-talk among three distinct habitats. This can be observed quantitatively based on the proportions of microbial interactions spanning body sites within each area, e.g. 69.57% among the vaginal sites and 53.19% among the oral sites, both exceeding the microbiome-wide baseline. The skin was further unique in that the large amount (57.65%) of its associations related microbes in corresponding left and right body sites (left and right antecubital fossae and retroauricular creases), reflecting consistent maintenance of bilateral symmetry in the skin microbiome.

8We began decomposing the network by categorizing microbial associations within each body area into body-site-specific relationships of two types: cross-site and within-site interactions. On average, these two classes make up 53.11 and 46.89 percent of the total edges, respectively (Table 2). First focusing on cross-site associations, a majority (66.10%) of such relationships were co-occurrences between the same or taxonomically related clades in proximal or bilateral body sites. This reflects coordinated community structure among ecologically related niches, such as similar dental plaques, vaginal sites, and bilateral skin sites. Body sites specifically connected by many positive associations were either in direct contact (e.g. tongue and saliva), proximal (e.g. sub- and supragingival plaques), or similar in terms of environmental exposure (e.g. bilateral skin sites), thus providing mechanisms to support comparable microbiota and exhibiting high levels of microbial co-occurrence. This pattern held true for the minority (33.90%) co-exclusions as well, with many occurring between bilateral skin sites or within subgroups of the oral cavity. This suggested that the first level of hierarchical co-occurrence structure in this network corresponded with groups of body sites representing distinct microbial habitats.

9Conversely, within-site relationships showed a much more balanced ratio of microbial co-occurrence (48.26%) vs co-exclusion (51.74%) interactions. Many of the negative within-site relationships were associated with the abundant signature organisms characteristic of each body site, for example Streptococcus in the oral cavity and Bacteroides in the gut. The relative abundances of these signature taxa varied greatly among individuals, in some cases (e.g. Bacteroides) spanning from 1% to 97% within a body site across the HMP population. It is generally very difficult to determine from relative abundance measurements alone whether these negative associations represent true anti-correlation (e.g. one organism out-competing another) or overgrowth of one organism while the rest of the population remains unchanged (resulting in a negative correlation due to compositionality of these data). This problem has a long history in quantitative ecology,. Our methods generally determine these relationships in the human microbiome to be stronger than what would be expected from compositionality alone (see Methods and Text S1), and the negative interactions detected here are thus likely biologically informative. This is supported by the fact that they are strongest in cases where distinct alternative dominant community members occurred among different individuals (e.g. Prevotellaceae vs. Lactobacillaceae in the vaginal area or Propionibacterium vs. Staphylococcus on the skin, ). The increase in negative interactions within habitats is also in line with the fact that most competitive mechanisms require proximity or physical contact, whereas positive interactions are likely to also occur from microbiome-wide shared environmental exposures.

10We further assessed several other measures of network community structure. Globally speaking, the network followed a scale-free degree distribution typical of biological systems, meaning that most clades possessed few interactions but a few clades possessed many (Figure 3A ), The network had a low average path length of three (contrasted with six in randomized networks), meaning that short paths existed between most clades, and it possessed a low average per-node cluster coefficient (0.1) measuring the local density of connections. Together, these values indicate that the microbial association network is structured to be scale-free and thus robust to random disruption, with only sparse local multi-organism clusters. Since these data only describe phylotypes at approximately the genus level, it remains to be seen whether a greater degree of locally clustered functional associations emerges among Operational Taxonomic Units (OTUs), species, or strains within these phylotypes. As the cluster coefficient distribution was not well described by the inverse node degree distribution, the network possesses no strong hierarchical modularity despite its scale-freeness, in contrast to the strong habitat-centric modularity.

11The diversity of microbial interactors (i.e. number of unique phylotypes) within each body site also proved to directly dictate its interaction density (Figure 3B). That is, communities with a greater number of different organisms had a proportionally greater number of positive and negative associations. Within these sites, the number of relationships scaled directly with the number of unique phylotypes (adjusted R2 of 0.75), the only body site with more interactions than expected for its diversity being the tongue dorsum (see also Table S2). This site also harbored the top-ranking hub phylotype (Firmicutes, see Figure 3A). In combination with the behavior of specific microbial hubs as discussed below, this might argue that most microbial taxa form strong metabolic or functional associations with adjacent taxa inhabiting the same body site habitat, allowing consortia to specialize within highly localized microbial niches.

12When randomizing between rather than within body sites, no body site pairs possessed more cross-site associations than expected (with the slight exception of tongue dorsum), whereas most body sites were significantly enriched for within-site relationships (the only exceptions being posterior fornix, mid-vagina, and antecubital fossae, which tended toward too few phylotypes to reach significance; see Figure 3D and Table S2), again confirming the microbiome's habitat-driven modularity. When calculating network properties in a body-area-specific manner, we found that the overall average path length between nodes in the oral cavity, which contributes most of the samples, was much larger (∼3.4) than those of the other body areas (ranging from ∼1.1 to ∼2.0). In addition to supporting the aforementioned degree of inter-site habitat formation in the oral cavity, this intriguingly suggests that other body sites in which fewer samples are currently available (see Table 1) have not yet exhausted the detection of microbial relationships in the human microbiome. More samples and greater sequencing depth may further improve detection power.

13We next examined the associations of individual clades with respect to interaction degree, observing highly connected “hub” clades to be found within each body area. Two classes of hubs appeared in the association network: clades highly connected within one body site, and clades acting as “connectors” between multiple body sites. Hubs included both specific taxa (e.g. Porphyromonas, see Figure 3A, Table S3) and larger taxonomic groupings (e.g. the phylum Firmicutes). Within-site hubs were often, although not always, abundant signature taxa (detailed below), high-degree exceptions including Atopobium on the tongue (28 total associations, 16 within-site) and Selenomonas on both tooth plaques (20 total/19 within and 7 total/3 within for supra- and subgingival, respectively). The latter provides a striking example of the niche-specificity of these low-abundance within-site interactors, as Selenomonas averages only 1.1% and 1.2% of the sub- and supragingival plaque communities, respectively, but associates preferentially (20 of 27, 74%) with members of the greater oxygen availability supragingival community. The clade's detection as a within-site hub thus corresponds with the ecology that might be expected of an organism known to be oxygen-sensitive, fastidious, and grown best in co-culture.

14Between-site hubs typically operated among body sites within the same area as described above, with two of the five most connected hub clades in the network falling into this connector category linking multiple body sites, Firmicutes and Proteobacteria on the tongue (see Figure 3A). The Firmicutes and Porphyromonas (phylum Bacteroidetes) hubs in the tongue also had the largest numbers of negative connections among all phylotypes, and all of these highly interactive clades centered on the tongue and spanned multiple related oral habitats. Signature clades such as the Firmicutes are of course highly functionally diverse, and this network suggests that the few abundant members in any one habitat might instead serve as “information processors” throughout a body area. In contrast to the low-abundance within-site hubs, this would allow them to provide baseline functionality complemented by distinct, less abundant clades with which they co-occur within differing body site habitats.

15Correspondingly, Firmicutes and other inter-site hub nodes showed a higher connectivity than the clades with highest intra-site degree (e.g. Bacteroidales in the subgingival plaque). Such clades with unusually frequent inter-site associations are thus outliers relative to the network's overall habitat-specific trend and suggest that inter-site hubs are particularly critical for associating similar sites within the same body area. In the oropharynx, for example, Streptococcus spp. with a modest degree of functional variation might be present throughout the habitat, interacting with distinct, more specialized clades within each body site. Almost all such high-connectivity hubs occurred among oral sites (e.g. Porphyromonas, Streptococcus, Veillonella, and others), the first notable exception being the Propionibacterium hub on skin sites (left and right retroauricular crease). All of these follow the same pattern, however, in which abundant phylotypes likely possessing within-clade functional diversity are distributed among related habitats within each individual.

16We additionally examined the phylogenetic rather than biogeographical distribution of these associations, testing whether clades tended to support more phylogenetically related (within-clade) or diverse (between-clade) interactions. We first investigated purely quantitative degree distributions by summarizing clades at the class level. Associations were summarized as the fraction of all possible interactions that were observed to occur, separated into positive and negative bins (Figure 3E). In addition, clade-specific over-representation of these bins was tested for significance by randomization (see Methods and Table S2). The only classes that showed significantly more negative (and, simultaneously, cross-clade) associations than expected were the Bacteroidia, Bacilli, and Fusobacteria. Most of the common classes in the human microbiome had more intra-clade edges than expected by chance (Actinobacteria, Bacilli, Bacteroidia, Betaproteobacteria, Clostridia, Epsilonproteobacteria, Fusobacteria, Gammaproteobacteria, Mollicutes, and Spirochaetes), most of which also have high cluster coefficients (Figure S3). Taken together with the biogeographical interactions assessed above, the enrichment for within-class associations likely indicates a phylogenetic aspect of the same behavior. Specifically, if one member of such a class is abundant in one body site within an individual, it (or closely related class members) also tends to be enriched in related body sites.

17We next considered relationships between class-level clades throughout the microbiome, summarized in Figure 4. Surprisingly, the Actinobacteria and Bacilli form only co-exclusion relationships with other classes, most strongly with Bacteroidia and Fusobacteria, and primarily within the oral cavity. These clades (which include the extremely abundant streptococci) might thus be largely self-sufficient in the functional diversity needed to maintain an oral community, excluding other clades when appropriately supported by e.g. environmental factors. Although a few classes were linked by positive as well as negative interactions (e.g. Clostridia and Bacteroidia), none of these reached significance on randomization. Classes connected by both positive and negative links might suggest either that the clades exhibit co-occurrence only in some environments or that some members of the two classes co-occur while others co-exclude. As the oral communities are both the most data-rich and the most alpha-diverse in the human microbiome, it is not surprising that most relationships are observed within and among them. For instance, 97% of the specific mutual exclusions between Bacilli and Bacteroidia members occur in oral sites, as do 81% of the members of the Clostridia and Bacteroidia. The second largest contribution to the latter exclusion (∼18%) comes from the gut, reflecting the frequently discussed Bacteroides/Firmicutes ratio observed in Western populations,, and similar tradeoffs (with few positive associations) were observed in other habitats such as the skin (e.g. Staphylococcus in the Bacilli and Propionibacterium in the Actinobacteria ).

18Co-exclusions such as these have previously been observed in the human microbiota to induce distinct alternative community configurations, which may differ across persons, as well as time points (e.g. early and late colonizers in community establishment or repopulation after disturbance). Although our methodology does not explicitly describe alternative community configurations, co-occurrence networks can in some cases capture them as extreme exclusion relationships between key microbial taxa. For instance, Ravel et. al reported five different vaginal communities in an independent cohort of healthy women, four dominated by Lactobacilli and the fifth diverse and featuring members of the Actinobacteria, Clostridia, Bacteroidia, and other classes. These alternative configurations occur as mutual exclusions in our genus-level phylotypes between Lactobacillus and members of this fifth diverse community (particularly anaerobes such as Anaerococcus and the Prevotellaceae). Furthermore, we see a strong negative correlation in stool samples between Bacteroides and members of the gut community, including the Ruminococcaceae and other Firmicutes. In other body sites, the clade relationship network (Figure 4) features a negative interaction between Bacilli and Bacteroidia classes that mostly occurs in the oral cavity, and oral Porphyromonas (a member of the Bacteroidia) is among the most highly connected negative hubs. Porphyromonas is abundant (avg. 3.3% s.d. 3.9%) in oral habitats but not in most cases the dominant clade; the clade also includes potential oral pathogens, and this may be one of the more striking examples of functional competition and co-exclusion occurring with a specific clade among several oral communities.

19The digestive tract is home to one of the most diverse and densely populated microbial communities in the human body. Oral sites made up half of the body sites surveyed here, as well as exhibiting the greatest within-subject microbial diversity. Correspondingly, associations between microbes within and among oral sites likewise comprised the majority (86.46%) of all edges in our co-occurrence network, also forming its largest connected component. This consisted of two clusters of organisms from the mouth soft tissues (gingiva, mucosa, and palate) and distal areas (tongue, throat, tonsils, and saliva); the oral hard surfaces (sub- and supra-gingival plaques) formed an additional isolated habitat that showed significantly fewer microbial associations with the remainder of the oral cavity (Figure 5). A complementary analysis of the HMP microbiomes has revealed evidence of three sub-habitats within the oral cavity based on overall similarity of their microbial communities, and these results demonstrate that the shared community structures of these habitats were to a lesser degree recapitulated in terms of specific microbial associations (see Figure 5 below).

20Although the current study is associative and does not by itself establish causative mechanisms of interaction for these microbial associations, many that we detect in the oral cavity in particular are supportd by known metabolic or biochemical interactions. For instance, in the context of cell to cell interaction, Fusobacterium species are known to be bridging organisms in the development of oral biofilms by co-aggregation through physical contact. This bridging occurs during biofilm maturation, allowing a more complex use of resources including sugars (the predominant carbon source for early colonizers) and proteins (used by late colonizers). In the hard palate, for example, positive associations were found between Fusobacterium and Capnocytophaga, Peptosptreptococcus, and Porphyromonas, which are in agreement with previously published cell-to-cell interactions,, and these predictions additionally implicate Leptotrichia and Parvimonas. Dental plaque associations included Parvimonas, Prevotella, and Treponema, also in agreement with existing evidence. However, those previously published aggregations are strain specific and, this study may be observing broader effects than the direct cell-cell contact preferences in previously described associations.

21Conversely, metabolic shifts may explain negative associations detected between other co-habiting microbes, e.g. Tannerella and Streptococcus in the subgingival plaque. The anaerobic Tannerella requires a much lower pO2 than Streptococcus and is proteolytic, while Streptococcus is a saccharolytic colonizer of the tooth surface that uses sugars as its primary source of carbon and is oxygen tolerant,. This continuous nutritional, metabolite (e.g. hydrogen peroxide), and oxygen gradient between the supragingval and the subgingival biofilms, along with differential exposure to host factors in saliva, is reflected through the gradual drop of the abundance of Tannerella as the streptococci increase (Figure S4). A similar example can be found in the Prevotella and species from the Flavobacteriaceae (represented here by Capnocytophaga; mean abundance 1.68±2.76%) in the tonsils. Less exposed surfaces of tonsillar crypts offer an anaerobic micro-environment favoring species like Prevotella, while other areas support the growth of carbon dioxide-dependant Capnocytophaga, a tradeoff that we detect here as a specific negative association.

22Phyla such as the TM7 and Synergistetes have only recently been characterized at the genetic level in the oral cavity,, and little is yet known about their roles in this microbial ecosystem. We identified a number of novel co-occurrences between members of these under-characterized phyla, including a positive association between members of the TM7 phylum (mean abundance 0.62±1.14%) and Moryella genus members (mean abundance 0.29±0.47%) in the tongue dorsum and a positive relationship between members of the Synergistetes phylum and Treponema in the subgingival biofilm. Since limited data on metabolic byproducts or requirements for these clades in the oral community are available, these newly identified putative interactors provide specific hypothesis for follow-up studies (e.g. by co-culture experiments).

23The degree to which microbial shedding from the oral cavity along the digestive tract might seed the distal commensal gut microbiota is as yet unclear. We found few (7) relationships between organisms in the two areas meeting our significance criteria, none of which were consistently supported by a majority of available data (Figure S5), suggesting no such direct microbial seeding within our level of detection in the healthy adult microbiome. Interactions detected within the gut itself consisted primarily of negative associations between Bacteriodes and Clostridia, especially members of the Ruminococcaceae family. These negative relationships reflect the tradeoff between Bacteroides (mean abundance 48.79±22.94%, range 1.47–97.14%) and Firmicutes (mean abundance 27.04±16.52%, range 1.49–91.78%), the two dominant gastrointestinal taxa and the subject of previous close study,. While oral microbial transit is clearly important during founding of the microbiome in infancy and in extreme cases such as illness,, these data suggest that it occurs at low levels in the normal adult microbiome. In such hosts, the naturally dense microflora of the lower gut may serve to further exclude the few bacteria that survive gastrointestinal transit.

24It is common practice to group microbial communities by ecological similarity,, and we extended this analysis method by summarizing relationships among similar habitats based on microbial cross-talk (Figure 5). Specifically, we organized pairs of body sites by the frequency with which they demonstrated co-occurring (or excluding) microbes (see Methods). Overall, this network recapitulates similarities in community structure among these microbial habitats as assessed by beta-diversity, with the added information of which microbes might drive these similarities. Conversely, co-exclusions spanning multiple habitats might represent cases in which competitive relationships or differing responses to host environment might bridge multiple habitats. Stool microbes (representing the gut microbiota), as above, did not demonstrate any detectable associations with inhabitants of the mouth; the airways microbiota (nares) likewise associated minimally with other body sites, although they were detectably structurally similar to the skin communities. The sub- and supra-gingival plaques were distinct from other mouth sites, and the vaginal communities and skin were again all highly similar. The sparsity of this body site network again illustrates that phylotypes rarely participate in detectable ecological relationships spanning distal body site habitats.

25We hypothesized based on previous findings in environmental communities that patterns of microbial co-occurrence and exclusion might be explained by their evolutionary relatedness and functional similarity. For example, closely related microbes might compete for limited resources, while functionally complementary bacteria would exhibit mutualism. To test this hypothesis, we compared two genomic properties of all microbial clades appearing in our network, their phylogenetic similarity (i.e. evolutionary relatedness) and a “functional” similarity score based on counting shared orthologous gene families (i.e. a measure of shared pathways and metabolic capacity). Phylogenetic distances were calculated as evolutionary divergence based directly on 16S sequence dissimilarity between all pairs of microbes. We compared this with a “functional” distance calculated as the Jaccard index of non-shared COG families between all pairs of microbial genomes (see Methods). For most pairs of microbes, these measures were highly correlated (Figure 6), not necessarily surprising in that both are influenced by gradual sequence change driven by molecular evolution.

26However, several exceptions to this pattern were apparent among the interacting organisms of our study. First, a dramatic separation of phylogenetic and functional distances occurred between positively and negatively associated clades (Figure 6, green lower left vs. red upper right): positive associations were enriched for both phylogenetic and functional similarity, while negative associations showed the inverse pattern. This was partially explained by the basic observation that similar organisms occupy similar niches, as most relationships among similar organisms occurred between clades at different body sites and often between the same clade at two proximal (e.g. oral) or bilateral sites (e.g. left and right retroauricular creases). Conversely, the preference for negative correlations to occur between phylogenetically and functionally different organisms (top right) suggests that the wide range of co-exclusion mechanisms, not only direct competition but also toxin production, environmental modification, and differential niche adaptation required substantial time to develop throughout evolution. Furthermore, interactions in the same body site were primarily negative, suggesting that competition or subniche differentiation were more prevalent in these data than were collaboration or niche sharing.

27Exceptions to both of these trends did occur, however, in that related organisms occasionally showed within-site competition, and phylogenetically distant clades sometimes co-occurred. A highlighted example of the former was the negative association between Bacteroides and Prevotellaceae family members (also phylum Bacteroidetes) in the gut, reflecting the recurrent tradeoff of this genus with the Prevotella as previously linked to enterotypes and/or dietary patterns. As these organisms are closely related, this might reflect alternative metabolic specializations in an otherwise fairly similar gut environment. Conversely, the Aggregatibacter were positively associated with members of the highly dissimilar Flavobacteriaceae family in the saliva. As mentioned above, the Capnocytophaga (dominant members of the Flavobacteriaceae in these data) are highly metabolically dependent, and positive correlations among organisms are enriched in oral biofilm associated organisms generally (see Figures 3 and 4).

28In addition to these on-diagonal outliers (Figure 6), several additional groups of organisms lay off the trend of functional and phylogenetic similarity. That is, some co-occurring/co-exclusive microbes were evolutionarily distant but functionally more similar than expected (below trend), while others were evolutionarily close but functionally distinct (above trend). Several relationships in the upper-left represent single functionally diverse clades that are also widely co-distributed among bilateral or related body sites, such as the Actinobacteria (skin) and Fusobacteriaceae (oral). Such clades' functional distances reflect a relatively high level of within-clade diversity, raising the possibility that a combination of environmental perturbation with a highly structured microenvironment might help to maintain a tension of high functional diversity within a limited phylogenetic range.

29In the oral cavity, co-occurrence of such outliers, appearing off the trend of functional and phylogenetic similarity (Figure 6) was limited to low-abundance community members, with some exceptions. Abundant signature taxa such as Streptococcus and Neisseria often excluded clades with more stringent environmental oxygen requirements regardless of their specific degree of relatedness. Prevotella - evolutionarily distant from the Neisseria but sharing much functional potential as defined by orthologous gene clusters - exhibited a negative association in the tonsil. Because of their functional similarity, particularly their shared metabolic requirements (both with varying degrees of saccharolytic and proteolytic activities at the species level, ), this strongly suggests a co-exclusion due to competition for resources, again in addition to their environmental oxygen requirement. At the family level, Pasteurellaceae (composed of Actinobacillus, Aggregatibacter, Haemophilus, and Pasteurella) negatively correlate with several other members of the microbiota. Co-occurrence is more common with pairs of functionally similar microbes able to co-exist through combinations of complementarity, commensalism, and cross-feeding of vitamins, amino acids, and other cofactors. Here, the Prevotella produce hydrogen, which influences the growth of Campylobacteria. Prevotella can also be supported by glycine and pyruvate produced from glutathione by Treponema species in the periodontal pocket. The most extreme case of organisms both related and correlated were the Bacteroides and Parabacteroides in the anterior nares. While these clades are not always well-resolved, this trend occurred in nine distinct samples (of 282 total), the only ones in which Bacteroides occurred nasally at >5% abundance. While a trivial explanation might be misclassification of a small portion of Bacteroides, this trend might instead suggest co-occurrence in a metabolic niche that rarely favors either organism but, in the rare occasion of favoring one, permits both.

30We analyzed ecological interactions among bacteria in the human microbiome using 16S marker gene abundance data from the Human Microbiome Project. Our methods for building a microbiome-wide microbial association network combined two complementary approaches: an ensemble of similarity/dissimilarity measures and a compendium of generalized boosted linear models. Relationship significance was assessed using a novel nonparametric approach to compositional data analysis, resulting in a network of co-occurrence and co-exclusion relationships representing potential microbial interactions and incompatibilities within and across body sites. Analysis of the network demonstrated strong organization of the human microbiota into body area niches, mostly among closely related individual body sites representing microbial habitats. A few “hub” microbes were observed to act as signature taxa driving the composition of each microcommunity. Many of these were also the dominant species within a body area, for example Streptococcus in the oral cavity and Bacteroides in the gut, and these highly abundant taxa also frequently co-occurred as connectors among multiple related body sites. In vivo mechanisms were available from prior work for many of these associations, and more generally the phylogenetic and functional relatedness of pairs of co-occurring microbes often explained their associations. In particular, phylogenetically related microbes tended to co-occur at proximal or environmentally similar body sites, while distantly related microbes with shared functional capacities tended to compete.

31This microbial association network was described from observational data, and the mechanisms underlying any of these putative interactions may be quite diverse. Positive co-occurrence association types could include nutritional cross-feeding, co-aggregation, co-colonization, signaling pathways, and co-survival in similar environments,. Negative exclusion interactions likewise might span toxin or small molecule production, environmental modification (to the detriment of microbial neighbors), immunomodulation, or gross overpopulation of a niche. Ecologically, these data alone do not resolve variations of mutualism, commensalism, amensalism, or predator-prey relationships,. Further, all of these ecological relationships, detected here based on microbial abundance patterns across many subjects, are themselves distinct from the biogeographical “co-occurrence” patterns observed by previous studies of individual microbes within subjects,. To distinguish between these, future work could include perturbation experiments (e.g. the removal of a species from a defined habitat such as the gut of a gnotobiotic mouse), as these are becoming less difficult to sustain technically. Analytic refinements might instead include defining directionality of relationships in higher-resolution (e.g. temporal) data; for instance, we expect a strict mutualistic relationship (where both partners cannot exist without the other) to be symmetric, whereas the relationship between a prey and a specialized predator is expected to be asymmetric (the prey can occur without its predator, but not vice versa). Negative co-exclusions may have fewer possible initial interpretations, comprising the types of competition outlined above, or they may indicate different, exclusive microbial community states occurring temporally or as linked to host environment,,.

32Methodologically, it is again important to emphasize that detecting significant co-occurrences among members of a population assayed as relative abundances can be surprisingly difficult due to compositionality. That is, an absolute increase in one organism's abundance can result in an apparent relative decrease of all other abundances, leading to spurious correlations. Extensive prior work has explained the problem in microbial and macroecological settings, and we have mitigated potential issues in these data through our ensemble approach and by principled calculation of significance thresholds using null distributions that incorporate the degree of similarity due solely to compositional effects (see Methods and Text S1). GBLMs were the most distinct method included in this ensemble and share some similarities with recently proposed genetic regulatory network (GRN) reconstruction techniques. GBLMs do provide methodology for discovering GRN-like higher-order interactions in microbial communities, but the accuracy needed to overcome the associated multiple hypothesis testing problems is not yet achievable from available 16S data,. We anticipate that future studies with species- or strain-level classification of deep shotgun metagenomic sequences may provide sufficient resolution for more detailed networks including such cooperative microbial associations.

33While it might be hoped that easily sampled microbiomes such as the saliva would serve as proxies for e.g. the broader oral microbiome, these results suggest that this is not generally the case. There are few strong correspondences among organisms even at closely related body sites, let alone distal sites, and very few cases where microbial abundance is quantitatively predictable from a proxy sample. In the HMP, this may be a feature of a healthy population, and additional relationships (or disruption of existing ones) might emerge in the presence of disease. Environmental factors that strongly impact the healthy microbiome may additionally not be captured for this population (e.g. diet) and can be further investigated by targeted methodology in future cohorts.

34This catalog of microbial co-occurrence and co-exclusion relationships thus provides an initial glimpse of potential mechanisms of community organization throughout the human microbiome. While this computational methodology can be applied to any communities assayed using marker gene sequencing, it is interesting to conclude by noting that the resolution of the resulting network is limited by the specificity of 16S sequence binning. The network discussed here, for example, leverages two specific hypervariable regions for taxonomic classification, each with strengths and weaknesses, and neither individually adequate for sequence classification at the species level,. Since it is likely that additional microbial associations will occur at the species or strain level, we anticipate that further community structure will emerge during analysis of metagenomic shotgun sequences taxonomically binned at a finer level of detail. Community shotgun sequencing will also provide functional information regarding metabolism, signaling, and, again, potential physical mechanisms of interaction, which can in turn be matched against complete reference genomes for co-occurring strains. Perturbation analyses in co-culture or, eventually, longitudinal studies in human cohorts will provide an intriguing means of investigating the impact of these microbial “wiring” diagrams on human health.

35Two complementary approaches, an ensemble of multiple similarity/dissimilarity measures and a compendium of generalized boosted linear models (GBLMs), were used to interrogate significant associations between microbial abundances. These were drawn from 18 body sites assayed by the Human Microbiome Project at two clinical centers using 16S rRNA gene sequencing. Simes method and Benjamini-Hochberg-Yekutieli false discovery rate (FDR) correction were used to combine the resulting networks. From this merged, global network, we summarized overall network properties (degree distribution, modularity, etc.), assessed patterns of microbial connectivity within and among body sites, and identified highly connected (hub) microbes. Phylogenetic and functional distances were calculated based on 16S rRNA gene sequence similarity and shared orthologous gene families, respectively, and combined with the network.

36Please see Text S1 for an extensive discussion of the methodology used to assess relationship significance in compositional data, which is presented in summary below.

37The 16S rRNA gene-based dataset of the normal (healthy) human microbiome was made available through the Human Microbiome Project (HMP) and is detailed in. Briefly, it consists of 454 FLX Titanium sequences spanning the V1 to V3 and V3 to V5 variable regions obtained for 239 healthy subjects enrolled at clinical sites in Houston, TX and St. Louis, MO. These cover 18 body sites covering five areas: the oral cavity (nine sites: saliva, tongue dorsum, palatine tonsils, keratinized gingiva, hard palate, buccal mucosa, throat, and sub- and supragingival plaques), the gut (one site: stool), the vagina (three sites: introitus, mid-vagina, and posterior fornix), the nasal cavity (one sample: anterior nares), and the skin (four sites: left and right antecubital fossae and retroauricular creases). Sequences of both 16S windows were processed separately using mothur into phylotypes using the RDP taxonomy as described in and, with full protocols also available on the HMP DACC website (http://hmpdacc.org/HMMCP). Genus level and above phylotypes were used for this analysis, for which the datasets from both windows were combined.

38This resulted in more than 5,000 samples comprising 910 taxa made available as part of the HMP (http://hmpdacc.org/HMMCP). These were further processed for this study by excluding any phylotype not supported by at least two sequences in at least two samples. Samples were removed as suspect if the most abundant taxon was detected by fewer than 1% of the sequences supporting it in the sample in which it was most abundant, and counts for the remaining 726 taxa were converted to relative abundances in each of the resulting 5,026 samples. Due to potential differences between clinical centers, the dataset was conservatively split into two subsets for further analysis, subjects recruited in Houston (3,296 samples) and those recruited in St. Louis (1,730 samples).

39For each resulting dataset, a compendium of generalized boosted linear models (GBLMs) was constructed by selecting all 324 combinations of source body sites ss and target sites ts. Each GBLM was fit using the abundances of all source taxa st within the source site to predict the abundance of each target taxon tt within the target site using a sparse linear model of the form:

40All additional non-leaf clades in the RDP taxonomy (i.e. families, orders, etc. up to the bacterial and archaeal domains) were included as source and target taxa. For ss equal to ts, i.e. predicting the abundance of a taxon tt when the abundances of all taxa in the same body site are known, the abundances of all parent and descendant clades of tt were excluded from the available source taxa st. That is, when ss = ts and tt was of the form domain|phylum|class|…|clade, all source taxa of the form domain, domain|phylum, domain|phylum|class, etc. or of the form domain|phylum|class|…|clade|subclade, domain|phylum|class|…|clade|subclade|subsubclade, etc. were excluded from the source taxa st. This prevented the abundances of xtt,ts from being predicted using abundances xst,ss on which they were directly dependent, while allowing the detection of predictive relationships between distinct clades within the same body site.

41The linear model was generalized to include binary categorical target taxa (in this case only gender and ethnicity) using standard logistic regression:

42As this is clearly an extremely high-dimensional problem, multiple a priori and post hoc steps were taken to enforce model sparsity and to avoid overfitting for each (ss, ts, tt) tuple. The first of these was to exclude from the available st any taxon not correlating with tt at a Spearman correlation of nominally p<0.05. The second was to boost linear model fitting rather than attempt to fit all βtt,ts,st,ss simultaneously. Boosted linear models retain the usual L2 least squares penalty, but are constructed in a manner similar to sparse forward variable selection or the LASSO. βs are considered for inclusion in the model one at a time and the parameter minimizing sum of squared error selected and included. However, each subsequent round of parameter fitting operates on the residuals of all previous rounds, thus “upweighting” poorly fit examples, and the inclusion of further non-zero βs stops after a fixed number of iterations.

43This tuning parameter and the model fitting process was 10-fold cross-validated and selected from the most accurate (by root mean square error for continuous tt and AUC for binary) of 50, 100, or 150 boosting iterations using the caret and mboost R packages. This resulted in a compendium of 7±4.9 non-zero parameters β, each evaluated with a 10× cross-validated R2/AUC and a nominal R2/AUC on the full dataset. Final model quality scores were assigned by A) subtracting AUCs below 0.5 from one, since caret does not calculate AUCs directionally, and B) retaining the minimum of the cross-validated and nominal R2/AUC. Any continuous model not achieving an R2 above zero after adjustment for the number of non-zero parameters, parameters p, training samples n) was discarded (55,424 retained).

44Even from cross-validated goodness-of-fit scores, the compositional structure of relative abundance data prevents straightforward assessment of model significance (see Text S1). We thus additionally fit twenty models per (ss, ts, tt) tuple after bootstrap re-sampling the values of tt across samples. The AR2/AUC of these bootstrap models provided a confidence interval around the observed AR2/AUC. Any model for which the 90% confidence interval failed to include the observed AR2/AUC was discarded. To construct the null distribution of associations due to compositionality alone, we fit twenty additional models after permuting the values of tt across samples and renormalizing them sample-wise, thus retaining compositional effects but breaking true associations. The GBLM was re-fit and the resulting null distribution of AR2/AUC values used to assess significance of the true model. The mean and standard deviation of the bootstrap distribution was z-tested against this null distribution using the jointly pooled standard deviation, providing one p-value per model per clinical center. Any model with FDR adjusted p-value greater than 0.05 was discarded (18,286 retained).

45The 16S data described above were first normalized by dividing each sample by its total phylotype sum. Mislabeled samples were removed and samples were again processed as two subsets, one per clinical center (Houston and St. Louis). These were encoded as a matrix in which each row represented a phylotype in a specific body site and each column represented an individual during one sampling visit. Rows with more than 2/3 zero counts were removed, leaving matrices of 1,217 (Houston) and 1,408 (St. Louis) phylotype-bodysite composite features collected for 248 and 144 subject-visit points, respectively.

46We built on composite co-occurrence scores as described in, for example, (for protein functions) and (graph clustering) to find groups of microbial lineages co-existing across a large number of environments. Specifically, we combined four diverse measures in order to overcome two major challenges in the inference of co-occurrence networks, particularly appropriateness of scoring measures to sparse count data and determination of statistical significance. The first is exemplified by the double-zero problem, in which a zero indicates either that an organism is absent or that it is below detection limit; Pearson and Spearman correlations are sensitive to this, while the Bray Curtis dissimilarity is not. The latter issue arises from the need to normalize across samples with unknown absolute abundances, either by relativizing or by downsampling; either procedure results in constrained sample sums, which introduce artificial correlations.

47We thus employed an ensemble approach combining four diverse measures: two measures of correlation (Pearson and Spearman) and two measures of dissimilarity (Bray-Curtis (BC) and Kullback-Leibler (KLD)). For BC and KLD calculations, rows were divided by their sum prior to computation. Additional measures were considered for our ensemble, including the Hellinger, Euclidean and variance of log-ratios, but these proved to be well-represented by the smaller final ensemble (see Figure S6).

48After running each of the above measures on the two 16S data matrices, one per clinical center, we set measure-specific thresholds as a pre-filter such that each measure contributed 1,000 top-ranking and 1,000 bottom-ranking edges to the network. Edge scores were computed only between clade pairs without parent-descendant relationship (e.g. without pairs of the type Actinobacteridae|Actinomycetales or Actinomycetales|Propionibacterineae) for clades in the same body site. To assign statistical significance to the resulting differently-scaled scores, we first computed edge- and measure-specific permutation and bootstrap score distributions with 1,000 iterations each. In order to address the compositionality issues discussed above, we re-normalized the data in each permutation, providing a null distribution that captures the similarity introduced by compositionality alone (see Text S1). We then computed the p-value as above by z-scoring the permuted null and bootstrap confidence interval using pooled variance. P-values were tail-adjusted so that low p-values correspond to co-presence and high p-values to exclusion. For BC and KLD, we did not compute re-normalized permutations, because these measures are intrinsically robust to compositionality. Instead, we calculated their p-values using the bootstrap interval compared to a point null value that was computed by permutation.

49Finally, to remove unstable edges, we removed all edges whose score was not within the 95% confidence interval (limited by the 2.5 and 97.5 percentiles) of the bootstrap distribution. Additionally, a number of BC-supported negative links were removed because they were due to abundance profiles including one extreme outlier. This affected the following clades for St. Louis: Actinomycetales in stool, Corynebacterium and Corynebacteriaceae in the tonsils, Lactobacillus and Lactobacillaceae in the anterior nares and for Houston: an unclassified Neisseria in the left retroauricular crease.

50Data summarized as relative abundances are referred to as compositional. Because they sum to one, their elements are not independent and may exhibit spurious correlations regardless of the true underlying relationship. To mitigate this effect, we assessed the significance of our ensemble scores and GBLMs as described above, each by comparing a bootstrap confidence interval around the observed score with a permuted null distribution that includes repeated renormalization to account for compositional effects alone. In each permutation, only the target taxon row was randomized and all samples subsequently renormalized to a constant sum. Because permutation breaks correlation structure while renormalization reintroduces compositionality, a null distribution coupling these elements induces correlation from compositional structure alone. Comparing this null distribution to a standard bootstrap confidence interval around the observed value provided a straightforward nonparametric test of association accounting for compositionality. Simulation studies showed that the bootstrap-renormalization scheme was successful in discounting compositional effects while preserving true correlations (see Text S1).

51Simes method was used to combine all ten networks (5 methods×2 study centers) into one final network, as it is robust against non-independent tests. A strict intersection of the two clinical centers' networks, rather than a p-value combination, was also examined and found to be over-stringent due to systematic differences in the data (Figure S8). After merging, p-values on each final edge were corrected to FDR q-values using the Benjamini-Hochberg-Yekutieli method and a q-value cutoff of 0.05 was applied. The positivity or negativity of each relationship was determined by consensus voting over all integrated data sources, ranging from −10 (most negative) to 10 (most positive, see Figure 2). Edges with indeterminate directionality (direction score of zero) were removed. Finally, only edges with at least two (out of ten) supporting pieces of evidence were retained.

52The formula by Clauset et al. compares the fraction of edges within input clusters with the fraction of within-cluster edges that would be expected for a randomized network. We clustered the network using the Markov cluster algorithm (MCL) and computed network modularities for a range of inflation parameters. The strongest modularity (0.28) was measured for an inflation of 1.3 and is slightly below the cut-off recommended by Clauset et al.. This modularity was, however, higher than any measured for 100 randomized networks (which preserved node and edge number, but in which edges were randomly re-assigned) and was therefore retained as significant.

53To assess whether specific body sites were more connected than expected by chance, we repeatedly (1,000 times) selected as many nodes as a body site contains from the global network at random and counted their edge number. This resulted in a distribution of edge numbers for random node sets. To retain only plausibly significant edges for further calculation, we corrected for multiple testing by multiplying the nominal p-value from this distribution with the number of tests carried out and retaining values below 0.05. We repeated this test separately for positive, negative, intra- and cross-edges. For visualization, the network itself was also separated into within-site, within-area, and between-area subsets for further inspection (Figure S9).

54To assess the interaction strength between body site and clade groups (Figures 4 and 5), the number of relationships between nodes in each group pair was counted and normalized by the group member product, which represents the number of links two groups can potentially form. This was repeated in 1,000 randomized networks (generated as described for the computation of network modularity). A p-value was computed from the count distribution and node group relationships with p-values above 0.05 were discarded, retaining only the fractions of total possible relationships which were significantly higher than those expected by chance.

55Genome sequences of 1,107 organisms were retrieved from NCBI (http://www.ncbi.nlm.nih.gov/genomes/lproks.cgi, December 2010) and 16S sequences were extracted. These 16S sequences were aligned using MUSCLE 3.8.31 and a full phylogenetic tree reconstructed by maximum likelihood using FastTree 2.1. A matrix of all pair-wise distances was created from this cladogram and distances between any two nodes (e.g. families, orders, etc.) calculated by taking the median of all distances (as provided in units from FastTree) between all pairs of leaf taxa descending from the two nodes.

56Functional complements of the same genomes were summarized using COG families as assigned by NCBI annotations. This resulted in an abundance matrix with 4,685 columns (corresponding to COG families) and one row per genome. Columns summing to less than 10% of the number of genomes were removed, resulting in 3,514 usable COG families. Pairwise scores between genomes were calculated using Jaccard index, with distances for higher-level clades computed using medians in the same manner as described above. Final functional distances were represented as the Jaccard index of non-shared COG families between pairs of genomes.