Functional host-specific adaptation of the intestinal microbiome in hominids

Abstract

Fine-scale knowledge of the changes in composition and function of the human gut microbiome compared that of our closest relatives is critical for understanding the evolutionary processes underlying its developmental trajectory. To infer taxonomic and functional changes in the gut microbiome across hominids at different timescales, we perform high-resolution metagenomic-based analyzes of the fecal microbiome from over two hundred samples including diverse human populations, as well as wild-living chimpanzees, bonobos, and gorillas. We find human-associated taxa depleted within non-human apes and patterns of host-specific gut microbiota, suggesting the widespread acquisition of novel microbial clades along the evolutionary divergence of hosts. In contrast, we reveal multiple lines of evidence for a pervasive loss of diversity in human populations in correlation with a high Human Development Index, including evolutionarily conserved clades. Similarly, patterns of co-phylogeny between microbes and hosts are found to be disrupted in humans. Together with identifying individual microbial taxa and functional adaptations that correlate to host phylogeny, these findings offer insights into specific candidates playing a role in the diverging trajectories of the gut microbiome of hominids. We find that repeated horizontal gene transfer and gene loss, as well as the adaptation to transient microaerobic conditions appear to have played a role in the evolution of the human gut microbiome.

Fig. 1. Community-level specificities of human and NHA fecal microbiomes.

a Phylogenetic Diversity across host groups at a sampling depth of 1 Mio. mapped reads per sample. b Ordination of unweighted UniFrac distances of all samples, colored by host subgroups. c Tanglegram of host (left) and microbiome (right) trees, the latter based on unweighted UniFrac distances. d SGB sharing coefficients between host group. Rows represent reference host groups; columns represent the groups with which they share overlap. Numbers in the tiles are P-values from the analysis for enrichment (↑) and depletion (↓) in the reference group based on random 1000 permuations (unadjusted, one-sided). All analyzes are based on n = 211 independent samples (n_Human= 71, n_Chimpanzee = 91, n_Bonobo = 12, n_Gorilla = 40). Boxplots show the following elements: center line: median, box limits: upper and lower quartile; whiskers: 1.5 × interquartile ranges.

Fig. 2. Taxonomic differences the microbiota of humans and NHAs.

a Effect sizes (t-values from univariate linear regression) of abundance differences of genera in fecal samples of NHAs (n = 143) and humans (x-axis) and humans within African (n = 23) and European (n = 48) populations separated as well (y-axis). Points are colored according to the association groups with European (dark blue) and African human populations (light blue), or to indicate enrichment in NHAs (orange) or humans (bright blue). Taxonomic groups not found to be associated with any of the groups (all Q > 0.05, two-sided) are shown in grey. Horizontal and vertical lines depict the t-value threshold ( | t-value | > 4.04) for statistical significance after Bonferroni-correction. b Per-sample and host group abundances of selected genera found with unchanged abundances across all groups (top), increased abundance in NHAs (n = 143) or humans (n = 71; rows 2 and 3, respectively), or in humans from Africa (n = 23) or Europe (n = 48; rows 4 and 5). Points are colored according to host genus: Gorilla = greens, Pan = reds, oranges, and yellows, and human = blues. c Cumulative abundance trajectories of taxa associated with human communities and NHAs. Shown are the per-sample cumulative abundances within each host group, grouped based on a taxon’s association with either NHAs (n = 143), all humans(n = 71), or one of the human population subgroups. All boxplots show the following elements: center line: median, box limits: upper and lower quartile; whiskers: 1.5 × interquartile ranges.

Fig. 3. Functional differences in the microbiota of humans and NHAs.

