Detection of stable community structures within gut microbiota co-occurrence networks from different human populations

Data aggregation

Given our objective to compare network community structures across datasets (Fig. 1A), we required OTUs that would be comparable between them. To this end, we carried out clustering of sequences across combined data from multiple sources. To maximise sequencing similarity, we selected two datasets with experimental approaches best matching those of the gut microbiota profiles obtained for TwinsUK, which used 16S rRNA gene sequencing of faecal samples. The LLDEEP and Israeli-PN datasets were selected as they carried out gut microbiota profiling by amplifying the V4 region of the 16S rRNA gene using the same PCR primers, and used paired-end sequencing on the Illumina MiSeq platform with read lengths sufficient to capture the whole V4 region. Notable differences between studies include faecal sampling and DNA extraction techniques. Both TwinsUK and LLDEEP utilised aliquots from faecal samples stored at −80 °C, whilst the Israeli-PN study utilised a mixture of faecal swabs stored at −80 °C and OMNIgene-GUT stool collection kits stored at −20 °C. All three studies used both chemical and mechanical lysis in DNA extraction but employed different protocols: TwinsUK utilised the MoBio PowerSoil HTP extraction kit, the LLDEEP cohort utilised the Qiagen AllPrep kit, and the Israeli-PN DNA was extracted using the MoBio PowerMag Soil DNA extraction kit.

Extraction of DNA from faecal samples, amplification and sequencing of the V4 region of the 16S rRNA gene, and demultiplexing of sequencing reads has previously been described for the TwinsUK cohort (Goodrich et al., 2014). The paired-end demultiplexed reads were joined using join_paired_ends within QIIME with an overlap of at least 200nt to form single reads covering the full V4 region (Navas-Molina et al., 2013). DNA extraction and 16S rRNA gene sequencing within the LLDEEP samples has been described in detail previously (Fu et al., 2015). Data from the LLDEEP cohort was provided in a similar format having used custom scripts to merge the paired end data to full length reads covering the V4 region and split data by individual (Gevers et al., 2014; Fu et al., 2015). Ethics approval for the TwinsUK study was given by the NRES Committee London—Westminster (REC Reference No. : EC04/015), and the Lifelines-DEEP study was approved by the institutional ethics review boards of the University Medical Center Groningen (ref. M12.113965).

Raw data accompanying the Israeli-PN publication was downloaded from the European Nucleotide Archive (ENA) (Accession: PRJEB11532) (Zeevi et al., 2015). This was processed similarly to the TwinsUK data and demultiplexed using published barcode mappings. Samples not listed in the accompanying metadata or with ambiguous read identifiers were removed.

OTU clustering, table filtering, and diversity analyses

The cleaned, joined, and demultiplexed data was concatenated to produce one sequencing file covering all three studies. After quality control, 4,426 samples (2,764 TwinsUK, 1,023 LLDEEP, 639 Israeli-PN) were included in the analysis. The complete sequencing from all three sets contained 381,767,528 reads. These were dereplicated, removing reads occurring only once, resulting in 5,728,288 unique reads. We chose to use a de novo approach to cluster these to OTUs, as de novo OTU clustering is not influenced by the reference database used, can capture novel diversity, and can produce more accurate clustering of OTUs than reference-based approaches (Westcott & Schloss, 2015; Jackson et al., 2016a). De novo clustering was carried out using the cluster_fast command in the VSEARCH package with a 97% similarity threshold (Rognes et al., 2016). The resultant 94,070 representative sequences were filtered to remove chimeric reads using the uchime_denovo command (within VSEARCH), producing a final set of 17,123 representative sequences. These were used to generate OTU counts across all samples using the VSEARCH usearch_global command (Rognes et al., 2016). Post-processing, average read counts were 82,695 ± 745 for TwinsUK, 49,962 ± 964 for LLDEEP, and 130,378 ± 5,534 for Israeli-PN (mean ± SEM).

