Functional convergence of Yunnan snub-nosed monkey and bamboo-eating panda gut microbiomes revealing the driving by dietary flexibility on mammal gut microbiome.

The gut microbiomes of non-human primates have received a great deal of attention due to their close relationship to humans. In recent years, these studies have mainly focused on the gut microbiome of wild primates, which will be helpful to understanding the evolution of primates and their gut microbiomes (e.g., gut microbiome plasticity and diet flexibility). However, there is still a lack of basic information on the gut microbiomes from wild populations. Here, we investigated the gut microbial composition (16S rRNA gene) and function (metagenome and metagenome-assembled genomes (MAGs)) of Yunnan snub-nosed monkey populations in Weixi County, Yunnan Province, China, that had diets either completely based on wild-foraging or were regularly supplemented with human provisioned food. We found a significant difference in the gut microbiome between these two populations: the gut microbiome of the wild-foraging (no food provision) population was enriched genes involved in the detoxification of bamboo cyanide (high proportion of bamboo shoot intake) and chitin (from insect diet) digestion, while the gut microbiome of the food provisioned (e.g., fruits) wild populations were enriched genes involved in carbohydrate metabolism. Moreover, the gut microbiome of the wild-foraging population shared a putatively functional convergence with the gut microbiome of wild bamboo-eating pandas: such as microbes and genes involved in the cyanide detoxification. Therefore, the gut microbiome of the Yunnan snub-nosed monkey displayed the potential plasticity in response to diet flexibility. Long-term food-provisioning of the wild population has led to dramatic changes in gut microbial composition, function, and even antibiotic resistance. The antibiotic resistance profile for the wild Yunnan snub-nosed monkey population could be considered the baseline and an important piece of information for conservation.

Introduction

The gut microbiome impacts host nutritional intake, health, development, and behavior [1], [2], [3], [4], [5]. Host diet, phylogeny, and physiology all play an important role in shaping the mammal gut microbial composition and function [1], [6], [7], [8], [9], [10]. Non-human primates (NHP), close phylogenetic relatives of humans, provide a good example for the study of host-microbiome interactions, which would aid in understanding the co-evolution of humans and their gut microbial communities [11], [12], [13], [14], [15], [16], [17]. In current NHP gut microbiome research, there are two basic approaches. One is inter-species, the investigation into the evolution of different primate species and their respective gut microbiomes, such as the changes in the gut microbiome since the diversification of humans and ape species [18], the impacts on the primate gut microbiome by host physiology across 18 NHP species [19], the convergence of human and Old World monkey gut microbiomes [10] across 18 NHP species, and the reconstructed MAGs (over 1000 novel species) from 203 metagenomic samples spanning 22 NHP species [17].

The other is intra-species, the exploration of factors shaping NHP gut microbiomes, such as gut microbiome plasticity (seasonal fluctuations) in response to dietary flexibility (e.g., wild black howler monkeys [12], wild great apes [20], wild white-faced capuchins [21], wild geladas [22]), environment factors (e.g., habitat fragmentation) leading the dissimilarity in the gut microbiome (e.g., wild black howler monkey [23], wild Udzungwa red colobus monkeys [24]), social behavior leading to gut microbiome transmission (e.g., wild baboons [25], wild chimpanzees [26], wild red-bellied lemurs [27], captive common marmosets [28]). As these different types of studies have progressed, an increasing number of researchers have focused on the gut microbiome plasticity of wild NHPs. However, the gut microbiome of many wild NHP still requires exploration, given the number of species, over 500, and dietary diversity in wild habitats. However, considering the high activity of many wild NHPs, one of the greatest difficulties is obtaining fresh fecal samples at the level of individual members of a species, which would allow fecal samples to be linked to individual’s health and dietary information. Captivity has a profound effect on NHP gut microbiome composition, leading to an increased abundance of Prevotella and Bacteroides, which are highly associated with a decrease in dietary fibers and plant contents [29]. The uncharacterized question here is what changes may occur in gut microbiome composition and function between wild populations from the same region that either forage for all of their nutritional requirements and those provided supplemental food during the same time period?

The Yunnan snub-nosed monkey (Rhinopithecus bieti), is a rare and endangered animal endemic to China. This species mainly lives in high-altitude forests, at elevations of 3000–4400 m, in southwestern China and southeastern Tibet [30] and is considered to be the highest altitude-dwelling non-human primate [31]. Habitat fragmentation and human disturbance have seriously impacted the population [32], [33], resulting in a sharp reduction in their population [33]. In recent years, with the establishment of protected nature reserves, there has been a trend of population recovery among wild Yunnan snub-nosed monkeys [33], [34]. At present, the total population of the species has been estimated to be approximately 3,000 individuals belonging to 17 natural groups in Yunnan and Tibet [34]. R. bieti in the wild mainly feed on a variety of plants, including lichens, bamboo shoots, mature leaves, fruits seeds, young leaves, buds, flowers, bark/petiole/stem, as well as insects and fungi [35]. To make up for the frequent food shortages in the wild, supplemental foods have been provided on a regular basis to these animals in many of the nature reserves. At the provisioning sites, artificial foods consisting of native lichen and high carbohydrate foods (carrots, apples, or pumpkin seeds) are provided each day [36]. After long-term habituation to this artificial diet, it was observed that some of the provisioned wild monkey groups have come to rely on the artificial foods [36], such that their main diet (about 70%) is now lichen and the provisioned foods, and relying less upon bamboo shoots and insects. Therefore, the Yunnan snub-nosed monkey is an excellent model to investigate the potential gut microbiome plasticity in the wild and food-provisioned populations at the same time within the same region.

