Ultra-deep sequencing of Hadza hunter-gatherers recovers vanishing gut microbes

Abstract

The gut microbiome modulates immune and metabolic health. Human microbiome data are biased toward industrialized populations, limiting our understanding of non-industrialized microbiomes. Here, we performed ultra-deep metagenomic sequencing on 351 fecal samples from the Hadza hunter-gatherers of Tanzania and comparative populations in Nepal and California. We recovered 91,662 genomes of bacteria, archaea, bacteriophages, and eukaryotes, 44% of which are absent from existing unified datasets. We identified 124 gut-resident species vanishing in industrialized populations and highlighted distinct aspects of the Hadza gut microbiome related to in situ replication rates, signatures of selection, and strain sharing. Industrialized gut microbes were found to be enriched in genes associated with oxidative stress, possibly a result of microbiome adaptation to inflammatory processes. This unparalleled view of the Hadza gut microbiome provides a valuable resource, expands our understanding of microbes capable of colonizing the human gut, and clarifies the extensive perturbation induced by the industrialized lifestyle.

Keywords: Treponema succinifaciens; human gut microbiome; hunter-gatherer; industrialization; metagenomic assembly; metagenomic metaanalysis; metagenomics; strain sharing.

Figure 1.. A vast resource of Hadza gut microbiome data.

(A) Overview of sample collection for shotgun metagenomic sequencing of Hadza fecal samples. (B) Summary of the computational workflow, tools used, and primary data generated from Hadza stool samples. (C) Number of samples versus the number of bases sequenced per sample for 21 human gut metagenomic data sets and the present study.

Figure 2.. The Hadza gut microbiota contains substantial multi-domain novelty.

(A) Phylogenetic tree of bacterial species-level representative genomes (SRGs) from Hadza and UHGG based on bacterial single copy gene alignment. Branch colors correspond to phyla. SRGs from species-level groups consisting of only genomes assembled from the Hadza or only UHGG are colored green and orange in the outer ring, respectively. The number of SRGs found in the Hadza, UHGG, or both is shown as a horizontal bar plot. Hadza genomes that correspond to organisms isolated in this study are marked with blue stars. (B) The percentage of bacteriophage species clusters assembled from the Hadza that are novel at the species level according to the MGV, categorized by phylum of the predicted host. Bacteriophages without a host prediction are labeled “Uncharacterized”. (C) A phylogenetic tree of eukaryotic genomes recovered from Hadza and Nepali gut metagenomes based on universal single copy genes. Public reference genomes are marked with blue text labels. The heatmap shows the prevalence of the individual eukaryotes in the Hadza, Nepali and Californian cohorts. (D) For each population, the percentage of predicted proteins from recovered genomes that are present in the UniRef100 and UHGP-95 protein databases. (E) The percentage of metagenomic reads mapping to various domains averaged across all metagenomic samples from each population. The phyla “Bacteriodota” and “Firmicutes_A” are shown separated from other bacteria. “Unmapped” depicts the percentage of reads that do not map to any genomes, and “Low confidence” depicts the percentage of reads that map to genomes with less than 50% genome breadth. (F) The Shannon diversity of bacteria, archaea, bacteriophage, and eukaryote genomes in metagenomes sequenced in this study. P-values from two-sided Mann-Whitney-Wilcoxon test with multiple hypothesis correction; *: p < 0.05, **: p < 0.01, ***: p < 0.001, ****: p < 0.0001, ns: p ≥ 0.05. N=332, 43, 14, and 12 metagenomes from Hadza, Nepal Ag., Nepal For., and California populations, respectively. “Nepal For.” includes the Chepang foragers, while “Nepal Ag.” includes Raute, Raji, and Tharu agrarians. (G) Collectors’ curves depicting the average number of genomes detected per sample in each population sequenced in this study after rarefaction to various sequencing depths. The vertical dotted lines indicate the average per-sample sequencing depth of this study (~23 Gbp) and the average depth of samples studied previously (~4 Gbp; ref.). Shaded areas around lines indicate 95% confidence intervals.

Figure 3.. Increased sequencing depth results in the detection of previously-uncharacterized and phylogenetically distinct taxa.

(A) The number of genomes detected in individual samples sequenced in this study when rarefying sequencing depth in 5 Gbp increments. Each set of points connected by a line represents an individual sample from which filtered reads were generated. Points and lines are colored by population. Dark lines are negative exponential curves fit to the rarefaction data for each population, shaded regions are the standard error. Vertical dotted lines indicate the average per-sample sequencing depth of this study (23 Gbp) and the average per-sample sequencing depth of samples used in Almeida et al. (4 Gbp). (B) Taxonomic distribution of organisms present at different ranges of relative abundance levels (horizontal stacked bar plots) and the percentage of species in each range that are absent from GTDB r95 (text percentages right of horizontal bars). Organisms detected at low relative abundance levels are more likely to be previously-uncharacterized than those that are more abundant.

