Strains, functions and dynamics in the expanded Human Microbiome Project

Abstract

The characterization of baseline microbial and functional diversity in the human microbiome has enabled studies of microbiome-related disease, diversity, biogeography, and molecular function. The National Institutes of Health Human Microbiome Project has provided one of the broadest such characterizations so far. Here we introduce a second wave of data from the study, comprising 1,631 new metagenomes (2,355 total) targeting diverse body sites with multiple time points in 265 individuals. We applied updated profiling and assembly methods to provide new characterizations of microbiome personalization. Strain identification revealed subspecies clades specific to body sites; it also quantified species with phylogenetic diversity under-represented in isolate genomes. Body-wide functional profiling classified pathways into universal, human-enriched, and body site-enriched subsets. Finally, temporal analysis decomposed microbial variation into rapidly variable, moderately variable, and stable subsets. This study furthers our knowledge of baseline human microbial diversity and enables an understanding of personalized microbiome function and dynamics.

Figure 1. Personalization, niche association, and reference genome coverage in strain-level metagenomic profiles.

a, Mean phylogenetic divergences between strains of species with sufficient coverage at each targeted body site (minimum 2 strain pairs). b, Individuals tended to retain personalized strains, as visualized by a principal coordinates analysis (PCoA) plot for Actinomyces sp. oral taxon 448, in which lines connect samples from the same individual. c, Quantification of niche association (Methods; only species with sufficient coverage in at least five samples at two or more body sites). Higher values indicate greater phylogenetic separation between body sites. d, PCoA showing niche association of Haemophilus parainfluenzae, showing subspecies specialization to three different body sites. e, PCoA for Eubacterium siraeum. f, Coverage of human-associated strains by the current 16,903 reference genome set (Methods). Top 25 species by mean relative abundance when present (>0.1% relative abundance) are shown (minimum prevalence of 50 samples). Sample counts in Supplementary Table 2, and distance matrices are available from Extended Data Table 1b. PowerPoint slide

Figure 2. Core and distinguishing functions of human body site microbiomes.

a, In total, 28 metabolic pathways were core at all 6 major body sites (‘supercore’ pathways). An, anterior nares; Bm, buccal mucosa; Pf, posterior fornix; S, stool; Sp, supragingival plaque; Td, tongue dorsum. Two supercore pathways and b, 17 additional pathways were core in multiple body areas and enriched among human-associated taxa (‘human microbiome-enriched’ pathways). c, 21 pathways were considerably more abundant at 1 body site than at sites from all other body areas (‘body site-enriched’ pathways). Heat map values reflect the first quartile of relative abundance (heat maps are expanded in Extended Data Fig. 3). In pathway bar plots, total (community) abundance is log-scaled, and the contributions of the top seven genera are proportionally scaled within the total. ‘Other’ encompasses contributions from additional, known genera; ‘unclassified’ encompasses contributions of unknown taxonomy. PowerPoint slide

Figure 3. Temporal dynamics of individual species and microbial pathways at each targeted body site.

a, Jaccard similarity is maximal between technical replicates and decreases with time, although within-subject similarity always exceeds between-subject similarity. b, Gaussian process decomposition of the variance in species abundances (each point is one species; filtering criteria in Methods) into three biologically relevant components based on their characteristic timescales (Methods). Technical noise was estimated (Supplementary Table 5) but not visualized. Species with high inference uncertainty (s.e.m. on the ternary diagram > 0.2) are grey and the inference is biased towards the centre of the diagrams (Methods). Labelled version in Extended Data Fig. 5. c, Same as b, but for abundances of all core pathways. d, Illustrative time series showing dynamics at different locations within the ternary plots (Extended Data Fig. 4d–f for real examples). Sample counts in Supplementary Table 1. PowerPoint slide

Figure 4. Assembly and annotation of body-wide human microbiomes.

a–d, Tukey boxplots of total assembly size, maximum and average contig lengths, and gene counts for single and co-assemblies (sample sizes in Supplementary Table 1). e, Rarefaction curves of gene families (open reading frame (ORF) clusters at 90% sequence similarity) from predicted genes generated from single assemblies for the targeted body sites (points), with a power-law fit (lines). The size of the HMP1 WMS dataset at each body site is also shown (open circles). The rarefaction trajectory was robust to changes in the sequence similarity threshold (for the 188 posterior fornix samples, the number of gene families ranged between only 1,131,796 to 1,271,891 for similarities 70–95%). Colouring is as in axis labels of a. Sample counts in Supplementary Table 1. PowerPoint slide

