Characterization of gut microbial community of rhesus macaques under the high-altitude extreme environments

Background: The mammal intestinal microbita involved in various physiological processes in host and play a key role in host environment adaption. However, for non-human primate(NHP), little is known about their gut microbial community in high-altitude extreme environment and much less to their adaption to high-altitude environment. In this study, we want to characterize gut microbial community of rhesus macaques from multiple high-altitude environment and by comparing it to low-altitude control group to reveal the differences between altitudes. Results: we collected the fecal samples of rhesus macaques from four high-altitude populations (above 3000m) and one low-altitude population (below 100m). We first analyzed the overlap of operational taxonomic units (OTUs) between populations and found 27.8% of OTUs (core OTUs) were shared by all five population.The majority of these OTUs have a higher abundance, whereas the unique OTUs have a lower abundance. By calculating alpha diversity index, we found high-altitude populations exhibited higher diversity. Statistical analysis of beta diversity indicated there were significant difference between high and low altitude population. Significant difference in composition were detected at phylum and family. At phylum level, high-altitude gut microbial community were dominanted by Firmicutes(63.7%), but low-altitude were dominated by Bacteroidetes(52.2%). At family level, high-alititude population were dominanted by Ruminococcaceae(36.4%), but low-altitude were dominated by Prevotellaceae(43.9%). Additionally, the abundance of Christensenellaceae are significantly higher in all high-altitude populations(3.33%) than low-altitude population(0.77%), despite a low abundance in two altitudes. Finally, function prediction indicated there was a significant difference in gene copy number of 29 level-2 pathway between high and low altitude population; and 26 of them are higher in high-altitude, especially in membrane transport and carbohydrate metabolism. Conclusions: We found

mechanism of NHPs to natural environments. So far evidence has accumulated that host adaption are associated with gut microbiota which arouse our interest: is the gut microbiota of Rhesus Macaque extend their adaption to various environment especially in extremely environment such as the Tibetan Plateau? However, researches focused on gut microbiotas assisting in high altitude adaptation are limited and mainly conducted on high-altitude herbivores such as pika, yak and tibet sheep. These animal have a large differences with NHPs or rhesus macaques in diet and phylogenetic relationship which may have different adaption mechanism on microbita. Although our preliminary work revealed that significant differences in diversity and the abundance of the core common microbiota between rhesus macaque populations with various altitude gradient, little is known about even the basic feature of high-altitude rhesus macaque microbial community [44].
Here, we collected fecal samples from four high-altitude (above 3000m) wild rhesus macaques populations. Thus we can explore gut microbial community characteristics of high-altitude macaque from multiple populations and get a general conclusion. To more explicitly characterize the gut microbiome in high-altitude populations, samples from one population living in low-altitude environment (below 100m) also collected.
In this study we are trying to answer the question: What is the tyical feature of the gut microbiota diversity, composition and function in high altitude populations?

Ethics Statement and Fecal Sample Collection
Before sample collection, all the animal work was approved by the Institutional Animal Care and Use Committee of the Sichuan Agricultural University (permit number SKY-S20171007). All Rhesus macaques fecal were collected at natural habitats with little disturbance from human and livestock. Samples were collected with sterile gloves immediately after each wild rhesus macaque had defecated, and were kept cool in a thermos with ice pack. The samples used for bacterial DNA preparation were taken from the inside of the feces under sterile conditions.Then samples were transported to the laboratory on dry ice and stored in -80°C before DNA extraction.
There are five Rhesus macaques populations involved in this study, four groups in highaltitude habitats and one in low-altitude habitats. The low-altitude fecal samples were collected at LingShui (shorted as LS). LS with an average elevation of 100 meters below(in Hainan province) belongs to tropical monsoon climate with a indistinct seasons, highoxygen and high-temperature and producing a considerable amount of high quality food (especially for mango and coconut in the study site) throughout the year, which is a hospitality environment for Rhesus macaques. All high-altitude habitats at the Tibet Plateau with a altitude above 3000m including JiangDa, BaiYu, YuShu and GongBujiangda shorted as JD, BY, YS and GB respectively. By the way, JD population is relatively closed to BY population (crow-fly distance 70km) but seperated by the Yangtze river. This is for the propose that to investigate if the river could be a potential geographical barrier to gut microbial community of rhesus macaques.