A phylogenetic tree was generated from the representative sequences using the default parameters of the make_phylogeny command in QIIME (Navas-Molina et al., 2013). Taxonomy of OTUs was assigned by matching representative sequences against the Greengenes v13_8 database using the default parameters of the assign_taxonomy command in QIIME (Navas-Molina et al., 2013). OTUs occurring in only one sample, and samples with less than 10,000 OTU counts were removed. Weighted and unweighted UniFrac beta diversity measures and subsequent principle co-ordinates analysis of them was carried out using the beta_diversity_through_plots script in QIIME (Navas-Molina et al., 2013). The combined OTU table was then split by dataset. For the purposes of alpha diversity calculations, the raw counts tables were rarefied to a depth of 10,000 reads. For each sample, Shannon and inverse Simpson diversity indices were calculated as the mean across ten rarefactions. Significant differences in alpha diversity between datasets were assessed using Mann–Whitney U-tests.

Calculating co-occurrence measures

Co-occurrence calculation was carried out on each dataset independently. Sub-tables were generated from raw (unrarefied) OTU tables that only contained OTUs found in at least 25% of the samples (Jackson et al., 2016b). All co-occurrence measures were calculated within these subsets as they are less sparse and hence more amenable to correlation measures (Weiss et al., 2016). The mean diversity was assessed on rarefied versions of the subset tables using the inverse Simpson index as for the full tables. OTU sparsity was assessed from the unrarefied table using the biom summarize_table command. Following the recommendations of Weiss et al. (2016), we then used these estimates of the mean inverse Simpson index (TwinsUK = 20.2, LLDEEP = 29.0, and Israeli-PN = 13.1) and OTU table sparsity (0.49 in all three) to select an ensemble approach to co-occurrence detection that combines four different correlation measures: CoNet, SparCC, Pearson’s, and Spearman’s.

P-value thresholding to generate co-occurrence networks

Intersected networks were generated by combining the edge tables from the CoNet, SparCC, Pearson’s, and Spearman’s methods. This was done independently for each dataset. Edges were retained in the final network if the direction of co-occurrence matched, and the edge p-values were below a given threshold in all four methods. This was carried out at multiple different p-value thresholds (0.05, then ranging in powers of ten from 0.01 to 10⁻⁸), to generate multiple intersect networks for each dataset with a gradient of stringency for edge inclusion. We then aimed to determine which was the most appropriate threshold to use to generate the final networks. Rather than make an entirely arbitrary selection, we chose to use fit to a scale-free network as a biologically motivated method to identify the optimum p-value threshold to use.

A scale-free network has a node degree distribution that follows a power law, i.e., there are few highly connected nodes and many more less connected nodes. This distribution has been observed in several biological networks, including microbiota co-occurrence networks (Jeong et al., 2000; Albert, 2005; Zhang & Horvath, 2005; Faust et al., 2012; Tong et al., 2013). It is also frequently used as a threshold for edge inclusion in the widely used WGCNA package (Langfelder & Horvath, 2008). Fit to a scale free network was calculated by first extracting the degree distribution of a network using igraph (Csardi & Nepusz, 2006). The fit of this distribution was then assessed using the scaleFreeFitIndex function from the WGCNA package in R (Langfelder & Horvath, 2008). This provides the R² of the fit to a scale-free model, which the creators of WGCNA suggest should be >0.8, and the slope of the fit, which they suggest should be close to −2 to indicate a good fit. The optimum p-value threshold was selected based on these criteria across all three datasets (see Results) and the resultant intersect network at this threshold was used in all further analyses. Visualisation and generation of descriptive statistics from the final networks was carried out using Gephi (Bastian, Heymann & Jacomy, 2009).

Detecting communities within co-occurrence networks

Between OTU adjacency matrices were generated that represented the final intersect networks for each dataset. Negative correlations were removed (considered zeroes in the adjacency matrix) to generate unsigned networks as they represented a considerably small proportion of edges (<1% in all datasets). Community detection was carried out on each dataset’s network independently using the genLouvain 2.0 package within MATLAB (Mucha et al., 2010), which implements the Louvain approach to modularity maximisation (Blondel et al., 2008). This defines community partitions by assigning nodes to unique communities, then iteratively combining neighbouring nodes into communities if it results in an increase in modularity across the whole network. Modularity is a measure of the number of edges within communities relative to between communities and a higher value represents better community definition (Newman & Girvan, 2004).

