Resequencing and comparison of whole mitochondrial genome to gain insight into the evolutionary status of the Shennongjia golden snub‐nosed monkey (SNJ R. roxellana)

Abstract Shennongjia Rhinopithecus roxellana (SNJ R. roxellana) is the smallest geographical population of R. roxellana. The phylogenetic relationships among its genera and species and the biogeographic processes leading to their current distribution are largely unclear. To address these issues, we resequenced and obtained a new, complete mitochondrial genome of SNJ R. roxellana by next‐generation sequencing and standard Sanger sequencing. We analyzed the gene composition, constructed a phylogenetic tree, inferred the divergence ages based on complete mitochondrial genome sequences, and analyzed the genetic divergence of 13 functional mtDNA genes. The phylogenetic tree and divergence ages showed that R. avunculus (the Tonkin snub‐nosed monkey) was the first to diverge from the Rhinopithecus genus ca. 2.47 million years ago (Ma). Rhinopithecus bieti and Rhinopithecus strykeri formed sister groups, and the second divergence from the Rhinopithecus genus occurred ca. 1.90 Ma. R. roxellana and R. brelichi diverged from the Rhinopithecus genus third, ca. 1.57 Ma. SNJ R. roxellana was the last to diverge within R. roxellana species in 0.08 Ma, and the most recent common ancestor of R. roxellana is 0.10 Ma. The analyses on gene composition showed SNJ R. roxellana was the newest geographic population of R. roxellana. The work will help to develop a more accurate protection policy for SNJ R. roxellana and facilitate further research on selection and adaptation of R. roxellana.

forests in Nujiang Prefecture, China (Geissmann et al., 2011;Liedigk et al., 2012). Rhinopithecus roxellana was distributed widely in China during the Pleistocene; however, today, wild R. roxellana populations only exist in three isolated mountainous regions, including the Minshan and Qionglai mountains in Sichuan and Gansu provinces (SG), the Qinling mountains in Shaanxi province (QL), and the Shennongjia National Nature Reserve (SNJ) in Hubei province. Rhinopithecus roxellana is classified as an endangered species by the World Conservation Union and is protected at Level I by law of the People's Republic of China on the protection of wildlife (Le, 2015;Pan et al., 2009;Zhou et al., 2014). Due to the distinct appearance of R. roxellana, a beautiful golden coat and snub nose, and its rarity, R. roxellana is an endangered arboreal primate icon in China (Figure 1). R. roxellana has a population of approximately 22,000, including approximately 16,000 individuals from the SG population, approximately 5,500 individuals from the QL population, and approximately 1,000 individuals from the SNJ population Yang et al., 2012).
Nevertheless, neither of these two hypotheses can perfectly explain the evolutionary relationships among these three populations. The Shennongjia National Nature Reserve is the most eastern habitat of R. roxellana in China, and SNJ R. roxellana has the smallest R. roxellana population. To protect the SNJ R. roxellana population, it is particularly important to study the genetic evolution and the genetic diversity of SNJ R. roxellana. Analyses of the complete mitochondrial DNA control region suggested that the QL and SNJ populations originated from the SG population, but there was no gene flow between the SNJ and SG populations (Luo, Liu, Pan, Zhao, & Li, 2012). In contrast, the structure analyses of 16 microsatellite loci showed that the SNJ population consisted of two groups and that the SNJ population originated from the SG and QL populations . Moreover, some studies showed SNJ population had the lowest diversity (Chang, Luo, Liu, et al., 2012;Luo et al., 2012), other studied showed SNJ population had rich genetic diversity (He, Zhang, Peng, Li, & Li, 2010). However, until now the phylogenetic relationship and genetic diversity of SNJ R. roxellana is still unclear.
Mitochondrial DNA (mtDNA) has been widely used as a molecular marker in phylogenetic and phylogeographic studies of the Rhinopithecus genus and other primates because it is high mutation and substitution rates, low effective population size, high copy number, rare gene recombination, maternal transmission, and easy accessibility (Finstermeier et al., 2013;Kolleck et al., 2013;Liedigk et al., 2012;Yang et al., 2012;Yu et al., 2011;Zhang, 1998). Mitochondrial genes/fragments such as the D-loop, Cytochrome b (Cytb), tRNA, and NADH gene have been used for phylogenetic studies on genetic divergence and phylogenetics of the R. roxellana species (Kolleck et al., 2013;Le, 2015;Luo et al., 2012;Sterner et al., 2006;Zhang, 1998).
The complete mitochondrial genomes were used to analyze primate phylogenetic relationships and divergence date (Pozzi et al., 2014).
In this study, we used the Illumina Hi-seq2000 platform to resequence the complete mitochondrial genome of SNJ R. roxellana.
Furthermore, we analyzed the mitochondrial genomic divergence and phylogenetic relationships within the R. roxellana species and the Rhinopithecus genus to determine whether R. roxellana is a monotypic or polytypic species and assess the divergence ages and evolutionary status of the SNJ population within R. roxellana.

| Ethics statement
The research complied with protocols approved by the Forestry Ministry of Hubei Province, China, and abided by the legal requirements