DNA Extraction and Sequencing
DNA was extracted from fecal samples using a TIANamp Stool DNA kit (Tiangen, Beijing, China), following the manufacturer's instructions. The integrity of the extracted genomic DNA was verified by 1.0% agarose gel electrophoresis. V3-V4 regions of bacterial 16S rRNA gene (from 341 to 806) were amplifed from extracted DNA using barcoded primers 341 F (5′-CCTACGGGNGGCWGCAG -3′) and 806 R (5′-GGACTACNVGGGTATCTAAT-3′). PCR was performed in a 50 μL reaction system containing 1.5 μL of each primer, 100 ng template DNA, 5 μL10× KOD Buffer, 5 μL 2.5 mM dNTPs and 1 μL KOD polymerase. The PCR conditions consisted of a denaturation step at 95 °C for 2 min, the amplifcations were carried out with 27 cycles at a melting temperature of 98 °C for 10 sec, an annealing temperature of 62 °C for 30 sec, and an extension temperature of 68 °C for 30 sec.
Finally, an extra extension step at 68 °C for 10 min was performed. The barcoded PCR products were purified using a DNA gel extraction kit (Axygen, China) and quantified using QuantiFluorTM real-time PCR. Then,next-generation sequencing was performed by Illumina Hiseq 2500 PE250, which was conducted by Genedenovo Inc. (Guangzhou, China).

Data Analyses
Low quality sequences that containing the ambiguous base(N) more than 10% were removed, and sequences with high-quality base (quality score above 20) less than 60% were also removed. Then tags were merged using the FLASH program (version1.2.11) with default parameters. Low-quality contigs were removed using Qiime(V1.9.1)[45]. After removing chimera, the sofware MOTHUR was used to remove the redundant tags to get unique tags [46,47]. The obtained unique tags were then used to calculate the abundance.The demultiplexed reads were clustered at 97% sequence identity into operational taxonomic units (OTUs) using the UPARSE pipeline [48]. OTUs were classified using RDP classifier and the Silva version 123 and Greengenes version gg_13_5 reference set with a 80% confidence threshold [49][50][51]. After annotation, we find three samples(YS-1 YS-2 LS-9) were abnormal higher in Proteobacteria and Actinobacteria which may reflect the dysbiosis in gut microbiota [52]. So we removed these samples in the subsequent analyses.
The alpha diversity index Shannon, Simpson, Chao1, ACE index, and rarefaction curves were calculated using Qiime. The Beta Diversity metrics ,the weighted and unweighted UniFrac distance matrices were calculated and visualized with R statistical software.
Statistics of welch's t-test, wilconxon rank sum test, adonis and anosim test was also calculated using R.The function of gut bacteria communities were predicted using phylogenetic investigation of communities by reconstruction of unobserved states

Sequencing Profiles
After quality filtering, 6873729 16S ribosomal RNA gene sequences were obtained from 50 rhesus macaque feacl sample (137474±21752 per sample). Then sequences were clustered at 97% sequence identity and 3005 otus were generated(1234±181 per sample ) ( Table S1). The rarefaction curves had already reached a plateau at this sequencing depth( Fig.S1), suggesting that the sequencing depth had met the demand for subsequence analysis.

Gut microbiotas diversity in high altitude populations
Firstly, we calculated the alpha diversity indexs to evaluate the gut microbiotas diversity in high altitude populations. The Chao1 (1519±147), ACE (1494±147), sob (1158±124), shannon (6.670±0.2280) and simpson(1-D=0.9704±0.005958) index have been profiled. In detail, GB and YS showed a similarly higher Chao1 index are 1659 and 1658; JD was in the middle(1444) and BY group was relatively lower in 1314( Fig.2; specific data of each samples are presented in table S2 ). Specially, the Chao1 index showed significant difference between all high-altitude groups(two-tailed t-test P<0.05) except between GB and YS(two-tailed t-test P=0.98). The results of ACE, sob and shannon indexs were similar to Chao1 indexs. In contrast, the differences in simpson indexs between all groups were more slightly. Simpson indexs were highest in the GB group and only reached significant level with the lowest group(two-tailed t-test P<0.05).

Fig.2
When compared with low-altitude populations, the high-altitude populations exhibited a series feature of high alpha diversity, for example, shannon and simpson indexs are higher in all high-altitude populations, although lack a significant difference when compared with BY and JD separately(comparison of alpha diversity between all groups can be founded in Table S3 ). But if we merged all high-altitude populations as a group and compared with low-altitude population, the difference reached a significant level (Fig.3 two-tailed t-test for shannon p=0.01 and for simpson p=0.005). This result also was validated by Wilcoxon rank sum test even with a stricter confidence level(p=0.007 and p=0.0003 for shannon and simpson respectively). Other alpha diversity indexs such as Chao1 index (Fig.3) also indicated that the high-altitude population with a relatively high alpha diversity.