Figure 4.. VANISH and BloSSUM taxa have distinct global prevalence, phylogeny, and functional capacity

(A) A heatmap depicting the presence of 524 SRGs (columns) within metagenomic samples from populations living different lifestyles (rows). Darker blue indicates SRG presence, lighter blue indicates SRG absence. SRGs with >30% prevalence among all samples in any lifestyle category were included. (B) SRGs were classified as “BloSSUM” or “VANISH” based on their prevalence across lifestyles (see methods for details). Colored bars correspond to columns in the heatmap. (C) A phylogenetic tree of VANISH and BloSSUM species-representative genomes. Outer ring is colored based on VANISH or BloSSUM status (magenta and blue, respectively). Branches are colored based on the phylum that tree tips belong to (“Firmicutes” includes Firmicutes, Firmicutes_A and Firmicutes_C). Black bands around the ring indicate genus or family membership of notable VANISH/BloSSUM taxa. (D) The prevalence of VANISH (magenta), BloSSUM (blue) and non-enriched taxa (gray) in the Hadza, transitional lifestyle populations, and industrial lifestyle populations. Dashed lines connect median prevalence across the taxa in each category surrounded by standard deviation (color shaded regions). Solid lines show the median prevalence for 6 representative taxa in each of these lifestyle groups. (E) Volcano plot showing enrichment of KEGG Orthologies (KOs) in either BloSSUM or VANISH genomes. Each point represents a different KO. KOs are colored based on functional categories. The x-axis is the average number of gene copies per VANISH genome minus the average number of gene copies per BloSSUM genome. The y-axis displays the p-value resulting from Fisher’s exact test with multiple hypothesis correction. The dashed horizontal line indicates the threshold of significance (p < 0.05). P-values are corrected with the Benjamini-Hochberg method. (F) A gene set enrichment analysis (GSEA) comparing VANISH and BloSSUM genomes. The 15 pathways with the highest normalized enrichment scores (NESs) are depicted, all of which reach statistical significance (p < 0.001 with false discovery rate correction; GSEA). For each pathway (row), vertical lines indicate the position of individual genes in that pathway in a ranked list of genes sorted by (prevalence in VANISH taxa - prevalence in BloSSUM taxa). The color of KEGG Pathway indicates whether the pathway is enriched in VANISH or BloSSUM taxa.

Figure 5.. Spirochaetota that are highly abundant in the Hadza are absent in industrial samples.

(A) A heatmap showing the relative abundance of the 10 most prevalent Spirochaetota species in the Hadza, Nepali, and American cohorts. All samples are sequenced to approximately the same sequencing depth. (B) A phylogenetic tree of all Spirochaetota species using genomes from NCBI, the UHGG and the species-representative genomes added in this study. Clades of commensal organisms in the genera Brachyspira, Spirochaeta, and Treponema are highlighted. (C) A phylogenetic tree of all Treponema succinifaciens MAGs in the UHGG in addition to new MAGs recovered in this study (annotated in outer ring). The inner ring is colored based on the country of origin of the individual contributing the MAG. (D) World map showing locations of populations from which T. succinifaciens MAGS were recovered as nodes (TZA = Tanzania, MDG = Madagascar, NEP = Nepal, FIJ = Fiji, PER = Peru, ELS = El Salvador). Arrows indicate the detection of transition events between populations as detected by stochastic character mapping. Thickness of the arrow indicates frequency of the transition event (thickest arrow is Tanzania to Fiji = 17.1%). The top 7 most frequent transition events are shown, accounting for 65.7% of all transitions.

Figure 6.. Microdiversity, replication rates, and patterns of strain sharing among Hadza gut bacteria.

(A) Pfams with high or low pN/pS values in Hadza fecal metagenomes. The x-axis displays the mean pN/pS value of all genes annotated with each Pfam within Hadza fecal metagenomes. The y-axis displays the probability that the number of times genes annotated as each Pfam were in the top 10% or bottom 10% of all genes on detected genomes was due to random chance (binomial test with multiple hypothesis correction). The 30 Pfams with the lowest p-values for low and high pN/pS were manually annotated with broad functional categories. (B) In situ replication rate measurements of all taxa detected in Hadza adult metagenomes across seasons. Error bars indicate 95% confidence intervals. (n.s. P > 0.05; **** p ≤ 0.0001; Wilcoxon rank-sum test). (C) The mean number of strains shared between Hadza adults broken down by various types of familial relationships. Exact p-values shown from Wilcoxon rank-sum test. (D) Rectangles along the circumference represent Hadza individuals and each link drawn between boxes indicates a shared strain. Links between members of the same bush camp are colored based on the bush camp; links between bush camps are colored black. The mean number of strains shared between members of the same bush camp and the p-value comparing strains sharing among members of that bush camp vs members from different bush camps are shown (Wilcoxon rank-sums test).