The genLouvain algorithm includes a γ parameter, which controls the size and number of detected communities (Mucha et al., 2010). A smaller gamma value promotes the detection of a small number of larger communities, while larger gamma values promote the detection of a high number of smaller communities. To find an optimal value for the γ parameter, we carried out community detection using a range of γ values for each dataset (0.1–1, increments of 0.01), and assessed the stability and statistical significance of the resultant communities.

Community properties and host associations in TwinsUK

The relative abundance of each community in a sample was found by summing the read counts of its constituent OTUs and dividing by the total number of reads observed in the sample. Association between the mean abundance of an OTU in a dataset and the number of OTUs in its parent community was assessed using Spearman’s correlation. Taxonomy within communities was investigated by counting the number of OTUs assigned to different taxa at each taxonomic level. The mean identity between the representative sequences of the OTUs in each community was assessed for each dataset. All pairwise comparisons between the OTUs within a community were carried out using global alignments performed using BLAST within the pairwise2 command of Biopython, with a score of 1 for a match, −1 for a mismatch, −2 for opening a gap, and −1 for extending a gap (the defaults used by the Mothur align.seqs command) (Cock et al., 2009; Schloss et al., 2009).

The heritability of TwinsUK communities was estimated using the mets package in R. Log transformed relative abundances were used for each community. For each community ACE, CE, E, and AE models were fit using data from complete twin pairs. Co-variates included age, BMI, sample collection method, sex, and sequencing depth. For each community, the best fitting model was determined as that with the lowest Akaike information criterion.

Association analyses between community abundances and BMI and age was carried out using log transformed relative abundances of the communities (log10 (relative abundance +10⁻⁶)). These phenotypes were selected as they were also available in both the LLDEEP and Israeli-PN datasets. To investigate associations, linear models were fitted for each community using the lm function in R. Community abundance was the dependent variable with BMI, age, gender, and sequencing depth as independent variables, as these were the maximal set available across all datasets. The coefficient and significance of associations were extracted for BMI and age. P-values were then adjusted for multiple testing using the FDR method in R.

Comparison of community structure between populations

Integrating microbiota data from different populations

The aim of this study was to compare gut microbiota community structures across populations (study overview in Fig. 1A). To this end comparable OTUs were generated by combining gut microbiota sequencing data from three geographically diverse populations. These were the British TwinsUK and Dutch Lifelines-DEEP (LLDEEP) cohorts and data from an Israeli personalised nutrition study (Israeli-PN) (Table 1) (Fu et al., 2015; Zeevi et al., 2015). Principle coordinates analysis of overall microbiota composition between samples found that whilst there was some grouping by dataset, there was also significant overlap between them (Fig. 1B). Comparison of median alpha diversities between the populations found significant differences between all three (Fig. 1C).

Across the 15,361 OTUs detected in the complete data 48% were shared across all three datasets, with 79% being found in at least two datasets (Fig. 1D). There was also considerable overlap between TwinsUK and LLDEEP of OTUs not found in the Israeli-PN data. Similar patterns were observed when considering more abundant OTUs (found in >25% of the respective populations). Although LLDEEP and TwinsUK shared more OTUs, examining the mean taxonomic distributions at the phyla level the TwinsUK and Israeli-PN cohorts were most alike (Fig. 1E). The LLDEEP cohort contained relatively lower levels of the phylum Bacteroidetes, and a higher abundance of Firmicutes bacteria. Further differences were observed at the genus level, where the Israeli-PN study had a higher average abundance of Prevotella.

Scale-free p-value thresholding of edges in co-occurrence networks