Fig.3
We measured beta diversity index (bray, jaccard, weighted and unweighted Unifrac distances) to untangle inter-individual variation among all Rhesus macaques gut microbial communities. Based on weighted Unifrac distances, We found inter-individual variation were much higher in high-altitude populations GB and JD and the other two populations YS and BY were relatively lower. And this result was maintained when using bray distances and jaccard distances (Fig.S2). Given that we detected a higher alpha diversity in highaltitude populations, we compared beta divierity between two altitudes and found it also significant higher in high-altitude population based on weighted Unifrac distances (Fig.S2 P<0.05 between all group, except between BaiYu and LingShui P=0.08).

Fig.4
Based on weighted and unweighted Unifrac distances, UPGMA cluster analysis have been conducted to visualize the results of beta diversity. We conduct low-altitude population LS as the outgroup, then the unweighted Unifrac UPGMA clustering tree divided into two clade, I and II. CladeI consisted of high-altitude population BY; other three high-altitude populations were on clade II (Fig.4a). Although the geography distance between BY, JD and YS were all similar and close(about 100km), JD had closer relationship with low-altitude population LS while samples from YS clustered tightly with GB's which had a farther distance. We also found that all samples were clustered by regional distribution clearly in unweighted Unifrac UPGMA clustering tree, but it can't distinct low-altitude population from high-altitude population.
However, patterns in the weighted Unifrac UPGMA tree were different. The most obvious difference was all samples were clustered by altitudes (Fig.4b). In the high-altitude clade, some high-altitude samples were deviated from their region and clustered into other population which indicating the pattern that samples distincted by regional distribution showed a trend of dispersion when counted the sequence abundance.
Next we performed PCoA to directly visualize the relationship of beta diversity distance among four high-altitude populations. In the weighted Unifrac PCoA polt, samples from high-altitude population GB , YS and JD exhibited the closer relationship while BY distributed below the X-axis and seperated with above-mentioned high-altitude samples slightly (Fig.5). On the other hand, high-altitude population mainly concentrated on the left of Y-axis and Seperated from low-altitude population in the X-axis. And Pcoa1 account main difference for 42.65% while Pcoa2 account for 14.14%. Suggesting the gut microbial community structure of high-altitude populations differed from low-altitude's. To further validate these differences, we conducted a analysis of similarity (ANOSIM) test on weighted Unifrac and unweighted Unifrac distance results. The result of ANOSIM test proved that there were significant difference between low-altitude and high-altitude populations(unweighted Unifrac distance r = 0.4895, P < 0.001; weighted Unifrac distance r = 0.6406, P < 0.001). We also performed a permutational multivariate analysis of variance (PERMANOVA) on weighted Unifrac and unweighted Unifrac distance results.The PERMANOVA results coincided with the ANOSIM(unweighted Unifrac distance r 2 = 0.1143, P = 0.001; weighted Unifrac distance r 2 = 0.2955, P= 0.001). In short, these data suggested that the gut microbiotas of high-altitude population exhibted a feature of relatively high alpha diversity and beta diversity. And it clearly distincted from low-altitude population in beta diversity distances when the information of abundance was included(weighted).

Fig.6
To figure out the microbita component feature in high-altitude populations, we conducted analysis on the OTU and gut microbiota composition. In this study, sequences more than 99.95% has been annotated into phylum level, the predominant Phylum of high-altitude rhesus monkey populations were Firmicutes(63.16±14.71%) and Bacteroidetes (Fig.6 a; 24.25±12.07%). Firmicutes were highest in GB group accounting total sequecnes for 75.45%; this number were little lower in YS group(72.59%); and for BY were 55.45%; and even in the lowest group JD(53.05%) this number were significant higher than low-altitude group (35.98% T-test P<0.001). There were a negative correlation between Firmicutes and Bacteroidetes, which means the proportion of Bacteroidetes in GB(16.95%) were less than other high-altitude populations such as BY(20.12%) and YS(22.82%). Notably, the Bacteroidetes of JD were much more higher than other high-altitude groups but still lower than low-altitude group. Christensenellaceae, Rikenellaceae and Veillonellaceae they contributed sequences more than 80% (Fig.6 b). Particularly, the abundance of Prevotellaceae were highest in JD while Ruminococcaceae were lowest. Specially, when compared with low-altitude group, highaltitude populations contained more Firmicutes(two tailed t-test P<10 -9 ) but less Bacteroidetes(two tailed t-test P<10 -5 ) than low-altitude population. Thus the ratio of Firmicutes to Bacteroidetes in high-altitude populations was more than three times as large as low-altitude population.
To further explore the distinguishing feature in species composition in high altitude populations, t-test was used in family level. There were some significant differences in the abundances of the gut bacterial communities between samples from different altitude.
Christensenellaceae(two tailed t-test P<0.01) and Ruminococcaceae(two tailed t-test P<0.01) were significant high in high-altitude population, and Prevotellaceae was significant high in low-altitude population( fig.7). In addtion, Ruminococcaceae and Prevotellaceae were the largest family to high-altitude and low-altitude population, respectively. It is also worth noting that these discrepancy has been detected in all four high-altitude populations when compared with low-altitude population seperately.