a Effect sizes (t-values from univariate linear regression) of the abundance differences of KEGG orthologs in fecal samples of NHAs (n = 143) and humans (n = 47, x-axis) and humans within African (n = 23) and European (n = 48) populations separated as well (y-axis). Points are colored according to the association groups with European (dark blue) and African human populations (light blue), or according to general enrichment in NHAs (orange) or humans (mid-blue). Taxonomic groups not found associated with any of the groups (all Q > 0.05, two-sided) are shown in grey. Horizontal and vertical lines depict the t-value threshold ( | t-value | > 4.77) for statistical significance after Bonferroni-correction. b KEGG ortholog (KO) effect sizes from the previous analysis for differential abundance between African vs. European human population and NHA vs. Human associated taxa, respectively. Shown KOs are ordered into functional higher-level KEGG categories that were found enriched (Q_Fisher’s<0.05, two-sided) among KOs with significantly different abundances between groups. Horizontal bars indicate median t-values of all KOs in a KEGG category as an estimate for the direction of the enrichment. c KOs found enriched in the pangenomes of humans from Africa or Europe across four bacterial families shared across continents. Shown are the Z-values of the fixed-effects meta-analysis of KOs for enrichment across microbial families. KOs are sorted into functional higher-level KEGG categories that were found enriched (Q_Fisher<0.05, two-sided) among KOs with significantly different prevalence (Q_Meta<0.05, two-sided) between African and European pangenomes. Horizontal bars depict median Z_Meta-values of KOs within higher-level KEGG categories. d Prevalence of the urocanate hydratase gene (K01712) in clades found higher abundant in humans from Europe across four microbial families. P-values (Fisher’s exact test, two-sided): P_{Bacteroidaceae} = 4.55 × 10⁻⁷, P_{Lachnospiraceae} = 0.017, P_{Oscillospiraceae} = 0.052, P_{Ruminococcaceae} = 0.030. Stars indicate per-clade differences in gene prevalence (Fisher’s exact test, two-sided): *P < 0.05, **P < 0.01, ***P < 0.001.

Fig. 4. Cross-microbial clade functional associations with NHA and human hosts.

a KO terms consistently enriched in human- (blue diamonds) and NHA-associated (orange diamonds) taxonomic clades. Individual genus-level P_Fisher -values are shown in points colored by phylum. Shown P-values are unadjusted for taxon-level tests and Bonferroni-corrected for multiple testing for the meta-analysis. All tests were two-sided. b Prevalence patterns of KO terms in human- (blue) and NHA-associated (orange) Prevotella SGBs. Filled shapes represent KOs with significant (Q < 0.05, Fisher’s exact test, two-sided) differences in the statistical test. c Results of the tree reconciliation analysis for the cydA gene in Prevotella SGBs found in humans (blue) and NHAs (Pan: orange; Gorilla: green) demonstrate a history of frequent transfer events across 1000 reconciliations with random seeds. Filled and empty shapes represent cydA-positive and -negative SGBs, respectively. Red arrows depict gene transfer events from a donor to a recipient node found in at least 50% of reconciliations and are weighted by frequency. Red triangles mark nodes that were identified as gene transfer recipients with >50% frequency independent of the donor node. Black circles mark speciation events with >50% frequency. The ten highest abundant Prevotella species with established names are shown for orientation. The Prevotella tree was rooted using Paraprevotella clara as the outgroup (not shown).

Fig. 5. Cophylogeny across humans and non-human African great apes.

a Subtree phylogenies of groups with significant results in the Mantel-test based analysis for co-phylogeny. Tip colors and shapes correspond to the host subgroups. Trees were rooted on a randomly selected outgroup from a related family (not shown). b Enrichment and depletion of co-phylogeny patterns across microbial families with at least 10 SGBs in the analysis. Bars are colored by phylum, corresponding to the colors in Fig. 1. Filled bars denote significant (Q_Fisher<0.05, two-sided) enrichment and depletion. The dashed line represents the average co-phylogeny ratio across all SGBs. c Per-sample proportion of SGBs with co-phylogeny patterns across and colored by host subgroups. Group sizes are given in the x-axis labels. d Enrichment (blue) and depletion (red) of 43 in-silico inferred microbial traits, genome size and gene count in association with cophylogeny signals. Effect sizes and P-values from mixed effects logistic regression accounting for phylogenetic relatedness of SGBs. The horizontal line marks the threshold of significant Bonferroni-adjusted P-values (two-sided). e SGB-level gene counts across nine phyla, grouped by the presence (blue) and absence (red) of a cophylogeny signal. Within-phylum differences were assessed by two-sided Wilcoxon rank-sum test. f Prevalence of inferred D-Xylose utilization by SGBs across phyla, grouped by the presence (blue) and absence (red) of a cophylogeny signal. Within-phylum differences in prevalence were assessed using a two-sided Fisher-test. Group sizes of SGBs within phyla negative and positive for cophylogeny signal are given in the x-axis labels (n = neg/pos). g Tanglegram of Bifidobacterium maximum-likelihood phylogenies based on 120 GTDB marker genes (left) and gyrB sequence (right). Tip colors and shapes correspond to the host subgroups. Across all panels, stars indicate level of significance: * P < 0.05, ** P < 0.01, *** P < 0.001. Exact P-values can be obtained from Supplementary Data 14–16. All boxplots show the following elements: center line: median, box limits: upper and lower quartile; whiskers: 1.5 × interquartile ranges.