The majority of bamboo shoots have a high proportion of cyanide compounds [37]. The bamboo-eating giant pandas were historically distributed in the Yunnan region [38], and their gut microbiomes display several adaptive features in response to their diet. These features include a high proportion of Pseudomonas associated with the detoxification of cyanide compounds found in bamboos [39]. Here, another question was whether the wild Yunnan snub-nosed monkeys shared a potential convergence in the function of the gut microbiome with wild giant pandas. Currently, among snub-nosed monkey species, the gut microbial composition of the Guizhou snub‐nosed monkey (R. brelichi) and Sichuan snub-nosed monkey (R. roxellana) have been identified by using the 16S rRNA gene [40], [41]. Captive R. brelichi exhibit decreased gut microbial diversity and a reduced number of microbes within families that assist in the digestion of complex plant materials [40]. Therefore, the function of the Yunnan snub-nosed monkey (R. bieti) gut microbiomes is still unclear.

In this study, we investigated the gut microbial composition (16S rRNA genes) and function (metagenome and MAGs) in wild-foraging and diet-provisioned Yunnan snub-nosed monkey populations in Weixi County, Yunnan Province, China (Fig. 1). Our group has worked at this site for approximately 11 years, and we were able to collect fresh fecal samples that could be connected to individual animals. Thus, we tested two hypotheses: [1] The significant difference in the gut microbial composition and function between these two populations; [2] Did the wild-foraging population share a putatively functional convergence with the gut microbiome of wild bamboo-eating pandas due to similar diets (e.g., high bamboo diet intake)?

Fig. 1

The study area and the diet of Yunnan snub-nosed monkeys. A, Map of research sites; The red box represents the distribution range of the food-provisioned wild group (Xiangguqing group, FP); The blue box represents the distribution range of wild group (Anyi group, W). B, C, D, and E are photos of wild Rhinopithecus bieti feeding on lichen, bamboo shoots, bamboo leaves, and Sorbus thibetica. Photo F is a R. bieti feeding on insects in a dead tree, and photo G is of the insects being fed upon. Photos H and I are of the semi-provisioned wild group R. bieti feeding on provisioned food (carrot and apple). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Results and discussion

The significant differences in the gut microbial composition between the wild foraging and food-provisioned population

We found significant dissimilarity in the gut microbial community between the wild foraging (W) and food-provisioned (FP) Yunnan snub-nosed monkey (YSM) populations using 16S rRNA gene sequences. The relative abundance of Proteobacteria, Bacteroidetes, and Spirochaetes was significantly higher in the W population (20 individuals) than in the FP population (28 individuals) (Fig. 2A and 2B (Lefse analysis)). The relative abundance of Firmicutes, Actinobacteria, and Fibrobacteres were significantly higher in the FP population than in the W population (Fig. 2A and 2B (Lefse analysis), Fig. S1). At the genus level, the abundance of 14 described genera (Proteobacteria: Pseudomonas, Comamonas, Acinetobacter, Yersinia, Serratia, Massilia, and Duganella; Bacteroidetes: Flavobacterium, Chryseobacterium, Sphingobacterium, and Dysgonomonas; Firmicutes: Paenibacillus and Anaerocolumna; Spirochaetes: Treponema) were significantly enriched in the W population (Fig. 2B and Fig. S1), and 13 genera were significantly enriched in the FP population, including eight genera from Firmicutes (Roseburia, Faecalibacterium, Eubacterium coprostanoligenes group, Ruminococcus 1, Ruminococcaceae UCG-002, Ruminococcaceae UCG-005, Ruminococcaceae UCG-010, and Ruminococcaceae UCG-013), three genera from Bacteroidetes (Pedobacter, Prevotella, and Bacteroides), one genus from Proteobacteria (Parasutterella), and one genus from Fibrobacteres (Fibrobacter) (Fig. 2B and Fig. S1). The second major difference was found in the alpha diversity. The Shannon index of the food-provisioned population was significantly higher than that of the wild foraging population. Normally captive lifestyles or built environments will lead to a decrease in gut microbial diversity in mammals [42]. However, animals that rely on a single source of food (or few sources) usually present a low diverse microbiome [43]. Here, we found that providing a foreign, more diverse food supply may have led to a more highly diverse gut microbiome (Table 1). The PCoA ordination and PERMANOVA, using unweighted Unifrac distance, further supported the significant dissimilarity in the gut microbiome community between the W and FP populations (Fig. 2C; PERMANOVA test: p < 0.01). Also, the microbial community was more similar in the PCoA plot in the W animals compared to FP animals. In this study, there is a profound difference in the gut microbiome composition between the wild foraging and food-provisioned populations within the same natural region.

Fig. 2