| Sample information
For sequencing, we collected a frozen liver sample of a three-monthold male R. roxellana who was from Shennongjia National Nature Reserve in China and killed during a male takeover in August 2012.
His corpse was immediately stored at −20°C after death.

| DNA extraction, short-read sequencing, de novo assembly and annotation
In the study, we resequenced the mitochondrial genome of SNJ R. roxellana. Firstly, we extracted pure liver mitochondria using a tissue mitochondria isolation kit (Beyotime Biotech Inst., Beijing, China), and then, mtDNA was extracted using the blood DNA kit (E.Z.N.A.; Omega, USA). The extract was prepared for both Sanger sequencing and next-generation sequencing (NGS). One part of the mtDNA was used to construct a paired-end sequencing library with an insertion size of approximately 500 base pairs (bp) and was sequenced by NGS with Illumina Hi-seq2000. We prepared the libraries as follows: fragments. Then, the PCR products are selected by running another 2% agarose gel to recover the target fragments. Purify the gel with QIAquick Gel Extraction kit (Qiagen). The remainder mtDNA was used to perform a PCR amplification by the primers MT (15,543) and Mt (7-535) ( Table 1) and general sequencing to link the gaps between scaffolds. The PCR products were preliminarily confirmed on a 1.0% agarose gel, sequenced directly in two reactions with forward and reverse primers by Sanger sequencing and then spliced with contigs or scaffolds using the ContigExpress software package.
We used SOAP de novo to assemble the short reads from NGS data (http://soap.genomics.org.cn). In order to successfully de novo assemble the short reads, we filtered the low-quality reads (When the quality value of a base was less than or equail 7, we defined the base as low-quality base. If there was more than 10% low-quality base in a read, the read was referred as low-quality read.), and used the highquality clear data, and used -K, -R, -F, and -u as the program parameters in the SOAP to finish assembly. Each read that was identified as mtDNA was aligned to mtDB. These alignments were then merged, and each alignment column was examined to determine the majority base, yielding large-assembled contigs or scaffolds of mtDNA. To evaluate accuracy and completeness, the assembled mitochondrial genome sequence homology analysis was conducted using online BlastN in GenBank. Mitochondrial genomes were annotated in reference to the R. roxellana mtDNA sequence (JQ821835.1). If the reading frame of protein-coding genes was disrupted, the original read assembly was revised and corrected manually.

| Genetic divergence
To estimate the genetic divergence of the SNJ R. roxellana in base composition, we performed the genomic analyses including base composition, base substitution, amino acid replacement with respect to reference mtDNA sequence, and the number of diversity base and the gaps to R. brelichi mtDNA sequence using the complete mitochondrial sequence data. Hence, alignment between SNJ R. roxellana mtDNA and QL R. roxellana mtDNA (reference sequence) was conducted. Alignment between QL R. roxellana mtDNA and R. brelichi mtDNA, alignment between SG R. roxellana mtDNA and R. brelichi mtDNA, alignment between SNJ R. roxellana mtDNA and R. brelichi mtDNA were performed.
All pairs' alignment was conducted using online BlastN in GenBank, and amino acid replacement was analyzed according to gene annotation in GenBank. We further analyzed the genetic divergence between SNJ R. roxellana and reference mtDNA based on 13 functional mtDNA genes.

| Phylogenetic reconstruction of genus Rhinopithecus
To assess the phylogenetic position of SNJ R. roxellana among the Rhinopithecus genus and the relationships within the R. roxellana species, multiple alignments and phylogenetic analyses were conducted using the neighbor-joining (NJ), maximum-likelihood (ML), and Bayesian inference (BI) methods. NJ and ML phylogenetic trees were constructed using the MEGA6.0 program (Kimura, 1980 with random trees and were sampled every 1 × 1,000 generations.
The first 5 × 10,000 trees were treated as burnt in and discarded. The Bayesian consensus tree was constructed using the remaining trees.
Internodes with posterior probabilities of 95% were considered statistically significant. In all analyses, positions containing gaps and missing data were eliminated in the phylogenetic tree construction. In this study, eleven complete mtDNA sequences of seven Rhinopithecus, one Pygathrix nigripes, one Nasalis larvatus, Trachypithecus germaini, and Presbytis melalophos were retrieved from GenBank Database (Table 2); the latter four species were used as outgroups.

| Divergence age analysis
To estimate R. roxellana divergence ages, we inferred the evolutionary history based on complete mtDNA, employing a complete molecular clock approach as implemented in MEGA6 (Tamura et al., 2013). In order to ensure a relatively accurate genetic divergence time, we first analyzed the information site and genetic distance of all sequences and then used the NJ and ML methods to infer the evolutionary history (Tamura, Nei, & Kumar, 2004;Tamura et al., 2012). All positions containing gaps and missing data were eliminated, and evolutionary analyses were conducted in MEGA6. In the NJ method, the evolutionary distances are computed using the ML method (Felsenstein, 1985) and The tree is drawn to scale, with branch lengths measured in the relative number of substitutions per site. We selected the proposed split for the divergence between SG R. roxellana (DQ355300.1) and P. nigripes of 6.57 (6.69-6.45) million year ago (MYA, Ma) as the calibration points in two models (Perelman et al., 2011).

| Assembly and annotation on mtDNA of SNJ R. roxellana
We resequenced the mtDNA of SNJ   (Table 3). These results confirmed the strong transitional bias (transitions > transversions) in primate mtDNA sequences as reported in the previous study (Zhang, 1998

| Phylogenetic analyses
In this study, after all the positions containing gaps and missing data were eliminated, 16,460 positions were used to construct the phylogenetic tree. Regardless of the phylogenetic method employed, the topology structure of three phylogenetic trees was highly congruent, so we only present the NJ tree here (Figure 2

| Divergence age analysis
Divergence times for all branching points in the topology were calculated with the RelTime method, and all branching points in the topology were calculated using the NJ method and the ML method based on the Tamura-Nei model. Codon positions included 1st + 2nd + 3rd, using the branch lengths contained in the inferred tree. In the NJ methods, the optimal tree with the sum of branch length was 0.47216046. In the ML methods, the estimated log like- The most recent common ancestor of R. roxellana lived approximately 0.10 Ma (Figure 4).
F I G U R E 2 Evolutionary relationships of Rhinopithecus genus was inferred using the neighbor-joining method. The evolutionary distances were computed using the maximum composite likelihood method and are shown as the number of base substitutions per site SNJ-R.roxellanae(KM504390.1).

| Evolutionary relationship of geographic groups of R. roxellana
The phylogenetic relationships of the Rhinopithecus genus have been studied Karanth et al., 2008;Li et al., 2010;Osterholz et al., 2008;Yang et al., 2012); however, the evolutionary relationships of R. roxellana are still in dispute (Pan et al., 2009). Some researchers insisted that the QL and SNJ populations originated from the SG population Li et al., 2007;Luo et al., 2012), while others believed that the SG population was a heterozygosis of the QL and SNJ populations (Chang, Luo, Liu, et al., 2012;Liedigk et al., 2012). The Shennongjia National Nature Reserve is the most eastern habitat of R. roxellana in China, and SNJ R. roxellana has the smallest population (Pan et al., 2009). Analyses of the complete mtDNA control region suggested that the QL and SNJ populations originated from the SG population, but there was no gene flow between the SNJ and SG populations . In contrast, the structure analyses of 16 microsatellite loci showed that the SNJ population consisted of two groups and that the SNJ population originated from the SG and QL populations .
Moreover, some studies showed SNJ population had the lowest diversity (Chang, Luo, Liu, et al., 2012;Luo et al., 2012), another study showed SNJ population had rich genetic diversity . Therefore, the origin of R. roxellana is controversial in previous studies because of different genetic markers. In this study, we used three complete mitochondrial genome sequences, report (Zhou et al., 2016). The result, to some extent, could explain why the early results were controversial.

| The divergence ages of geographic groups of R. roxellana
The divergence ages of endangered species motivate us to determine their protected status and inform protection policy. The early study on the base of autosomal genome clusters showed the divergence ages between the northern species (R. roxellana and R. brelichi) and the "Himalayan" species (R. bieti and R. strykeri) was 1.60 Ma, and the divergence ages between R. roxellana and R. brelichi was 1.69 Ma (Zhou et al., 2014). But the result on the base of mitochondrial haplotype clusters showed the divergence ages of geographic of R. roxellana were 1.17 ± 0.70 and 0.53 ± 0.30 Ma (Chang, Luo, Liu, et al., 2012). In the study, although the evolutionary tendency of the Rhinopithecus genus and R. roxellana was consistent in the NJ and ML methods, the variation in divergence times was greater in ML than in NJ, and the divergence times occurred later in ML than in NJ. In the NJ method, the divergence ages between the northern species and the "Himalayan" species were 2.41 Ma, and the divergence ages between R. roxellana and R. brelichi were 2.23 Ma, and the divergence time of the R. roxellana was from ca. 0.14-1.16 Ma, and that of SNJ R. roxellana was ca. 0.14 Ma. In the ML method, the divergence ages between the northern species and the "Himalayan" species were 1.9 Ma, and the divergence ages between R. roxellana and R. brelichi were 1.57 Ma, and the divergence time of the R. roxellana was from ca. 0.08-0.
The ND 2 and ND 6 of R. roxellana presented significantly positive adaptation to high altitude and cold weather stress. In the study, 1 gene (ND 1 ) showed positive selection, three genes (ND 2 , COXI, and COXII) displayed neutral selection, and the others exhibited purifying selection in the complete mitochondrial genome sequence of SNJ R. roxellana. These results obviously indicated more purifying selections on mitochondrial proteins in SNJ R. roxellana, consistent with strong purifying selections on mitochondrial proteins in primates (Felsenstein, 1985). Purifying selections was also known as negative selections, that is, more unfavorable environmental factors acted on SNJ R roxellana. The more unfavorable environmental factors resulted in more abundant genetic divergence and less similarity.

| CONCLUSIONS
We obtained the complete mitochondrial genome sequence of SNJ