We split the OTU data and generated co-occurrence networks for each dataset using an ensemble approach—taking the intersect of four different correlation measures (TwinsUK example in Fig. 2A). This required selection of a p-value threshold at which correlations were considered significant and retained as edges in the networks. We generated ensemble networks for each dataset at a range of p-value thresholds and selected a cut-off where the resulting networks’ degree distributions best fit a scale-free distribution (Fig. 2B). We observed no consistent trend across the p-values, but this might be expected given the differences in the number of OTUs and samples in the datasets and varying levels of p-value precision between the co-occurrence measures used. Indeed, investigating network properties at each threshold for each method individually (Fig. S1), we observed that the p-value accuracy of SparCC meant that at at all thresholds below 10⁻² there was no change in the number of edges in the SparCC networks (only edges with p = 0 remained). We also observed that the drop in fit to a scale-free distribution at lower p-value cut-offs in the TwinsUK and Israeli-PN networks was driven by trends in the CoNet networks, which had the lowest number of edges and hence were the principal determinant in the intersects at these lower thresholds. However, overall, the trends of the ensemble networks did not simply reflect those of any one of their constituent methods, displaying properties only emergent upon intersection of all four. From examination of the ensemble networks (Fig. 2B), we found that including edges with a p-value <0.01 produced intersect networks with a good fit to the scale-free model in all three datasets. As such, this threshold was selected to generate the networks used for community detection.

From the 6,761 unique edges observed across in all three networks, 166 were observed in all three, 882 in the TwinsUK and LLDEEP, 677 in TwinsUK and Israeli-PN, and 229 in the LLDEEP and Israeli-PN datasets. Summary statistics for the resultant networks are shown in Table 2. The TwinsUK and LLDEEP network structures were most alike. The Israeli-PN data, which contained fewer nodes (OTUs), produced a smaller and more connected network.

Detection of communities in microbial co-occurrence networks

Selecting a γ parameter for community detection

After establishing co-occurrence networks we used the Louvain modularity maximisation algorithm to detect communities within them (Mucha et al., 2010; Blondel et al., 2008; Newman & Girvan, 2004). The version used includes a constant (γ) that can be used to manipulate the number and size of resultant communities (see ‘Materials and Methods’). To determine an appropriate value for γ, we carried out repeated community detection on the three co-occurrence networks at various γ values and assessed the stability (as variation of information between runs) (Fig. 3A) and statistical significance (mean modularity of real compared to randomised networks) (Fig. 3B) of the community definitions at each γ.

Community definitions were stable (low variation of information between runs) at γ values approaching zero (Fig. 3A). This would be expected as in this instance most OTUs fall into a few large communities (Fig. 3C). Similarly, at these low γ values, a networks edges are most likely to be within communities, hence we observed high modularity values for both the true and randomised networks. As γ increased we found modularity estimates for the random networks dropped significantly lower than those of the true networks (Fig. 3B). This shows that the real co-occurrence networks have significantly more community structure than would be expected by chance, and provides confidence that the communities identified at these γ values reflect the biological relationships of their member OTUs. The stability of definitions (mean variation of information) fluctuated across the γ range, however it should be noted that the lowest estimates (1—variation of information) were above 0.9 (range 0–1, least to most similar), showing that even the least stable γ produced very similar communities between runs.

Final community definitions

From Fig. 3, we selected a γ value of 0.4 for community detection. This γ provided both high modularity estimates (that were significantly higher than the random networks) and good stability in all three datasets. We then used the Louvain algorithm with γ = 0.4 to detect the final OTU community definitions in the three co-occurrence networks (Fig. 4A, Table S1). Communities in each network are hence referred to using arbitrary numbers ranging from 1 to the number of communities in the network.

The number of communities detected within the networks (TwinsUK = 96, LLDEEP = 105, Israeli-PN = 31) reflected the number of OTUs within them, but was more similar when only considering communities containing at least five OTUs (TwinsUK = 35, LLDEEP = 36, Israeli-PN = 11). There was a positive correlation between the mean relative abundance of an OTU and the total number of OTUs in the community it was assigned to in the TwinsUK and LLDEEP networks (rs = 0.1, p = 0.002 and rs = 0.27, p < 0.0001 respectively). However, overall each community in all three networks contained a range of both high and low abundance OTUs (Fig. S2).

Taxonomy, heritability, and host associations of communities in TwinsUK

Comparison of communities across populations

Similarities in overall community structure

Having identified communities within each dataset we aimed to determine how similar the segmentation of OTUs was between them. We again used the normalised variation of information measure to compare groupings, and found that the similarity of community assignments was significantly higher in all three pair-wise comparisons than would be expected by chance (Fig. 5A). This shows that OTUs are forming similar communities across the three geographically diverse datasets.