Fig.7
Based on OTUs data, we conducted a venn diagram to explore the common OTUs between populations. To improve the reliability of our data and reflect the gut microbiota community of each group we filtered OTUs which appeared in less than two individual at the same group, and 1642 OTUs were remained. Furthermore, 458 OTUs were shared by all five populations commonly, namely they were rhesus macaques core OTUs, when calculated by abundance they constituted total sequence for 90%( Fig.S3 and Table S4).
Besides core OTUs, we found there were 82 OTUs present in all four high-altitude populations constituting total sequences for 2%. Most of these(65 of 82) OTUs were belonged to Firmicutes at the phylum, and most microbes in this phylum belong to the Ruminococcaceae(40 of 82) at family level. This result suggested that differences in unique OTUs might be a small part in abundance, but it was highly consistent in all highaltitude population indicate the important in environment adaption.
Functional profiling of microbiotas in high altitude populations.

Fig.8
To evaluate metabolic function differences between gut microbiotas communities associated with high altitude, we used Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) to infer microbial community function.
Microbial community function has been assigned to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway including 35 level 2 and 221 level 3 ortholog groups. NSTI number used to estimate reliability of this method are at a acceptable value (Fig S4). The gene copy number were highest in the high-altitude population GB(65445) and YS(65163), and BY(57831) and JD(57145) were at a lower copy number. And which was significant higher than low-altitude group(52049) (P=0.0005). In order to verify if this high gene copy is affected by sequence abundance, we rarefied sequence into same reads(30000), and got a similar result. Thus, there were pervasive differences pathway between two altitude populations including 29 level 2 and 160 level 3 ortholog groups(t-test P<0.05). Strikingly, almost all differences pathway were higher in high-altitude population. Functions enriched in high-altitude populations were involved in a wide range of processes, for instance, Membrane Transport, Carbohydrate Metabolism, Amino Acid Metabolism, Replication and Repair, Translation, Energy Metabolism and so on (Fig.8). Specially, only 3 of 29 differences pathway in level 2 were higher in low-altitude. They were Glycan Biosynthesis and Metabolism and Digestive System and Immune System Diseases, but only Glycan Biosynthesis and Metabolism reached a worthy of notice level in abundance(3%), others were in trace abundance (0.04%-0.1%). We changed the species annotation database from SILVA to greengene, replaced t-test with wilcoxon rank sum test to further confirm this phenomenon and came up with a very similar result.
In conclusion, these data indicate that pervasive biological processes were elevated in intestinal microbes of high-altitude Rhesus Macaques when compared with low-altitude population, and mainly contributed by abundance difference in common OTUs.

