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 sub-sequence 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).
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.
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 high-altitude 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).
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 r2 = 0.1143, P = 0.001; weighted Unifrac distance r2 = 0.2955, P= 0.001).
Fig.5 PCOA plot based on unweighted Unifrac distances(a) and weighted Unifrac distances(b) .
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).
Microbial community structure in high-altitude populations
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. In general, general Firmicutes and Bacteroidetes constituted about 90%. The subordinate phylum were Spirochaetae(5.29±7.40%), Verrucomicrobia(2.16±2.90%), Proteobacteria(1.82±1.68%), and Actinobacteria(1.43±1.46)%; sequences from other phylums were lower than 1%. Notablely, compared with other groups, BY was enriched in Spirochaetae(15.10%).
In family level, the rate of annotation ranged from 90.5% to 97.9% (on average 95%), Ruminococcaceae(32.24±8.68%), Lachnospiraceae (12.75±7.32%) and Prevotellaceae(10.46±9.74%) were major family, plused with other represented family such as Bacteroidales_S24-7_group(7.37±6.60%), Spirochaetaceae(5.25±7.42%), 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, high-altitude 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.
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 high-altitude population indicate the important in environment adaption.
Functional profiling of microbiotas in high altitude populations.
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.