We also repeated community structure comparisons between networks using both higher and lower p-values to threshold edges, and higher and lower values for γ in community detection (Fig. S5 and Table S5). In all instances the community structure was higher than expected by chance, and similar communities were observed both within and between datasets when using alternate parameters. This shows that the approach is robust to differences in parameterisation, and that the similarity observed between datasets is not solely manifest at the selected values. However, it should be noted that in all comparisons the variation of information can only be used to compare matching sets, so these results only show communities are similar within the OTUs shared by networks.

Differences between datasets

In combining these datasets we observed significant differences in the diversity and taxonomic profiles between datasets. These likely reflect a mixture of both the differences in the study populations and the experimental approaches used. Environmental differences between cohorts, in terms of genetics and other lifestyle factors are known to influence the gut microbiome (Goodrich et al., 2014; Falony et al., 2016; David et al., 2014). For instance, low diversity Prevotella dominant microbiomes have been previously observed in cohort studies at differential levels and are potentially linked to dietary intake (Arumugam et al., 2011; Wu et al., 2011). The larger proportion of Prevotella dominance, and associated reduced diversity, in the Israeli-PN gut microbiota might therefore reflect geographical and dietary differences. However, technical disparities such as faecal sampling and extraction method, which were known to be different between the present studies, are also known to influence taxonomic composition and cannot be discounted (Walker et al., 2015; Kennedy et al., 2014; Sinha et al., 2017). The increased overlap in the TwinsUK and LLDEEP data in comparison to the Israeli-PN dataset in both the OTU and network structures might also result from the higher sample sizes in these studies. However, it is notable that even with the intrinsic population and technical differences between the three studies we still observed consistent elements across them. Most OTUs were found in at least two studies and there was overlap of samples from all three in beta diversity PCoA plots. It is within this common ground that similarities in community structure were observed.

Generating co-occurrence networks

To generate the co-occurrence networks in this study we used an ensemble approach as suggested by Weiss et al. (2016). However, we found that CoNet was the main edge filter when intersecting networks and also produced the best fit to a scale-free distribution when considered alone, and thus might be sufficient. This might be expected given that it is itself an ensemble approach. However, metrics such as SparCC are not implemented within CoNet, and different correlation methods can produce better results depending on the format and properties of the input dataset (Weiss et al., 2016). These observations may therefore only apply to the current study. Furthermore, there are several alternative approaches that were not considered in the comparisons of Weiss and Van Treuren et al. such as SPIEC-EASI and CCLasso (Fang et al., 2015; Kurtz et al., 2015). Further exploration of community detection using these methods is warranted. An approach to generate quantitative microbiota profiles using a sample’s total microbial cell count has also recently been described (Vandeputte et al., 2017). This negates issues of data compostionality and was shown to reduce the number of spurious correlations observed between microbiota. Using such approaches in future studies could enable the use of more traditional correlation metrics and simplify downstream community detection.

We selected edges for inclusion in co-occurrence networks based on their fit to a scale-free topology. Comparing community structures using different p-value thresholds we observed similar community structures between different datasets when using different thresholds. This suggests different p-value cut-offs could be used for each dataset if a single threshold does not produce a reasonable scale-free fit across all datasets as in this study. However, we also observed similar community structures within datasets when using thresholds producing networks with a poorer fit to a scale-free distribution. This is likely due to the use of a more stringent ensemble approach to co-occurrence estimation, which limits edges in the network to the strongest and most consistent associations (arguably most likely to reflect biological phenomena) prior to further thresholding. Use of a scale-free approach to p-value parameterisation might therefore be unnecessary where similar stringent co-occurrence methods are used. It is also not possible to determine if the true network has a scale-free distribution, and indeed the true network could be defined in a number of different ways depending on what relationships are considered to constitute an interaction or edge.