Discussion
Gut microbial community of high-altitude Rhesus macaques are dominated by Firmicutes  [60]. It has also been validated that Prevotellaceae have a positive correlation with fruit consumption in NHPs like western lowland gorillas and Verreaux's sifakas [59,61]. And this result is conincident with our observation that high altitude habitat belonging to plateau climate with a limited fruit production. However, a research focused on high-altitude ruminants found Prevotella is higher than their relatives in lowaltitude [41]. But there are also difference in phylogenetic relationship, physiology and food habit between ruminants and NHPs which are all influence factors to gut microbial community [62]. This discrepancy implies various species may have diverse changes in gut microbita composition to high-altitude environment adaptation.
In addition, our previous work based on elevational gradient found abundance of Therefore, we could infer that Christensenellaceae are pervasive higher in high-altitude rhesus macaques populations and may play a critical role in high-altitude adaption.
Previous research found Christensenellaceae is significantly correlated with lean body shape (body mass index (BMI)<25), and the abundances of Christensenellaceae were highly associated with lower triglyceride [19,63,64]. But these research didn't illustrate the specific mechanism that how Christensenellaceae interacte with host Physiology and how Christensenellaceae reduce weight and total adiposity gain. In our study, high-altitude Rhesus Macaque are living a cold environment with high energy requirement which is proved by increased in abundance of Firmicutes. Therefore, Christensenellaceae are less likely to reduce weight by suppressing energy harvest. We speculate the mechanism might be that Christensenellaceae can stimulate host metabolic rate and increase energy consumption. This high metabolic rate could help Rhesus Macaque maintain their body temperature during the cold winter. This hypothesis is partial implied by other study that find the relative abundance of Christensenellaceae has a markedly positive correlation with food intake, body mass change and athletic ability[65].
Specially, we find every population contains Treponema (belong to Spirochaetae) at the proportion of 1%-3% except BY group which contains Treponema up to 15%. Usually, Treponema is known for T. pallidum, the cause of syphilis and yaws. This genus also includes proficient cellulose and xylan hydrolyzers, and it is pervasive in NHPs but rare in industrialized human populations [66,67]. In a study focused on western lowland gorillas, Treponema was the only genus associated with low fruit yield [61]. Treponema could be an adaptation to diets riched in fiber, and it is possible that Treponema help the host to extract nutrients from the fibrous foods by degradating fibre. High Treponema content in BY group may imply less desirable food and high fibre intake, but the specific factor needs more investigation on surrounding areas.
In contrast to widespread in NHPs, Treponema is found only in a few rural communities with a substantial content, like Hadza whose diet is high in fiber riched plant foods[67].
Fibre degraders can help release nutrition from fibre riched foods. On the other hand, Fibre riched food could promote the growth of these species. Western diet that consisted of high fat and low cellulose may lead to the exhaustion of fibre degraders which plays an important role in evolutionary history. Study focused on gut microbiota of wild NHPs may provide a excellent view on microbiome-associated human disease.

High-altitude diversity
In the Venn diagram, all rhesus macaques groups share 458 otus(core otus) and they constitute most of microbiome of a individual which means they may be critical important to the essential function of rhesus macaques. But we also found 82 OTUs shared by four high-altitude populations with far geographical distance. Although these high-altitude OTUs are less in number, the consistency in all high-altitude populations indicate they may play a important role in high-altitude environment adaption. We also obversed that high-altitude populations exhibit a feature of high alpha diversity in multivariate alpha index, and this situation is similar to high-altitude pika [40]. Gut microbiotas composition are shaped by many factors such as genetic background, diet, environmental differences, and social behavior [20,25,70,71]. In the same population, animals are similar in genetic background and diet share and interact with the same environment, it's expectable that gut microbiotas composition of Rhesus Macaque are more similar in the same group than different group. The result of unweighted Unifrac UPGMA cluster dendrogram is at expectation that samples are clustered as geographical populations, and also indicate river has a stronger isolation effectiveness than mountain to gut microbiota community of rhesus macaques. However high-altitude populations can't be distinct from low-altitude population outgroup in this dendrogram. Interestingly, in the weighted dendrogram high-altitude populations are separated from low-altitude population, but individual doesn't clustered as geographical populations strictly. These data suggest that high-altitude environment may shape gut microbiotas composition in a similar way which beyond differences in genetic background, geographical position, and social relationships.

Pervasive enrichment of gut microbiota function in high-altitude Rhesus Macaque
All mammals have a symbiotic relationship with their microbial community which contains 10 13 bacterial cells and more than 3 million genes which enable the community to carry out multiple function [72,73]. Unlike the host genome, the microbiome can change the composition of the microbial community or evolve rapidly in individual microbial genes, resulting in modified transcriptomic, proteomic and metabolic profiles [74,75] Glycan Biosynthesis and Metabolism pathway is significant higher in low-altitude population and are concentrated in Lipopolysaccharide biosynthesis. Lipopolysaccharide is the major component of the outer membrane of Gram-negative bacteria and also a typical endotoxin [76]. It's in coincidence with our hypothesis that excessive nutrition in lowaltitude population result in reproduction of specific bacteria which is the premonition of disorder.

Conclusions
We found high-altitude rhesus macaques possesses a distinct gut microbial community The authors report no conflict of interest.

Authors' contributions
Huailiang       The compositions of high-altitude rhesus macaques gut microbiota at phylum level and family level. Abundance out of top 10 are be classified as other in this plot.

Figure 8
The comparison of gene copy number between high-altitude and low-altitude population. The significance is indicated by *P<0.05, **P<0.01, ***P<0.001

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.