The gut microbiome community using 16S rRNA genes in wild (W) and food-provisioned (FP) populations. A, The relative abundance of the dominant phyla in each fecal sample. B, LEfSe (Linear discriminant analysis Effect Size) was used to determine the significant difference in the abundance of gut microbiomes between W and FP populations. C, The Shannon index in each fecal sample (Wilcoxon test, p < 0.001). D, Principal Coordinates Analysis (PCoA) plot was built using unweighted UniFrac distances to assess beta diversity.

Table 1

Food contribution rate of provisioned and wild groups.

Groups	Type of food	Food source	Contribution rate %
Wild (food-provisioned) group(FP)	Lichen	wild	25.06 ± 2.19
	Bamboo shoots	wild	0
	Mature leaf	wild	14.48 ± 2.99
	Fruit/seed	wild	11 ± 2.81
	Flower	wild	2.3 ± 1.62
	Young leaf	wild	5.27 ± 0.84
	Peanut	Provisioned	7.25 ± 0.53
	Egg	Provisioned	0.58 ± 0.17
	Seed of Toxicodendron vernicifluum	Provisioned	5.61 ± 0.29
	Pumpkin seeds	Provisioned	4.13 ± 0.55
	Apple	Provisioned	9.18 ± 0.59
	Lichen	Provisioned	15.15 ± 2.12
Wild (no food-provisioned) group(W)	Lichen	wild	49.32 ± 5.05
	Bamboo shoots	wild	32.63 ± 4.34
	Mature leaf	wild	10.1 ± 1.15
	Fruit/seed	wild	4.93 ± 0.53
	Flower	wild	2.01 ± 1.44
	Young leaf	wild	1 ± 0.4

This divergence might be associated with the differences in diet composition. The proportion of bamboo shoots intake was high in the W population but rare in the FP population (Table 1). The majority of bamboo shoots have a high proportion of cyanide compounds [37]. The gut microbiomes of bamboo-eating pandas display a high proportion of specific Proteobacterial groups (e.g., Pseudomonas), which were related to the detoxification of cyanide compounds in bamboos [39], [44]. Thus, the significantly enriched abundance of these Proteobacterial genera (e.g., Pseudomonas, Comamonas, and Acinetobacter) might be a response to bamboo shoot intake in the W population. Meanwhile, the significantly enriched abundance of specific Firmicutes (e.g., Ruminococcaceae and Prevotella) and Fibrobacteres groups (e.g., Fibrobacter) might be involved in the digestion of the high carbohydrate food in the FP (food-provisioned) population (e.g., a high proportion of fruit intake). Ruminococcaceae groups and Prevotella are associated with carbohydrate metabolism [29], [45], [46]. In addition, we observed insectivorous behavior in the W population, which was rare in the FP population (Fig. 1). This discrepancy may be due to the seasonal food shortages experienced by the W population. Insects have a high proportion of chitin. Thus, we speculated the genes coding for the enzymes involved in chitin degradation would be significantly enriched in the W population as compared to the FP population. Therefore, we decided to explore the metagenomes (functional analysis) of the two populations to reveal the gut microbial plasticity (e.g., the difference in the cyanide, carbohydrate, and chitin degradation between W and FP populations) driven by dietary flexibility.

Functional divergence in the gut microbiome between the foraging and food-provisioned populations driven by dietary flexibility

We analyzed 24 fecal metagenomes (8 individuals from W and 16 individuals from FP populations) and found significant functional dissimilarity in the gut microbial functions between the foraging and food-provisioned populations. The taxonomic assignment of these 24 metagenomes confirmed high proportions of genera from Proteobacteria and Bacteroidetes (e.g., Proteobacteria: Pseudomonas, Comamonas, Acinetobacter, Yersinia, Serratia, Massilia, and Duganella; Bacteroidetes: Flavobacterium, Chryseobacterium, Sphingobacterium, and Dysgonomonas) in the W population, and a high proportion of Firmicutes (e.g., Roseburia, Faecalibacterium, Eubacterium, Ruminococcus, and other groups from Ruminococcaceae), Bacteroidetes (e.g., Prevotella), and Fibrobacteres genera (e.g., Fibrobacter) in the FP population (Fig. 3A). Interestingly, the PCA ordination and PERMANOVA analysis using Bray-Curtis distance (based on KEGG level 4: the gene abundance table) found significant divergence in the gut microbiome between the W and FP (food-provisioned) populations (Fig. 3B; PERMANOVA test: p < 0.01). Lefse analysis using the KEGG pathway abundance revealed that the abundance of many carbohydrate metabolism and energy metabolism pathways were enriched in the FP (food-provisioned) population’s gut microbiome (Fig. 3C). For example, six KEGG level 3 pathways from carbohydrate metabolism (e.g., fructose and mannose metabolism, galactose metabolism, starch, and sucrose metabolism, and pentose phosphate pathway) were significantly enriched in the FP population. The predicted taxonomy of these specific KEGG pathway analyses confirmed that the main gut microbial genera involved in these functions in FP population included Prevotella, Ruminococcus, Alistipes, Clostridium, Roseburia, Faecalibacterium, and Bacteroides (Fig. 4). Considering the high proportion of fruit intake (food-provisioned) and the high proportion of carbohydrates (e.g., sugars) [29], the divergence of gut microbial function demonstrated the response of the gut microbiome to the differences in diet between these populations.

Fig. 3