Extended Data Figure 1. Extended body-wide metagenomic taxonomic profiles in HMP1-II.

a, The combined HMP1-II datasets include a total of 2,355 metagenomes (724 previously published and 1,631 new, including 252 technical replicates). These span the project’s six targeted body sites (anterior nares, buccal mucosa, supragingival plaque, tongue dorsum, stool, and posterior fornix) in addition to at least 20 samples each from 3 additional sites, of the 18 total sampled sites: retroauricular crease, palatine tonsils, and subgingival plaque. Metagenomes are now available for at least one body site for a total of 265 individuals. b, PCoA using Bray–Curtis distances among all microbes at the species level. c, Relative abundances of the most prevalent and abundant microbes (bacterial, viral, eukaryotic, and archaeal) among all body sites, as profiled by MetaPhlAn2. Prevalent eukaryotic microbes are shown at the genus level. d, Taxonomic profiles do not vary more between sequencing centres, batches, or clinical centres than they do among individuals within body sites. Ordinations show Bray–Curtis principal coordinates of species-level abundances at each body site. Within-site ecological structure is as expected, with no divergence associated with technical variables along the first two ordination axes.

Extended Data Figure 2. Geographic, temporal, and biogeographic strain variation.

a, Mean distance (Kimura two-parameter) among strains from subjects either within or between three-digit zip codes (the finest degree of geographic information available). Data and sample sizes are in Supplementary Table 2. b, Mean strain divergences between different visits for the same subject and body site compared to the mean distance between the same visits for the same subject and body site (technical replicates) for each species. c–u, PCoA plots based on the Kimura two-parameter distance are shown for Escherichia coli (c), Actinomyces johnsonii (d), and all species (that is, those shown Fig. 1b ; e–u), sorted in descending order of their niche-association score (Methods). Distance matrices used to generate these PCoAs are publicly available (Extended Data Table 1b).

Extended Data Figure 3. Core and distinguishing functions of human body site microbiomes.

This figure extends Fig. 3, with additional details and examples. a, 28 metabolic pathways were core (>75% prevalent) at all major body sites. We refer to these as ‘supercore’ pathways. b, Pathways core to greater numbers of body sites tended to have broader taxonomic ranges, with supercore pathways among the most broadly distributed (Tukey boxplots). c, 19 pathways (including two supercore pathways, starred in a) were core in multiple body areas and specifically enriched among taxa inhabiting the human microbiome (annotated to <10% of non-human-associated genera). Human-microbiome-enriched pathways include specific MetaCyc-defined variants of more broadly defined or distributed processes, for example, peptidoglycan biosynthesis (PWY-6471). d, ‘Site-enriched’ pathways are considerably more abundant at one body site than at sites from all other body areas. Black dots indicate the site where each site-enriched pathway achieved its peak abundance. Heat map values reflect the first quartile of relative abundance in a particular body site (coordinated with the percentile cutoff for a core pathway). e–g, Additional examples of the three pathway classes enumerated in a, c, and d, respectively. In each example, total (community) abundance is log-scaled, and the contributions of the top seven genera are proportionally scaled within the community total. ‘Other’ encompasses pathway contributions from genera outside of the top seven, and ‘unclassified’ encompasses pathway contributions from unidentified members of the community.

Extended Data Figure 4. Sampling interval distribution, parameter fits for simulated samples, and examples of microbial species abundance dynamics and corresponding Gaussian process fits.

a, Distributions of differences in time between samples at each targeted body site. Technical replicates are shown as Δt = 0. b, Parameter fits for simulated samples with U = 0, B = 0, T = 0.95, N = 0.05, and varying l (see Methods). Simulated samples were drawn with the real sample distribution and count from each site, to show how limitations in sampling at certain sites alter the fidelity of the fits. c, Parameter fits for five simulated samples with each of the three pure components (coloured red, green and blue), as well as all even mixtures of pairs of them (for example, yellow points are even mixtures of U and T), and even mixtures of all three (black), for differing levels of technical noise (N) and fixed l = 0.5. Uncertain inferences are more desaturated. d–f, Three examples of taxonomic profiles fit with the Gaussian process model are shown on plots designed to allow a direct comparison between the data and the fit Gaussian process, and allow the different dynamics to be visualized despite the limit of only up to three time points per person. Each example was chosen as an exemplar of one of the three non-technical components in the model. Insets denote confidence deciles of the MCMC samples. The abundance of Fusobacterium periodonticum in the tongue dorsum shows strong time-varying behaviour (d), Bacteroides stercoris in stool shows mostly inter-individual differences (e), and Gemella haemolysans in the buccal mucosa is dominated by biological noise (f). The plots show the absolute difference in arcsine square-root transformed microbial abundance (|Δx|) between pairs of samples from the same person against the difference in time between samples (points). A Gaussian-smoothed estimate of the standard deviation of the points is also shown (blue line, bandwidth three months), along with the expected difference from the fit Gaussian process (red line). The standard deviation of differences between technical replicates (points with Δt = 0 months) is also shown as the line stub at the origin, directly visualizing the level of technical noise. Biological noise is visible here as the difference between technical noise and the variance of the remaining points extrapolated to the origin. The time-varying component is visible as a gradual increase in the variance of the differences over time (that is, gradually increasing red and blue lines). Finally, inter-individual differences are visible by comparing the limit of the variance of the data with the variance of differences between subjects (green line).