However, it has previously been shown that microbial interaction network structures can approximate a scale-free distribution (Chaffron et al., 2010; Faust et al., 2012; Tong et al., 2013). Scale-free networks are also well conserved in other aspects of biology across different domains of life (Wolf, Karev & Koonin, 2002; Albert, 2005), and this approach is used by existing methods such as WGCNA (Langfelder & Horvath, 2008). Using typical co-occurrence approaches a p-value cut-off will be required to determine edges for inclusion when generating a network. Using a scale-free fit, or similar data driven approach, avoids using entirely arbitrary selections and provides a uniform rationale that can be used to generate comparable network structures across different datasets.

Community detection

Following network creation, we utilised a modularity maximisation algorithm to detect communities within the network. Comparing modularity to randomised networks ensured we were detecting structure manifest as a result of non-stochastic processes. Similarly, our application of variation of information enabled identification of the most stable communities across the range of γ values tested. Such parametrization steps are often neglected in network analyses but should be considered in future microbiota studies, particularly where the aim is to identify the optimal community configuration. However, overall the community definitions were stable and similar across the range of γ values considered here, suggesting that the resultant OTU communities were robust to small changes in the γ parameter.

Modularity maximisation assigned OTUs to single communities. This may not represent the true nature of microbiota interactions but provides some analytical benefits. Unambiguous definition of communities facilitates comparison between network structures. This has applications in comparison between disease and control groups (Baldassano & Bassett, 2016) or, as we have shown, replication of community associations across datasets. Unambiguous assignments of OTUs to individual communities might also be desirable for studies aiming to design synthetic communities to replicate the benefits observed with treatments such as faecal microbiota transplants (De Roy et al., 2014). For example, these methods could be applied to identify common microbial communities in successful donor samples, and generate more controlled, rationally designed synthetic communities for use in place of current less specific approaches. There is also evidence that a community-centric approach to studying the gut microbiome in relation to human health can be more relevant than investigating individual taxa. Several studies have shown that associations between individual gut microbiota and host health can depend on their wider community context (Gevers et al., 2014; Goodrich et al., 2014; Baldassano & Bassett, 2016; Ridaura et al., 2013).

Beyond identifying communities, understanding their formation would require extension of the described approaches to quantification of taxa using metagenomic sequencing. This, especially used in tandem with metatranscriptomics and metabolomics, might provide indications to the mechanisms of interaction. Metabolic modeling from metagenomic data can also be used to improve the inference of interactions from co-occurrence observations by predicting interspecies metabolic dependencies and/or shared niche specialization (Levy & Borenstein, 2013). However, even with such improvements, cross-sectional approaches are inherently limited to inference of interactions from co-occurrence across samples. Time-series data and in vitro studies will also be required to delineate directional effects and validate individual interactions (Faust et al., 2015). For instance, there are existing methods that have been used to infer interactions from covariation between microbiota across time series (Steele et al., 2011), and it has been shown that interactions observed in pair-wise species co-cultures can be used to predict outcomes in more complex multi-species cultures (Friedman, Higgins & Gore, 2017). Combining in vitro observations with community interactions observed in the host environment could be a powerful tool for the design of synthetic bacterial communities for use as gut microbiome targeting medicines (Lindemann et al., 2016).

The stability of communities across populations

Although 16S rRNA gene data cannot elucidate interaction mechanisms, we were able to determine several biological phenomena from the microbial approximations provided by OTUs. Most notably, that OTUs formed similar communities with similar host associations in the co-occurrence networks of three different populations. This shows that the variation that exists between the populations (e.g., the OTUs unique to each dataset) does not significantly alter the interactions between the OTUs shared across all three.