The functional analysis using 24 metagenomes (8 from W population and 16 from FP population). A, Heatmap displaying the top genera in the gut microbiome using metagenomes. B, PCA plot built using Bray-Curtis distance (the relative abundance table of the genes at KEGG level 4). C, LEfSe was used to determine the significant differences in the relative abundance of KEGG pathways between W and FP populations.

Fig. 4

Putatively predicted taxonomy of the KEGG pathways using 24 metagenomes (8 from W population and 16 from FP population). A, The main microbial genera involved in Carbohydrate metabolism in the metagenomes of the W and FP populations. B, The main microbial genera involved in fructose and mannose metabolism (belong to Carbohydrate metabolism category). C, The main microbial genera involved in starch and sucrose metabolism (belong to Carbohydrate metabolism category). D, The main microbial genera involved in galactose metabolism (belong to Carbohydrate metabolism category).

CAZy analysis (Carbohydrate-Active enzymes) further confirmed the above finding (Fig. 5). For example, the abundances of 22 GH families (Glycoside hydrolases) were significantly higher in FP population than that in W population, and the abundances of 22 GH families were enriched in the W population (Whitney U test, Bonferroni-corrected, Fig. 5A). For example, GH2 and GH43 were the two dominant GH families, and they were mainly involved in sugar metabolism (e.g., galactose, mannose, and arabinose) [47], [48], [49]. The predicted taxonomy analysis found that these two GH families putatively came from Prevotella, Alistipses, Bacteroides, Ruminococcus, Clostridium, Roseburia, and Faecalibacterium in the FP population (Fig. 5B). However, GH23 (including lysozyme type G (EC 3.2.1.17) and chitinase (EC 3.2.1.14)) was the top GH family in the W population and significantly higher than that in the FP population (Fig. 5A). The predicted taxonomy analysis found that GH23 putatively came from Cornamonas, Pseudomonas, Acinetobacter, and Rahnella in the W population (Fig. 5C). Next, we focused on chitinase (EC 3.2.1.14) in the KEGG level 4. Again, we found the abundance of genes coding for putative chitinase significantly higher in the W population than in the FP population (Fig. 6A). The predicted taxonomy analysis found that GH23 putatively came from Bacteroides, Janthinobacterium, Cellulomonas, Flavobacterium, Rahnella, Sphingobacterium, Dysgonomonas, and Pseudomonas in the W population (Fig. 6A). This finding might be related to our field observations of insectivory being common in the W population but rare in the FP population.

Fig. 5

The CAZy analysis (Carbohydrate-Active enzymes) of these 24 metagenomes (8 from W population and 16 from FP population). A, Heatmap of the relative abundance of GH families (Glycoside Hydrolases). A blue asterisk indicates significantly higher (Whitney U test with Bonferroni corrected) relative abundance of this GH family in the metagenomes of the FP population. A brown asterisk indicates significantly higher (Whitney U test with Bonferroni corrected) relative abundance of this GH family in the metagenomes of the W population. *, p < 0.05; **, p < 0.01; ***, p < 0.001. B, Circos was used to visualized the contribution of bacterial taxa (at the genus level) regarding the GH families based on the TPM of bacterial genera for the annotated GH families and the TPM of GH families in all GH families in the W population. C, Circos was used to visualized the contribution of bacteria taxon (at the genus level) regarding the GH families based on the TPM of bacteria genus for the annotated GH families and the TPM of GH families in all GH families in the FP population. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6

The genes coding for putative enzymes involved in chitin and cyanide compound degradation. A, The relative abundance of gene coding for a putative chitinase (EC 3.2.1.14) involved in chitin degradation and its predicted taxonomy (top genera). B, The relative abundance of genes coding for putative enzymes (including thiosulfate/3-mercaptopyruvate sulfurtransferase, nitrilase (TST), thiosulfate sulfurtransferase (glpE), cobalamin adenosyltransferase (EC 2.5.1.17), and nitrilase (EC 3.5.5.1)) involved in cyanide compound degradation and their predicted taxonomy (top genera). Welch's t-test (with Bonferroni-corrected) was used to test the significant difference in these genes between W (8: eight metagenomes) and FP (16: sixteen metagenomes) populations. “Others” included the low and no-rank microbial groups.

