Haplotype profiles
All 334 samples were successfully amplified, the ATP6 gene of Chinese rhesus macaque was 681 bp and 50 haplotypes were identified (GenBank accession numbers: MT039741-MT039790; Figure 1; additional file 1: table S1). Except for Hap1, Hap12 and Hap17, all other 47 haplotypes were unique to each population. Hap1 was widely distributed in high altitude populations and was the only haplotype detected in Xizang jiacha (XZJC), Yunnan jianchuan (YNJC) and Xizang leiwuqi (XZLWQ) populations, Hap1 also had the highest frequency of all haplotypes (accounting for 41.02% of all 334 samples). Hap12 was detected in Sichuan muli (SCML) and Sichuan panzhihua (PZH) populations (Hap12 was the only one haplotype detected in PZH). And Hap17 was detected in Sichuan seda (SCSD) and Sichuan xiaojin (SCXJ) populations. Seven haplotypes were detected in the Hainan lingshui (HNLS) population, which was the population with the highest number of haplotypes. Although Sichuan yajiang (SCYJ) had only five samples, four haplotypes were detected (Hd = 0.900±0.161, Table. 1).
We further analyzed ATP6 haplotype sequences, total 586 conserved sites and 95 variable sites were detected, including 25 singleton sites and 70 parsimony informative sites (Figure 2A). The haplotype sequences were translated into amino acid sequences (Figure 2B) and a total of 36 amino acid mutation sites were detected. The variation of ATP6 was significantly correlated with elevation. Haplotypes from high altitude populations were highly similar, especially 7 haplotypes from the Qinghai-Tibet Plateau (QTP), in addition to Hap3 detected M52L, the remaining six haplotypes detected M52F, and encoded the same amino acid sequence. In contrast, haplotypes of low altitude populations showed higher amino acid variation rate and some of the variation sites were shared interspecific: acid variant site M124S were detected in all haplotypes, M88M was detected in all haplotypes except for Hap35 and Hap36.
Genetic diversity
The haplotype diversity (Hd), nucleotide diversity (Pi), and the average number of nucleotide difference (K) were calculated to examine the genetic diversity of rhesus macaque populations (Table. 1). The overall Pi, Hd, and K values were 0.02332 ± 0.00226, 0.802 ± 0.022 and 14.982, respectively. Since only one haplotype was detected in XZJC,XZLWQ, XZGBJD and PZH, all three indices were 0. The Hd ranging from 0.100±0.088 (XZJD) to 0.900±0.161 (SCYJ). As for Pi, which was ranging from 0.00015±0.00041 (XZJD) to 0.00800±0.00163 (HN), and the K was ranging from 0.100 (XZJD) to 5.450 (HN). Tajima's D neutral test showed that except for SCML and GX populations, the test results of fourteen populations were not significant; conform to neutral mutations, reflecting these rhesus macaque populations sizes were relatively stable, and did not experience historically population contraction or expansion (XZJC, YNJC, XZLWQ and PZH did not conduct neutral test due to only one haplotype was detected). SCML (Tajima's D =-1.86311, P<0.05), indicated that SCML population had experienced population expansion, probably due to founder effect. GX (Tajima's D = 1.00927, P < 0.05), indicated population historical of contraction.
We calculated the pairwise gene flow (Nm), coefficient of differentiation (F-Statistics, FST) and genetic distance (D) among populations to reveal the differentiation populations (Additional file 2: Table S2 and additional file 3: Table S3). Results showed that within the most QTP populations, the pairwise Nm were greater than 1, except for XZJD - QHYS (0.81100), QHYS - XZJC (0.49479) and YNJC - SCBY (0.78167). As for the low altitude populations, the pairwise Nm of PZH with XZJC, YNJC and XZLWQ were 0, and the maximum value was 1.20022 (SCXJ-SCSD),indicating lower gene flow levels. The FST also got similar results, among populations from the QTP, only QHSY-XZJD, QHYS-XZGBJD, XZJC-XZGBJD, XZJC-QHYS, YNJC-XZGBJD and YNJC-QHYS were significant. For other rest populations, all the pairwise FST tests were significant, which indicated large genetic differentiation among these populations. As for D, the genetic distances among the QTP populations were 0 or 0.01, which indicating very low genetic distance. Gene flow test results showed that the Nm among PZH with XZJC, YNJC and XZLWQ were 0, however, including this three populations, the genetic distance among PZH with the QTP populations were low (0.004 - 0.007), while the genetic distance among other populations were higher, even SCML with a shared haplotype with PZH (D = 0.021).
Phylogenetic analysis
All 50 haplotypes of rhesus macaque ATP6 genes yielded a well-resolved Bayesian (BI) tree (Figure 3) by using M. thibetana (EU294187) and M. fascicularis (FJ906803) as outgroups, which clustered into two major clades (clade A and B). Clade A could be further divided into three sub-clades (sub-clade Ⅰ, Ⅱ and Ⅲ). Sub-clade Ⅰ contained haplotypes from SCSD, Sichuan heishui (SCHS), SCXJ and Sichuan hanyuan (SCHY), sub-clade Ⅱ was consisted of Guizhou xishui (GZXS) and Sichuan gulin (SCGL), and sub-clade Ⅲ contained haplotypes from PZH, SCML, SCYJ and haplotypes from the QTP. Clade B contained haplotypes from Guangxi (GX), HNLS and he’ nan (HN).
To further analyze the clustering and distribution of rhesus macaque ATP6 gene, A Median-joining (MJ) network was diagrammed (Figure 4), which had an analogical topology to the BI tree. Consistent with the BI tree, Hap12 appeared to be the common ancestor of SCML, PZH, SCYJ and Tibetan plateau populations. Hap1 appeared to be at the center of haplotypes from the QTP, namely, the rest of the haplotypes from the QTP appeared to evolve from Hap1 radiation. Interestingly, the six haplotypes from HN were significantly divided into two branches, one branch contained Hap46 and Hap50, and the other contained Hap45, Hap47, Hap48 and Hap49, which were clustered with haplotypes from HNLS.
In this study, we used two combinations for analysis of molecular variance (AMOVA; Table. 2): according to the elevations of the sampling areas, we divided 19 populations into three groups of high, middle and low altitude populations (test 1); seven groups according to subspecies differentiation based on previous research(27) and phylogenetic results of this study (test 2) to speculate the most suitable isolation type. In test 1, variation among group proportioned 43.13% (P<0.01), while in test 2, variation among group proportioned 90.88% (P<0.01).
Population demography
To further analyze the population history of rhesus macaque, we performed a mismatch distribution analysis of seven population groups based on the results of phylogenetic tree (Figure 5). On the whole, the mismatch distribution of rhesus macaque was unimodal, indicating that rhesus macaque experienced population expansion. Mismatch distribution of seven populations groups showed that except for HNLS (Figure 5F) and HN (Figure 5G), which were multimodal, the other five groups were unimodal, indicates that these two populations are stable.
Estimation of divergence time
Similar to previous studies (26, 28), our research of the divergence time among rhesus macaque populations (Figure 6) had also reached a phylogenetic tree that could be divided into eastern haplogroup (HNLS, HN and GX) and western haplogroup (populations from the QTP, Sichuan, Yunan and Guizhou). The divergence time of the eastern and western haplogroup was 1.07 Ma. The western haplogroup further differentiated into two sub-haplogroups and the divergence time among these two sub-haplogroups was 0.92Ma.
In relation to the divergence time between haplotypes and ATP6 haplotypes, we focused on the divergence time of two time nodes, Hap1 and Hap12. It can be seen that the divergence time of Hap1 was only 0.04Ma. Hap12 was on the outermost side of the whole branch, which was consistent with the results of phylogenetic analysis, and the divergence time was 0.25Ma.