Eco-evolutionary dynamics of gut phageome in wild gibbons (Hoolock tianxing) with seasonal diet variations

Abstract

It has been extensively studied that the gut microbiome provides animals flexibility to adapt to food variability. Yet, how gut phageome responds to diet variation of wild animals remains unexplored. Here, we analyze the eco-evolutionary dynamics of gut phageome in six wild gibbons (Hoolock tianxing) by collecting individually-resolved fresh fecal samples and parallel feeding behavior data for 15 consecutive months. Application of complementary viral and microbial metagenomics recovers 39,198 virulent and temperate phage genomes from the feces. Hierarchical cluster analyses show remarkable seasonal diet variations in gibbons. From high-fruit to high-leaf feeding period, the abundances of phage populations are seasonally fluctuated, especially driven by the increased abundance of virulent phages that kill the Lachnospiraceae hosts, and a decreased abundance of temperate phages that piggyback the Bacteroidaceae hosts. Functional profiling reveals an enrichment through horizontal gene transfers of toxin-antitoxin genes on temperate phage genomes in high-leaf season, potentially conferring benefits to their prokaryotic hosts. The phage-host ecological dynamics are driven by the coevolutionary processes which select for tail fiber and DNA primase genes on virulent and temperate phage genomes, respectively. Our results highlight complex phageome-microbiome interactions as a key feature of the gibbon gut microbial ecosystem responding to the seasonal diet.

Fig. 1. Experimental design and the performance of VMs and MMs.

a Locations of the two family groups of wild and female skywalker hoolock gibbons (group NK in red and NK in blue) and the number of associated fecal samples per month from October 2017 to December 2018. China boundary was obtained from the National Catalog Service for Geographic Information. b Overview of the experimental methods and bioinformatic pipelines used for VMs and MMs. c Distribution of phage genome size colored by the two metagenomic approaches (VM: green, MM: purple). d Distribution of phage taxonomic families recovered by the two metagenomic approaches. Inset: Comparison of the phage genome completeness (left panel) and the proportion of prophages (right panel) between the two approaches. The prophage proportion for the three replicates of ten randomly amplified VM samples were averaged. e The number of phage genomes recovered by the two different metagenomic approaches across fecal samples of the six gibbons. In boxplots, boxes represent the interquartile range (IQR), and the lines inside show the median. The lower and upper whiskers correspond to the lowest and highest values within 1.5 times the IQR. Statistical significance is based on non-parametric Wilcoxon t-test (unpaired and two-sided), and the n number (i.e., the sample size used to derive statistics) are provided for each group. Source data are provided as a Source Data file.

Fig. 2. Comparison of gibbon gut phageomes recovered through VMs and MMs.

a Venn diagram shows the number of specific and shared viral populations detected by the two approaches. The relative abundance (x-axis) and prevalence (y-axis) of specific and shared phage populations detected from VMs (b) and MMs (c). P-values in brackets indicate the statistical significance of the relative abundance or the prevalence between groups based on non-parametric Wilcoxon t test (unpaired). The density plots parallel to the axes show the distribution of the prevalence and relative abundance. The VM-specific vOTUs, shared vOTUs, and MM-specific vOTUs are colored in green, gray, and purple, respectively. d The number of specific and shared VM-vOTUs in VMs, and the number of specific and shared MM-vOTUs in MMs. Correlations between the number (e) or the normalized abundance (f) of phage populations recovered from VMs and MMs colored by the six gibbons. The adjusted R² values and best-fit lines for the linear regressions are presented. g Proportion of the virulent or temperate phage populations recovered from VMs and MMs in each sample. h Venn diagram shows the number of phage populations with different lifestyles. The overlapped section indicate phage populations that included both virulent and temperate types and were thus split into Vir-vOTUs and Tem-vOTUs. i The number of virulent and temperate viral populations in each fecal sample of the six gibbon individuals. In boxplots, boxes represent the interquartile range (IQR), and the lines inside show the median. The lower and upper whiskers correspond to the lowest and highest values within 1.5 times the IQR. Statistical significance is based on non-parametric Wilcoxon t test (unpaired and two-sided), and the n number (i.e., the sample size used to derive statistics) are provided for each group. Source data are provided as a Source Data file.

Fig. 3. Phage community variation across gibbons and their feeding seasons.