Furthermore, the abundance of some pathways from amino acid metabolism, lipid metabolism, metabolism of cofactors and vitamins, and xenobiotics biodegradation and metabolism were significantly enriched in the W population (Fig. 3C). These functional features in the gut microbiome of the W population might be associated with the high proportion of bamboo shoot intake (Table 1). Bamboo shoots are rich in proteins, minerals, and secondary compounds (e.g., cyanogenic glycosides), and are low in sugars [50], [51], [52]. Considering the previous finding in the bamboo-eating pandas (high proportion of Proteobacteria groups (e.g., Pseudomonas) related to cyanide compound degradation) [39], we then compared the genes coding for putative enzymes mainly involved in the degradation of bamboo. We found that genes coding for enzymes putatively involved in the degradation of this material (including thiosulfate/3-mercaptopyruvate sulfurtransferase, nitrilase (TST), thiosulfate sulfurtransferase (glpE), cobalamin adenosyltransferase (EC 2.5.1.17), and nitrilase (EC 3.5.5.1), were significantly enriched (Welch's t-test, Bonferroni-corrected) in the gut microbiome of the W population as compared to the FP population (Fig. 6B). The predicted taxonomy analysis found the genes coding for these enzymes mainly originated from Proteobacteria genera, such as Comamonas, Pseudomonas, Acinetobacter, Rahnella, and Variovorax in the W population (Fig. 6B). Thus, this may also reflect the gut microbial plasticity in the Yunnan snub-nosed monkey.

Dietary flexibility has been shown to impact NHP and human gut microbiome composition and function, such as wild black howler monkeys [12], wild great apes [20], wild white-faced capuchins [21], wild geladas [22], and Hadza hunter-gatherers [53]. For example, during the fruit foraging season, the gut microbiome of wild white-faced capuchins is dominated by carbohydrate metabolism, while the proportion of chitin carbohydrate-binding modules is high when eating insects [21]. The proportion of Prevotella in wild gorillas and chimpanzees is high during the fruit feeding period and may be associated the carbohydrate degradation [20]. Our study provides another example of the relationship between gut microbiome plasticity (composition and function) and dietary flexibility in a wild NHP population occupying the same region at the same time. From our results, it appears that diet might profoundly impact the gut microbiome of wild Yunan snub-nosed monkeys (YSM). However, we didn’t explore the potential effect on the host health after the profound changes in the gut microbiomes of the W and FP YSM populations. We did perform a prediction on the antibiotic resistance genes (ARGs) profiles using these metagenomes. We found a significant dissimilarity between the two populations (Fig. 7A and 7B; PERMANOVA test using Bray-Curtis distance: p < 0.05). The gut microbiome of the W population was enriched in multidrug resistance subtypes (e.g., mul_acrB, mul_mdfA, mul_mdtB, mul_mdtC, mul_mexF, mul_mexT, mul_ompR, and mul_oprM) that may originate from Pseudomonas, Acinetobacter, Rahnella, Comamonas, Variovorax, and Janthinobacterium (Fig. 7C). The gut microbiome in the FP population was enriched in subtypes of tetracycline (e.g., tet_tetQ, tet_ykkD) and vancomycin resistance (e.g., van_vanU, van_vanY, van_vanG, van_vanR, van_vanS), which may originate from Bacteroides, Prevotella, and Ruminococcus (Fig. 7D). We speculate that long-term food-provisioning to the wild YSM population had led to dramatic changes in gut microbial composition, function, and even antibiotic resistance.

Fig. 7

The antibiotic resistance genes (ARGs) profiles in the 24 metagenomes from W and FP populations. A, The relative abundance of the top types of ARGs in W and FP populations. “Others” included low abundance subtypes. B, PCoA ordination of these 24 metagenomes using Bray-Curtis distance (based on the relative abundance of the ARGs subtype). C, Circos showed the contribution of top microbial genera regarding the ARG subtypes based on the TPM of microbial genera for the annotated ARGs and the TPM of ARG subtypes in all annotated ARGs in the W population. D, Circos showed the contribution of top microbial genera regarding the ARG subtypes based on the TPM of microbial genera for the annotated ARGs and the TPM of ARG subtypes in all annotated ARGs in the FP population.

The analysis of metagenome-assembled genomes (MAGs) further confirmed the enriched genes were highly associated with degradation of the specific diet

We then investigated the gene composition in 88 high-quality MAGs derived from these 24 metagenomes using the metagenomic assembling method. The advantage of strain level (MAGs) analysis was that we could determine whether genes coding for putative enzymes involved in a specific metabolism were present in the strains. Comparative genomic analysis showed that the relative abundance of these stains was different between W and FP populations (Fig. 8), with the relative abundance of most strains from Proteobacteria, Actinobacteria, and Bacteroidetes being high in the W population but rare in the FP population, except for Helicobacter sp. MIT 14-3879, Proteobacteria bacterium CAG:495, Prevotella_multisaccharivorax, and Bacteroidetes bacterium. The relative abundances of most Firmicutes stains were high in the FP population and rare in the W population. We further confirmed that the genes coding for these putative enzymes (including thiosulfate/3-mercaptopyruvate sulfurtransferase, nitrilase (TST), thiosulfate sulfurtransferase (glpE), cobalamin adenosyltransferase (EC 2.5.1.17), and nitrilase (EC 3.5.5.1)) were mainly distributed in Proteobacteria and Bacteroidetes stains. For example, the genes coding for putative nitrilase (EC 3.5.5.1) was mainly found in Proteobacterial strains, such as Rhizobiales bacterium, Variovorax sp., Bin262: g_Delftia, Acidovorax sp. MR-S7, Oxalobacteraceae bacterium CAVE-383, Massilia_timonae, and Janthinobacterium sp. S3-2. Moreover, the gene coding for putative chitinase (EC 3.2.1.14) was mainly distributed in the W population among strains such as Dyadobacter koreensis, Sphingobacterium sp. JUb20, Pedobacter antarcticus, Janthinobacterium sp. S3-2, Cellulomonas timonensis, Sanguibacter suarezii, Carnobacterium maltaromaticum, and Vagococcus salmoninarum. Here, we found the difference in the relative abundance of the strain level (MAGs) between W and FP populations, which might be associated with the different diets. However, MAGs analyses such as this one should be taken with some caution due to the sequencing depth and the capability to assemble all MAGs in the metagenomes [54], [55].

Fig. 8

Phylogenetic analysis of the putative strains (metagenome-assembled genomes, MAGs) by binning from the 24 metagenomes from W and FP populations. The left panel displays the maximum likelihood tree using the MAGs. The number below each branch the confidence value. The middle heatmap shows the relative abundance1 of each bin in each metagenome. Relative abundance1, (the number of reads (in one metagenome) mapping to the bin)/(total number of reads in this metagenome). The histogram shows the relative abundance2 of each strain based on its mean relative abundance1 in W and FP populations. The right heatmap showed the number of genes coding for putative enzymes involved in the specific KEGG pathway in this study. These putative enzymes include chitinase (EC 3.2.1.14), thiosulfate/3-mercaptopyruvate sulfurtransferase, nitrilase (TST), thiosulfate sulfurtransferase (glpE), cobalamin adenosyltransferase (EC 2.5.1.17), and nitrilase (EC 3.5.5.1).

Functional convergence of Yunnan snub-nosed monkey and bamboo-eating gut microbiomes

One of the main characteristics in the gut microbial functions of wild foraging Yunnan snub-nosed monkeys was the enrichment in genes coding for putative enzymes involved in the degradation of the cyanide compounds present in bamboo. A similar result had been found in our previous studies on bamboo-eating pandas [39], [44], [56], [57], which were historically distributed in the Yunnan region [38]. Thus, based on our previously published 57 metagenomes (19 CA (meat-eating carnivorans), 12 HE (herbivore) [58], 10 OC (omnivorous carnivorans) [58], and 10 GP (giant pandas) [39], and 6 RP (red pandas) [39]), we tested our second hypothesis, that the wild foraging (high bamboo shoot intake and not provided with supplemental food) population shared a putative functional convergence in the gut microbiome with the wild bamboo-eating pandas. The PCoA ordination and hierarchical clustering using Bray-Curtis distance (based on the relative abundance of KEGG level 4 genes) showed the high similarity in gut microbial functions (Fig. 9A and 9B). Thus, we deduced that the potential functional convergence (high similarity) of Yunnan snub-nosed monkey and bamboo-eating gut microbiomes might be partially caused by similar diets (e.g., high bamboo diet intake), although they belong to different mammalian orders (Primates vs. Carnivora). Diet played an important role in shaping the gut microbiome in the Yunnan snub-nosed monkeys and bamboo-eating pandas. Diet drives convergence in the gut microbiomes of different mammals (e.g., opposing directionality for amino acid metabolism between carnivorous and herbivorous microbiomes [29]; composition similarity in myrmecophagous mammals [60]). The diet of the wild mammal may be more complex or diverse than we have known.

Fig. 9

Putatively functional convergence of Yunnan snub-nosed monkey and bamboo-eating panda gut microbiomes using 81 metagenomes. The PCoA ordination (A) and hierarchical clustering (B) using Bray-Curtis distance (based on the relative abundance of KEGG level 4 genes). 81 metagenomes: 24 metagenomes (W, wild population; FP, wild (food-provisioned) population) from the current study, and 57 metagenomes from previously published data (19 CA (meat-eating carnivorans), 12 HE (herbivore) (58), 10 OC (omnivorous carnivorans) (58), and 10 GP (giant pandas) (59), and 6 RP (red pandas) (59)). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Conclusion

In this study, we found a significant difference in the gut microbiome between wild foraging and food-provisioned Yunnan snub-nosed monkey populations, which might be associated with the different dietary compositions. The wild foraging (no supplemental food provided) population shared a putative functional convergence in the gut microbiome with the wild bamboo-eating pandas: such as microbes and genes involved in cyanide detoxification. Long-term food-provisioned to a wild NHP population led to dramatic changes in gut microbial composition, function, and even antibiotic resistance. The antibiotic resistance profile for the wild NHP population would be basic and important information for conservation.

Materials and methods

Study areas

Fecal samples of wild R. bieti were collected at AY (Anyi, 99°09′E, 27°27′N) and the XGQ (Xiangguqing, 99°21′E, 27°39′N). The two sampling sites were located in Weixi County, Yunnan Province, China, approximately 34 km apart (Fig. 1). The annual mean daily temperature is 9.8 °C with a maximum of 27.7 °C in July and a minimum of −9.3 °C in January. Annual rainfall is 1,371 mm(61). AY group (W foraging population) is the wild R. bieti population which is free to eat a natural diet, while XGQ group (FP, food-provisioned) is the wild R. bieti population, but supplement foods were provided regularly to them. In the wild R. bieti mainly feeds on a variety of plants, including lichens, mature leaves, fruits seeds, young leaves, bamboo shoots, buds, flowers, bark/petiole/stem, insects, and fungi [35]. The supplement foods at the provisioning site consisted the lichen and other foods (carrots, apples, or pumpkin seeds) each day [36].

Diet collection and analysis

The monkey groups were followed as close as possible without interfering with their normal activities. If conditions did not permit this, observations were continued with a binocular. The feeding behavior of the monkeys was observed by scanning instantaneous sampling at 15 min intervals to ensure the independence of the samples [62]. Between 3 May 2020 to 18 June 2020, a total of 11,807 instances of feeding individuals were recorded in the wild foraging group (Anyi group), and 19,697 instances of feeding individuals were recorded in the provisioned wild group (Xiangguqing group) from 4 May 2020 to 27 June 2020. Scanning records include the species, parts, and types of plants that the monkeys feed on. Individuals picking or placing food into their mouths were identified as feeding behavior. Feeding individuals were observed for at least 5 s to confirm the type/part of the food. Food types were divided into (1) Lichen; (2) Buds (referring to leaf buds); (3) Young leaf; (4) Mature leaf; (5) Flower (6) Fruit or seed; (7) Bamboo shoots; (8) Fungi; (9) Bark, petioles, stems; (10) Others, including insects, vertebrates, bird eggs, soil, etc. The following rules are used to determine the species of food: (1) Observe which part of the plant R. bieti was feeding on at close range, try to take pictures, and pick the leaves, flowers, and fruits of the plant for identification of species; (2) After the monkey group has left, enter the feeding site and gather the traces left to confirm the species of feeding further; (3) The act of peeling the bark off of deadwood and placing it into the mouth was considered to be feeding upon insects.

When calculating the contribution rate of different food species to the diet of R. bieti, each scanned individual was regarded as an independent sample, and the contribution rate of different food types and species was calculated using the following formula:

R = ΣC /ΣFwhere R is the contribution rate; k is the types or species of the food; ΣC is the sum of the k food types or k species; ΣF is the total number of feeding samples.

Fecal sample collection

A total of 48 fresh fecal samples were collected in June 2020 (Table S2), twenty in AY (W group) and twenty-eight in XGQ (FP group). R. bieti has evolved the habit of defecation after feeding in noon. The sampling time was arranged at 13:00–15:00. AY group samples were collected using the following process: (1) locating the group of wild R. bieti and follow them at a distance of about 200 m; (2) waiting for R. bieti to leave the feeding area; and (3) sampling the fresh fecal samples that are not mixed to make sure they belong to one individual. Fecal samples from XGQ group were collected near the provisioning site. When a group of R. bieti came to eat the provisioned food, they were observed at distances about 20–30 m and then wait for group members to defecate. After the group had left the provisioning site, the fresh stool samples were quickly collected. We had been traced food-provisioned (FP) group for about 11 years, and we could identify each individual in the group. All samples were put into sterile centrifuge tubes, sealed, labeled, and immediately frozen in liquid nitrogen, then transferred to a −80 °C freezer. DNA extraction was performed within one week after sample collection.

DNA Extraction, PCR amplification, and 16S-rRNA gene sequencing

Total DNA was extracted using the QIAamp DNA Stool Mini Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s protocol. The integrity of the nucleic acids was determined visually by electrophoresis on a 1.0% agarose gel containing ethidium bromide. The concentration and purity of each DNA extraction were determined using a Qubit dsDNA HS Assay Kit (Life Technologies, Carlsbad, CA, United States). The extracted total DNA was stored at −80℃. Only samples that meet the following criteria were used for 16S-rRNA Gene sequencing: (1) DNA integrity and no contamination; (2) DNA concentrations were > 10 ng/ul; (3) DNA total quality was > 100 ng. The V3-V4 hypervariable region of the bacterial 16S-rRNA gene was amplified from extracted total DNA with the universal bacterial barcoded primers 343F − 5′ - TACGGRAGGCAGCAG- 3′ and 798R − 5′ -AGGGTATCTAATCCT − 3′). The PCR reaction (30 µL total volume) contained 15 μL 2 × Gflex PCR Buffer, 0.6 μL Tks Gflex DNA Polymerase (1.25U/μl), 1 µL (5 pmol/μl) forward primer, 1 µL (5 pmol/μl) reverse primer, template DNA (50 ng), and 12.4 µL PCR grade water. PCR amplification program was set up according to the following procedure: 94 °C for 5 min, followed by 30 cycles of 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 20 s, and a final extension at 68 °C for 5 min. PCR amplification and paired-end sequencing were performed by Oebiotech Company (Shanghai, China) using the Illumina MiSeq platform.

16S sequence processing and analysis

Raw paired-end reads were then preprocessed using Trimmomatic software based on a sliding window (5 bp bases) for quality control [64]. Quality control criteria were as follows: (1) detect and remove ambiguous bases; (2) discard low-quality sequences with an average quality score below 20. After trimming, clean paired-end reads were assembled using FLASH software [65]. Assembled sequences were further denoised as follows: reads with ambiguous homopolymers longer than eight bp or below 200 bp were removed. Then, the clean reads were analyzed using QIIME v1.9.0 [66]. Chimeric sequences were identified and discarded using the vsearch algorithm [67], chloroplast sequences also were removed. Clean sequences were clustered into operational taxonomic units (OTUs) based on 97% similarity using a closed reference OTU-picking approach. To reduce sequencing error, singleton OTUs were removed. A representative read from each OTU was selected using the QIIME package [66]. All representative reads were annotated and blasted against the SILVA reference database (release 128, http://www.arb-silva.de/).

Before calculating the alpha and beta diversity metrics, samples with different sequencing depths were normalized by rarefying the OTU table to the minimum number (32,848) of sequences observed in all samples. The Shannon index and observed species were computed to assess the alpha diversity. Wilcoxon test for the significant difference between W and FP populations was conducted in SPSS software. LEfSe (Linear discriminant analysis Effect Size) was used to determine the significant difference in the abundance of gut microbiomes between W and FP populations [68]. The Principal Coordinates Analysis (PCoA) plot was built using unweighted UniFrac distances to assess beta diversity [69]. PERMANOVA test for group-level differences in the microbial composition was performed in QIIME package [66].

Metagenomic sequencing and analysis

Metagenomic sequencing for 24 fecal samples (8 from W and 16 from FP populations, random selection) was conducted using the Illumina HiSeq-PE150 platform. The raw data for each metagenome was about 10G, and we finally gained about 240G raw reads. The raw reads were filtered using Cutadapt [70]. BWA [71] was used to remove the host contaminations using the published Rhinopithecus bieti genome (NCBI accession number: GCF_001698545). MEGAHIT [72] was used to assemble the clean reads based on the default parameters (minimum contig length 500 bp). All coding regions (CDS) of metagenomic contigs were predicted by Prodigal [73] and were clustered by CD-HIT with these parameters (identity>=95% & overlap>=90%) [74]. Thus, we gained the unigenes after this step. The transcripts per million (TPM) were used to estimate Unigene abundance based on the number of aligned reads by bowtie2 [75]. Diamond [76] was used to conduct the alignment of unigenes against the NCBI micro-NR database (including bacteria, fungi, archaea, and viruses) with these parameters (e-value<=1e-5 & score>=60). After this step, we obtained taxonomic information for the unigenes. The functional annotations of unigenes were performed against the KEGG (Kyoto Encyclopedia of Genes and Genomes) [77]. The TPM of gut microbiome communities, KEGG pathways, and GH families per metagenome were transformed to relative abundance using STAMP [78].

LEfSe (Linear discriminant analysis Effect Size) was used to determine the significant difference in the abundance of KEGG pathways between W and FP populations [68]. The Welch's t-test (with Bonferroni correction) in STAMP [78] was used to compare the significant difference in the relative abundance of the genes coding for the putative enzymes involved in the specific dietary degradation pathways between W and FP populations. The Whitney U test (with Bonferroni correction) was used to test the significant difference in the abundance of GH families between W and FP populations. We applied PCA ordination in the vegan package [79] based on the Bray-Curtis distance matrices [80] using functional composition tables (relative abundance). Circos [81] was used to visualize the contribution of bacteria taxon (at the genus level) regarding the GH families based on the TPM of bacterial genera for the annotated GH families and the TPM of GH families in all GH families.

Based on our previously published 57 metagenomes (19 CA (meat-eating carnivorans), 12 HE (herbivore) [58], 10 OC (omnivorous carnivorans) (58), and 10 GP (giant pandas) [39], and 6 RP (red pandas) (59)), we tested our second hypothesis that the wild foraging (with high bamboo shoot intake and no supplemental food provisioning) population shared a putative functional convergence with the gut microbiome of wild bamboo-eating pandas. The PCoA ordination and hierarchical clustering using Bray-Curtis distance (based on the relative abundance of KEGG level 4 genes) were conducted in PAST4 [82].

Antibiotic resistance genes (ARGs) analysis

We blasted the identified genes against the ARDB database using SARG2.0 with the default parameters (e-value<=1e-7 & identity>=60) [83]. We then obtained the putative ARG assignment of these genes per metagenome. Next, we used custom Perl scripts to gain the abundance (TPM) of ARG subtypes for each metagenome. The TPM of ARGs per metagenome was transformed to relative abundance using STAMP [78]. The putative sequences of ARGs were blasted against the NR database in NCBI using diamond (e-value<=1e-5 & score>=60) [76]. From this, we determined the putatively predicted taxonomy of the ARGs within each metagenome. We used Circos to calculate the contribution of top microbial genera regarding the ARG subtypes based on the TPM of microbial genera for the annotated ARGs, and the TPM of ARG subtypes in all annotated ARGs. Bray-Curtis distance based on the relative abundance of ARGs subtypes was used to generate PCoA in PAST4 [82]. PERMANOVA test using Bray-Curtis dissimilarities was used to test the significant difference in the ARGs subtype community between W and FP populations.

Binning analysis (metagenomic assembled genomes (MAGs): Strain level)

We combined the clean reads of these 24 metagenomes for metagenomic assembly analysis. We use BWA [71] and Samtools [84] to map the clean reads to contigs. MetaBAT2 [85] was used to obtain the contigs for each bin based on the mapping result per metagenome. We used CheckM [86] for the quality control of each bin. We selected the high-quality bins (coverage > 80%, contamination rate < 10%) for the strain-level analysis. We used Salmon [87] to map the clean reads to these high-quality bins and determine the TPM of the bins in each metagenome. PhyloPhlAn [88] was used to construct the maximum likelihood tree for these bins. These high-quality bins were searched against the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database by using diamond [89], and then we obtained the gene composition for these bins MAGs. Thus, we could make a comparative genomic analysis to obtain the number of genes coding for putative enzymes involved in the degradation of specific dietary compounds in this study.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.