Several community associations with age and BMI in TwinsUK replicated in the LLDEEP and Israeli-PN datasets. These associations with age and BMI, and the heritability results within TwinsUK, broadly reflected previous observations from studies investigating taxa in isolation. For instance, in TwinsUK communities negatively associated with BMI were enriched with butyrate producers Ruminococcaceae and Coprococcus eutactus (Louis & Flint, 2009). Members of the Ruminococcaceae family have previously been associated with visceral fat mass in members of the TwinsUK cohort (Beaumont et al., 2016). Negative associations with butyrate have also been previously observed with metabolic deficits such as type 2 diabetes that are also associated with obesity (Gao et al., 2009; Qin et al., 2012). However, short chain fatty acids have also been observed at higher levels in obese mice (Turnbaugh et al., 2006), and specific taxa associations with obesity were not found in a recent meta-analysis of human data using BMI (Sze & Schloss, 2016). We observed two Bifidobacterium rich communities that were negatively associated with age. This is in line with two previous studies (Yatsunenko et al., 2012; Odamaki et al., 2016), although it should be noted that these considered a wider range of ages and the principal loss of Bifidobacterium with age was observed in infants. We also observed F.prausnitzii communities that were negatively associated and a community of Enterobacteriaceae that was positively associated with age, similar to previous observations within TwinsUK with frailty (Jackson et al., 2016b). Furthermore, the most heritable communities contained OTUs belonging to highly heritable taxa (Goodrich et al., 2014).

Co-occurrence patterns between taxa

It would be expected that community level associations reflected those of taxa based studies as each community constituted of taxonomically similar OTUs. This is in agreement with a previous study that found that co-occurrence to be higher between genetically similar taxa (Chaffron et al., 2010). This indicates either that interactions might evolve within closely related taxa; or that co-occurrence mainly detects genetically related taxa (more likely to have similar functionality) responding to environmental stimuli in the same manner. Reflecting observations that niche differentiation can be a major driver of co-occurrence patterns (Levy & Borenstein, 2013). A further possibility is that co-occurrence communities are grouping reads from source taxa that are improperly captured by the heuristic sequence clustering used to generate OTUs. Whilst this is unlikely the main driver of these communities (most had mean sequence identities more divergent than the 97% threshold), further exploration is warranted to determine the influence of OTU clustering threshold on co-occurrence patterns and communities.

We also observed several more taxonomically diverse communities with higher levels of sequence divergence between OTUs. From those highlighted in the results: The Israeli-PN community 3 contained a diverse number of taxa that have been associated with short chain fatty acid (SCFA) production including the genera Coprococcus and Blautia and the Rikenellaceae family (Duncan et al., 2002; Vital et al., 2015; Louis & Flint, 2017). Similarly, numerous taxa in Community 1 from the Israeli-PN network have been associated with SCFA production, including F.prausnitzii, Roseburia, Blautia, and potentially Oscillospira (Sokol et al., 2008; Louis & Flint, 2017; Konikoff & Gophna, 2016). Community 4 from the LLDEEP network also contained a similar combination of OTUs assigned to these same SCFA associated taxa. The LLDEEP community 19 contained 29 OTUs all assigned to the order Clostridiales but with a mixture of finer taxonomic assignments including Oscillospira and the family Christensenellaceae, which, as well as being a highly heritable taxon, has been experimentally shown to influence Oscillospira abundance in mice (Goodrich et al., 2014). OTUs within Community 5 from the TwinsUK network also contained taxa related to SCFA production including Blautia and Coprococcus, and TwinsUK community 26 contained Christensenellaceae and Oscillospira OTUs, similar to LLDEEP Community 19.

These results suggest that the communities containing a more diverse range of taxa are largely driven by a shared involvement in short chain fatty acid production in the gut. This could be due to each taxa having similar responses to host environment and substrate supply, or due to metabolic inter-dependencies between them. Further work using metabolic and metagenomic data will be required to determine which is the case, and is warranted given the beneficial health effects associated with gut microbial SCFA production (Maslowski & Mackay, 2011). The conserved communities containing Christensenellaceae and Oscillospira are also of interest given their established interaction and the heritability of the former family (Goodrich et al., 2014). Stable observation of this community across human cohorts might reflect a direct host action on the gut that selectively promotes these communities.

We are unable to determine if the co-occurrence communities observed here reflect closely related taxa responding to environmental niches in a similar manner, or direct interactions between microbes (such as metabolic dependencies) driving inter-dependencies between them. Most likely it is a combination of these factors. Nevertheless, co-occurrence communities provide a useful method to reduce the dimensionality of 16S rRNA gene sequencing data to units reflecting biological phenomena. Further analyses using the described approaches should investigate the influence of host factors such as genetics, diet and health, on these communities and their reciprocal influences on the host.