a Principal coordinates analysis (PCoA) of phage community structure (using Bray-Curtis dissimilarity) for the six gibbons after separation into virulent (left panel) and temperate (right panel) fractions. The PERMANOVA consider samples grouped by gibbon families. b Venn diagram shows the number of virulent (left panel) and temperate (right panel) phage populations in the six gibbons. c Manhattan plots show the virulent and temperate populations with abundance significantly changed between the HF and HL seasons. Solid dots represent diet-responsive phage populations identified in both gibbons A2 and B2, while the hollow dots represent diet-responsive phage populations specifically found in A2 or B2. d Bar graphs indicate variations in the mean abundance of diet-response virulent (top panel) and temperate (bottom panel) phage populations colored by host families from HF to HL seasons. The variation value was calculated by subtracting the mean abundance of phage populations in HF season from that in HL season. Only populations mostly (top ten) varied in abundance are highlighted. e The x axis represents hosts of the six Lachnospiraceae-associated Vir-vOTUs significantly increased in abundance and the 30 Bacteroidaceae-associated Tem-vOTUs significantly decreased in abundance from HF to HL season. Barplots indicate the Pearson’s correlations between the abundance of these host populations and leaf proportion in gibbons A2 and B2. Dot plots indicate the Pearson’s correlations between the abundance of Lachnospiraceae populations and corresponding virulent populations, and the abundance of Bacteroidaceae populations and corresponding temperate populations. P-values for multiple testing were adjusted using the Benjamini and Hochberg false discovery rate controlling procedure. Source data are provided as a Source Data file.

Fig. 4. Phage functional diversity and response to season diet.

a Accumulation curves of virulent and temperate phage PCs. b Venn diagram shows the number of PCs from virulent phages (Vir-PCs), temperate phages (Tem-PCs), and prokaryotic genomes (Prok-PCs). c Distribution of the ratio of shared PCs with prokaryotes in virulent and temperate genomes. Dashed lines indicate the average ratio in the distribution. Statistical significance is based on non-parametric Wilcoxon t test (unpaired). Barplots show the proportion of shared PCs between pairs of virulent phages (d) or temperate phages (e) and prokaryotes. The blue and orange bars represent the proportion of shared PCs for virulent and temperate phages and their predicted hosts in each host family, respectively. The gray bars represent the proportion of shared PCs for virulent and temperate phages and the nonhosts in each family. The error bar represent the standard deviation of the proportions. The n number in each parenthesis indicate the number of phage-host pairs in each family. Host families with most (top ten) number of phage-host pairs are shown. f The toxin-antitoxin genes identified on phage genomes and their abundance variations with dietary changes. Barplots show the number of PCs annotated as toxin and antitoxin genes found on virulent and temperate genomes. Color gradient in the heatmap denotes the Pearson’s correlations between the abundance of PCs related to toxin/antitoxin genes and dietary proportions in gibbons A2 and B2. P values for multiple testing were adjusted using the Benjamini and Hochberg false discovery rate controlling procedure. Source data are provided as a Source Data file.

Fig. 5. Phage-host interactions between lifestyles and across the sampling dates.

a Distribution of virus-host correlations across the fecal samples of the six gibbons. Dashed lines indicate the average correlation coefficients in the distribution. b Correlations between the number of virulent or temperate populations and the number of prokaryotic populations across the samples of the six gibbons. The adjusted R² values and best-fit lines for the linear regressions are presented. c Mantel’s correlation coefficients used to represent the associations between Bray-Curtis dissimilarities of virulent or temperate phage community and prokaryotic community. d Comparison of the host range (left panel) and phage range (right panel) by phage lifestyles. e The number of anti-phage systems in genomes of prokaryotes assigned to virulent and temperate phages. Dots represent mean number and segments indicate the standard deviation. f The proportion of shared SNPs identified by comparing the SNPs on virulent and temperate phage genomes at each date with all subsequent dates in gibbons A2 (left panel) and B2 (right panel). The adjusted R² values and best-fit lines for the linear regressions are presented. g Detailed functional information of annotated genes with most (top ten) number of SNPs located on virulent (left panel) and temperate (right panel) phage genomes in gibbons A2 and B2. Statistical significance is based on non-parametric Wilcoxon t test (unpaired and two-sided), and the n number (i.e., the sample size used to derive statistics) are provided for each group. Source data are provided as a Source Data file.