Extended Data Figure 5. Gaussian process decomposition of temporal variance for metagenomic species abundances.

The posterior mean of the decomposition of variance (Methods) is shown for each species (Supplementary Table 5), coloured by phylum. Uncertainty in the estimate was assessed by the square root of the mean squared distance on the ternary plot of MCMC samples from the posterior mean, and is codified with larger points indicating more certain estimates.

Extended Data Figure 6. Assembly annotation specificity for single and co-assemblies.

a, Tukey boxplot of the percentage of proteins in each functional specificity category. b, An example Venn diagram for one set of single-sample supragingival plaque assemblies and their combined co-assembly (co-assembly on bottom left), showing counts of shared genes (computed via strict alignment) between all combinations of assemblies; the co-assembly by itself contains 96.9% of all detected genes. c, Tukey boxplot of the number of Gene Ontology (GO) terms (generated using a GO Slim with around 1,700 terms) shared between single and co-assemblies, unique to the co-assembly, or unique to one of the single assemblies, generated from a random selection of 250 assemblies across 6 body sites. Co-assemblies capture GO terms that are not in individual assemblies.

Extended Data Figure 7. Sequencing statistics and assembly quality assessments.

a, Tukey boxplots of total raw reads per sample among body sites upon retrieval from the SRA. b, Percentages of human reads marked by BMTagger per body site. c, Percentages of non-human (bacterial) reads aligned to assemblies showing assembly effectiveness (Methods; read and contig mapping to assemblies and reference genomes). d, Comparison of the number of unique Pfam domains detected in each sample by HUMAnN2 and in the assemblies, coloured by body site. Pfam domains in HUMAnN2 were considered ‘detected’ if UniRef50 sequences annotated with the domain were present in a sample at >10 reads per kilobase (around 1× coverage). Pfam domains were detected in the assemblies if they were found on a single contig by Attributor (Methods). e, Number of Pfam domains that were detected in at least 75% of samples (core’ domains) by each method, for each targeted body site. Pfams domains are stratified by unknown function.

Extended Data Figure 8. Metagenomic features abundances significantly associated with host phenotype.

a, b, Significant associations of nontrivial effect size (FDR < 0.1 and |β| > 0.01) in a multivariate linear model (significance and coefficients in Supplementary Table 5) between taxon abundances (a) and pathway abundances (b). All detected associations are independent of all other metadata, including whether the subject was breastfed, the subject’s broad dietary characterization, temperature, introitus pH, posterior fornix pH, gender, age, ethnicity, study day processed, sequencing centre, clinical centre, number of quality bases, percentage of human reads, systolic blood pressure, diastolic blood pressure, pulse, whether the subject had given birth, HMP1/HMP1-II, and BMI (group sizes in Extended Data Table 1; see Methods). Non-significant associations here should not be considered evidence of no association.

Extended Data Figure 9. Updated associations in HMP1-II.

a, HMP cohort subjects reported whether they were breastfed as infants. Remarkably, overall phylum Firmicutes abundances were lower even during adulthood (subjects’ current ages were 18–40) in individuals who had been historically breastfed. b, c, Differences associated with infant breastfeeding persisted in other clades and body sites, for example, oral Neisseria (b), although age-linked associations differed among taxa (for example, overall oral Neisseria decrease with age) (c). d–f, Examples of associations significant in the original HMP1 metagenome set that were retained in the larger HMP1-II dataset include: d, Bacteroides vulgatus in stool is significantly more abundant in Asian people compared to those of other ethnicities. e, Lactobacillus crispatus in the posterior fornix is negatively associated with vaginal pH. f, Bacteroides are significantly more abundant in individuals who have been breastfed as infants. Boxplot whiskers are defined by Tukey’s method.