Comparative mitochondrial genomics and phylogenetic relationships of the Crossoptilon species (Phasianidae, Galliformes)

Phasianidae is a family of Galliformes containing 38 genera and approximately 138 species, which is grouped into two tribes based on their morphological features, the Pheasants and Partridges. Several studies have attempted to reconstruct the phylogenetic relationships of the Phasianidae, but many questions still remain unaddressed, such as the taxonomic status and phylogenetic relationships among Crossoptilon species. The mitochondrial genome (mitogenome) has been extensively used to infer avian genetic diversification with reasonable resolution. Here, we sequenced the entire mitogenomes of three Crossoptilon species (C. harmani, C. mantchuricum and C. crossoptilon) to investigate their evolutionary relationship among Crossoptilon species. The complete mitogenomes of C. harmani, C. mantchuricum and C. crossoptilon are 16682 bp, 16690 bp and 16680 bp in length, respectively, encoding a standard set of 13 protein-coding genes, 2 ribosomal RNA genes, 22 transfer RNA genes, and a putative control region. C. auritum and C. mantchuricum are more closely related genetically, whereas C. harmani is more closely related to C. crossoptilon. Crossoptilon has a closer relationship with Lophura, and the following phylogenetic relationship was reconstructed: ((Crossoptilon + Lophura) + (Phasianus + Chrysolophus)). The divergence time between the clades C. harmani-C. crossoptilon and C. mantchuricum-C. auritum is consistent with the uplift of the Tibetan Plateau during the Tertiary Pliocene. The Ka/Ks analysis showed that atp8 gene in the Crossoptilon likely experienced a strong selective pressure in adaptation to the plateau environment. C. auritum with C. mantchuricum and C. harmani with C. crossoptilon form two pairs of sister groups. The genetic distance between C. harmani and C. crossoptilon is far less than the interspecific distance and is close to the intraspecific distance of Crossoptilon, indicating that C. harmani is much more closely related to C. crossoptilon. Our mito-phylogenomic analysis supports the monophyly of Crossoptilon and its closer relationship with Lophura. The uplift of Tibetan Plateau is suggested to impact the divergence between C. harmani-C. crossoptilon clade and C. mantchuricum-C. auritum clade during the Tertiary Pliocene. Atp8 gene in the Crossoptilon species might have experienced a strong selective pressure for adaptation to the plateau environment.


Background
Crossoptilon, belonging to Phasianidae in Galliformes, is a rare but important genus endemic to China. The four previous recognised Crossoptilon species are C. harmani, C. mantchuricum, C. crossoptilon and C. auritum [1,2]. C. harmani is only distributed in a typical alpine and taiga habitats at elevations of 3500-3900 meters of southeastern Tibet. C. mantchuricum is mainly distributed in Lvliang Mountains of Shanxi Prov., Xiaowutai Mountains of Hebei Prov., Dongling Mountains of Beijing and some other local areas [1]. C. crossoptilon is only present in western Sichuan Prov., northwestern Yunnan Prov., southeastern Qinghai Prov. and eastern Tibet at high mountain areas [3,4]. It is common at high altitudes (3000-4300 m) and exists in coniferous forests, mixed broadleaf-conifer forests and alpine scrubs [4]. C. auritum is only encountered in the mountainous regions of Qinghai, Gansu and Sichuan provinces and Ningxia Hui Autonomous Region [1].
Previous studies on Crossoptilon species have focused on their distribution, ecology and behaviour. Some studies have used single gene or partial sequences in reconstructing their phylogenetic relationships [15], but a comprehensive study of the entire mitogenomes has rarely been attempted. To reconstruct the phylogenetic relationships among Crossoptilon and other Phasianidae species, the entire mitogenomes of C. harmani, C. mantchuricum and C. crossoptilon were sequenced, and thus four species of the genus Crossoptilon were thoroughly compared. Furthermore, the estimated divergence times of the four Crossoptilon species were calculated, and a Ka/Ks analysis was used to estimate the adaptation of the Crossoptilon mitogenome to different environments.

Sample collection and DNA extraction
Samples of C. harmani were collected from the Nagqu area in Tibet, China, and voucher specimen was deposited in Institute of Zoology, Shaanxi Normal University; C. mantchuricum and C. crossoptilon were collected from the Beijing Zoo, and samples were obtained from bird specimen collections of the National Zoological Museum, Institute of Zoology, Chinese Academy of Sciences. The collection was under the permit from Forestry Department and conformed to the National Wildlife Conservation Law in China. No living animal experiments were conducted in the current research. All the samples were preserved in 100% ethanol and stored at −20°C. The total genomic DNA was extracted from the liver/muscle tissue using the standard phenol/chloroform method [16].

PCR amplification and sequencing
The PCRs were performed under the following conditions, 2 min initial denaturation at 93°C, 40 cycles of: 10 s denaturation at 92°C, 30 s annealing at 58-53°C, 10 min elongation at 68°C in the preliminary 20 cycles, and 10 s denaturation at 92°C, 30 s annealing at 53°C, and 10 min elongation at 68°C with 20 s per cycle added to the elongation step in the succeeding 20 cycles, and finally an extension for 7 min at 68°C. The amplifications were performed in 15 μL reactions containing 2.4 μL of 2.5 mM dNTPs, 2.1 μL of each primer at 10 μΜ, 1.5 μl of 10× LA PCR Buffer I (Mg 2+ -free), 1.5 μL of 25 mM MgCl 2 , 1 μL of DNA template, 0.18 μL of 5 U/μl LA Taq polymerase (Takara, Dalian, China) and 4.22 μL ddH 2 O.
The sequences of the primers used for PCR amplification and sequencing of the mitochondrial genes were obtained from Sorenson (2003) [17] with minor changes (Additional file 1). The mitogenome of C. harmani was amplified in seven parts and the gaps were bridged using other adjacent primers. The PCR products were purified using the DNA Agarose Gel Extraction Kit (Bioteke, Beijing, China) after separation by electrophoresis on a 0.8% agarose gel. After separation and purification, the PCR products were sequenced by Sangon Biotech (Shanghai) Co., LTD. using the Primer-Walking method. The C. mantchuricum and C. crossoptilon mitogenomes were amplified in 12 or 13 fragments, and sequencing was performed using the Illumina Hiseq2000 high-throughput sequencing system of Shenzhen Huada Gene Technology Co., LTD. The gaps in the assembly after high-throughput sequencing were filled in by direct sequencing using the ABI 3730 DNA sequencer by Sangon Biotech (Shanghai) Co., LTD., using adjacent PCR primers.

Gene identification and genome analyses
The Staden sequence analysis package [18] was used for sequence assembly and annotation of C. harmani mitogenome. The complete genome assemblies for the mitogenomes of C. mantchuricum and C. crossoptilon were performed using the SOAP de novo software. Most tRNA genes were identified using tRNAscan-SE 1.21 [19] under the 'tRNAscan only' search mode, with the vertebrate mitochondrial genetic code and 'mito/ chloroplast' source. The protein-coding genes (PCGs), rRNA genes and the remaining putative tRNA genes that were not identified by tRNAscan-SE were identified by sequence comparison with other Galliformes species. The rrnS secondary structure of Crossoptilon was predicted based on the structure of Gallus gallus and Anas platyrhynchos obtained from the Comparative RNA Web (CRW) [20], and Pseudopodoces humilis (now as Parus humilis) structure [16]. The rrnL secondary structure was predicted based on the structure of Xenopus laevis obtained from the CRW database, Bos taurus [21] and P. humilis structures [16]. The RNAstructure software was used to identify and draw potential secondary structures in the single-stranded control region. The nucleotide compositions of the mitogenomes and amino acid information were analysed using MEGA 4.1 [22].

Sequence alignments
Along with the entire mitogenomes obtained in this study, 42 Galliformes sequences were used in the phylogenetic analysis, including two outgroups (Numida meleagris and Alectura lathami). The DNA sequences of the other species used in the phylogenetic analyses were downloaded from GenBank (the accession numbers and key information are shown in Additional file 2). The tRNA and rRNA genes and the CR were individually aligned using ClustalX 1.83 [23] with the default settings. All 13 protein-coding genes were translated into amino acids, and then aligned using MEGA 4.1 [22] with default parameters for each gene, and finally retranslated into nucleotide sequences.

Phylogenetic analyses of Phasianidae
Datasets containing 13 protein-coding genes (PCG) and all 37 genes plus the control region (mitogenome) were used to study the phylogenetic relationships within Phasianidae. Phylogenetic analysis based on nucleotide sequences was performed using PAUP*4.0b10 for the Maximum parsimony (MP) method [24], RAxML-7.0.3 for maximum likelihood (ML) [25] and MrBayes 3.1.2 for Bayesian inference (BI) [26]. For ML and BI analyses, models of the concatenated nucleotide sequences datasets were assessed independently using AICc in MrMo-deltest2.2 [27]. The best fit model GTR + I + G was chosen for the likelihood and Bayesian analyses. A consensus tree was generated for MP analysis under the majority rule. The reliability of the clades in the phylogenetic trees was assessed by bootstrap probabilities (BSP) computed using 1000 replicates, with random addition for each bootstrap replicate. The 1000 replicates bootstrap support was also performed in the ML analysis. Bayesian analysis with Markov Chain Monte Carlo sampling was run for 1000000 generations saving a tree every 100 generations, with one cold and three heated chains, and the burn-in time was determined by the time to convergence of the likelihood scores. The Bayesian posterior probabilities (BPP) were estimated on a 50% majority rule consensus tree of the remaining trees.
The MEGA 4.1 [22] was used to calculate the pairwise genetic distance for four Crossoptilon species with default parameters. The mitogenome aligned data and four single genes (nad2, CR, cytb and rrnS) obtained from GenBank (Additional file 3) were used to calculate the genetic distances. These genes were aligned singly, and be adjusted to consistent sequence lengths manually.

Divergence time estimates focused on Crossoptilon
Along with the PCG dataset obtained in this study, 42 Galliformes sequences were used to estimate the divergence time of the Crossoptilon species. The divergence time of Crossoptilon species was calculated using the Bayesian procedure implemented in BEAST v. 1.7.2 [29,30]. A relaxed clock was used with rates complying with a log-normal distribution [31]. The GTR + I + G model and a Yule prior were used in the analysis. The calibration points were based on the fossil records showing that stem Numididae-Phasianidae split at 50-54 Mya (million years ago) [32]; Arborophila rufipectus diverged from the other lineages in the Galliformes around 39 Mya [33]; Coturnix-Gallus split at 35 Mya [34][35][36]. The results of runs of 10 million generations were used after a burn-in of 100.

Ka and Ks analysis
To better understand the evolution at the DNA level and the role of selection in the four Crossoptilon species, we calculated the nonsynonymous and synonymous substitution rates using Kaks_calculator 2.0 [37] for six groups [C. harmani-C. mantchuricum (C.har-C.man), C. mantchuricum-C. crossoptilon (C.man-C.cro), C. harmani-C. crossoptilon (C.har-C.cro), C. harmani-C. auritum (C.har-C.aur), C. mantchuricum-C. auritum (C. man-C.aur), and C. crossoptilon-C. auritum (C.cro-C. aur)]. The ratio of nonsynonymous substitution rate (Ka) to synonymous substitution rate (Ks) is widely used as an indicator of selective pressure at the sequence level among different species. It is commonly accepted that Ka > Ks, Ka = Ks, and Ka < Ks generally indicate positive selection, neutral mutation, and negative selection, respectively [38,39]. To calculate Ka, Ks and Ka/Ks, a model averaging method was selected. This method includes 14 different models for calculation and derived the average values for Ka, Ks, and Ka/Ks [37]. The genetic code selected was the 'vertebrate mitochondrial code'. To further study the selective pressure acted on each protein-coding gene in the genus Crossoptilon, CodeML in PAMLX software [40] was used to find sites under strong selective pressure. The secondary structure analysis of amino acid was performed by using an online software TOPCONS [41].

Results
Comparison of the two sequencing methods used in this study revealed that high-throughput sequencing has greater coverage and accuracy relative to standard sequencing, albeit with higher cost. Using standard sequencing, we obtained 62 effective sequences with 2.63-fold coverage. In contrast, the high-throughput sequencing yielded effective assembly data with sequence depths (X) of 7604.35 for C. mantchuricum and 7810.82 for C. crossoptilon after filtering out some reads, such as low-quality or adapter-sequence-polluted reads.

Protein-coding genes
The 13 protein-coding genes of Crossoptilon genomes are similar to most other Phasianidae species, with nad5 and atp8 being the longest and shortest genes, respectively. The total length of the PCGs in each Crossoptilon species is 11358 bp after removing termination codons, containing approximately 3786 codons. The A + T content of the 13 PCGs is 53.4% in C. harmani and C. crossoptilon, and 53.3% in C. auritum and C. mantchuricum. Analysis of the base composition at each codon position of the concatenated PCGs shows that the second codon position has a higher A + T content (57.9% in C. auritum and C. mantchuricum, 58.1% in C. harmani, and 58.0% in C. crossoptilon, respectively) than the first and third codon positions. The amino acid frequencies in the Crossoptilon PCGs are similar, with Leu significantly more frequent than other amino acids. The different codon positions have the same base distributions, with C and G as the most and least frequent bases in the third codon, respectively.
All PCGs in Crossoptilon species start with the typical ATG codon, except the nad5 gene in C. mantchuricum and the cox1 gene in four Crossoptilon species, which start with GTG. Four types of stop codons are used by the coding genes, including TAA and TAG for most genes, AGG for cox1 in Crossoptilon and nad6 in C. harmani and C. crossoptilon, and an incomplete stop codon T-for cox3, nad4 and nad2 in four Crossoptilon species, respectively.

RNA genes
Similar to previously sequenced mitogenomes, the genomes sequenced in this study contain 2 rRNA genes encoding the small and large rRNA subunits, which are located between trnF and trnL(UUR) and separated by the trnV gene. The lengths and A + T contents of the rrnS and rrnL in the Crossoptilon genomes are within the range observed in other Phasianidae species. The rrnS contains three domains with 46 predicted stems, and the rrnL contains six domains with 59 stems (Additional files 4 and 5). The secondary structures of the rrnS in C. crossoptilon and C. harmani are identical; the rrnL secondary structures differ by only 2 bp in length. However, there are many differences in the rrnL secondary structures among the four Crossoptilon species, specifically in the loop near stem 44, which contains several substitutions and indels.
The complete mitogenome sequence contains 22 interspersed tRNA genes. All the tRNA sequences have the potential to fold into typical cloverleaf secondary structures except for trnS(AGY), which lacks the DHU arm (Additional file 6). The secondary structures of trnF, trnL(UUR), trnQ, trnM, trnC, trnY, trnS(UCN), trnH, trnL(CUN), trnT, and trnP are relatively conserved. However, the structure of trnS(AGY) is different in C. auritum relative to the other three Crossoptilon species (Additional file 6); the structure in C. harmani, C. mantchuricum and C. crossoptilon contain two additional bases (G and A) in the amino acid acceptor arm. The most frequent mismatch is G-U; other mismatches are also present, including U-U in trnM and trnG, C-C in trnD and trnL(UUR), C-U in trnF and trnI, and A-C in trnS(AGY). The C. crossoptilon and C. harmani tRNA structures are almost identical, with the exception of 1 bp. The A + T content of the tRNA genes is 57.6% in C. harmani and C. crossoptilon, 57.4% in C. mantchuricum and C. auritum.

Control region
The nucleotide composition of the control region in the Crossoptilon species has a bias against G, which is common in the mitogenome sense strand of vertebrates [42]. The control region of Crossoptilon is located in the conserved position between trnE and trnF and the length (1146 nucleotides) is conserved in all four species. The control region contains three domains: the ETAS Domain I (nt 1-312), Central Conserved Domain II (nt 313-780) and CSB Domain III (nt 781-1146). The A + T content of the three domains is similar in all the Crossoptilon species; Domains I, II, and III have higher contents of C and A, T, and A, respectively, relative to other bases. Domain III contains a higher percentage of A than Domain II, and Domain II has the highest G content among the three domains. The distribution of variable sites and conserved sites suggests that Domain II has relatively more conserved sites and Domain I has more variable sites compared to other domains.
The ETAS Domain I can be divided into two parts: part A, from nt 1-163 and part B, from nt 164-312.
There are two conserved blocks in part A, ETAS1 (nt 64-126) and ETAS2 (nt 124-163), which are similar to motifs previously identified in other avian and mammalian species. The first block is perfectly conserved among the Crossoptilon species and has sequence similarity to the "goose hairpin" described in some phasianids and Anas species [43][44][45][46] (Additional file 7). The secondary structure of this hairpin is determined by a stem of seven complementary Cs/Gs and a loop containing a TCCC motif also present in mammalian control region [47], which was associated experimentally with termination of H strands [48]. There are two copies of TCCC located at nt 22-25 and 183-186. The second block, which is perfectly conserved among the Crossoptilon species, has sequence similarity to the mammalian TASs, including the highly conserved motif GTGCAT, which is present in all sequenced Phasianids and Anseriforms. The motif GYRCAT (Y = C/T; R = A/G) is widespread in Domain I of some mammalian control region [47], and it has been duplicated in the R1 repeats of many species, including the opossum [49] and several rodents [50,51]. Its functional importance is suggested by both comparative [47] and experimental [48] data.  [43,52], and a bird similarity box is also present.
The CSB Domain III is highly variable and has sequences similar to mammalian CSB1. A poly(C) sequence (nt 781-792), similar to the O H of mammals, maps just a few nucleotides downstream from the putative CSB1 (nt 803-828). However, it is difficult to identify sequences corresponding to mammalian CSB2 and CSB3. The secondary structures of CSB1 and the putative goose hairpin in Crossoptilon are consistent with those of G. gallus (Additional file 7). The bidirectional LSP/HSP promoters (nt 982-1003) [53] are almost perfectly conserved among the Crossoptilon species. A stable hairpin (nt 1004-1017) that is rich in poly (T) and poly (A) strings is located immediately upstream of the promoters.

Phylogeny and divergence time of the Crossoptilon species
The final combined PCG dataset has 11376 characters after alignment. For this dataset, parsimony analysis shows a length of 27412 steps, with CI =0.313, RI = 0.499. The topologies between ML and BI trees of PCG dataset were consistent (Figure 1). For the mitogenome dataset, there was no difference in topology between the ML and BI trees of Galliformes (Additional file 8); however, both trees differed from the MP tree (Additional file 8). According to the phylogenetic results, the monophyly of Crossoptilon was strongly supported in MP, ML and BI analyses. Within Phasianidae, the topology ((Crossoptilon + Lophura) + (Phasianus + Chrysolophus)) was formed in most trees. The sister-group relationship between Crossoptilon and Lophura was supported with bootstrap values 69 and 100 in MP and ML trees, and posterior probabilities 1.00 in BI tree of PCG dataset.
The genetic distances of the mitogenomes (Additional file 9) are identical (0.026) for C. auritum and C. harmani, C. auritum and C. crossoptilon, C. harmani and C. mantchuricum, and C. mantchuricum and C. crossoptilon. However, the genetic distances between C. harmani and C. crossoptilon, and C. auritum and C. mantchuricum are only 0.001 and 0.002, respectively. The genetic distances based on four single genes show that the genetic distances of C. harmani-C. crossoptilon and C. auritum-C. mantchuricum are much smaller. This value is similar to the intraspecific level and far less than the interspecific level (Additional file 9).
According to the estimated timescale obtained from the phylogenetic tree containing consistent topology with PCG-ML/BI trees, the C. harmani-C. crossoptilon and C. mantchuricum-C. auritum splits occurred at 3.

Nonsynonymous and synonymous substitution
The analysis of variable sites in each protein-coding gene showed that the group C.har-C.cro and C.man-C.aur contained less variable sites, while other four groups included more (Additional file 10). The percentages of variable sites were higher in nad3 and nad6 genes among the groups C.har-C.man, C.man-C.cro, C.har-C. aur and C.cro-C.aur, while the percentages were lower in atp8 and atp6 genes. The Ka/Ks value of atp8 in the five groups is far greater than 1 except for C.har-C.cro, and atp6 and cytb in C.har-C.cro (Table 2), which show a strong positive selection. The Ka/Ks values of the other genes analysed in Crossoptilon species were less than 1, which shows a purifying selection ( Table 2). Furthermore, the P-value (Fisher exact test) is much less than 0.001, except for atp8 in C.har-C.man, C.man-C.cro, C. har-C.aur and C.cro-C.aur (0.558274 or 0.289204, respectively, which are obviously higher than 0.05) and nad3 in C.har-C.cro (0.154), indicating that the difference is significant.
Both atp6 and cytb genes have one different base in the C.har-C.cro group, i.e., nt 245 in atp6 gene (T base or Phe in C. harmani, while C base or Ser in C. crossoptilon) and nt 709 in cytb gene (C base or Leu in C. harmani, while T base or Phe in C. crossoptilon). The atp8 gene in C.man-C.aur group also contains one different base (nt 124), G base in C. mantchuricum, while A base in C. auritum, and conrespondingly Val in C. mantchuricum, while Met in C. auritum. Further study showed that four amino acid sites in atp8 (T at 11 position, I at

PBS analyses
PBS analyses were performed to better understand the contributions of different parts of the mitogenomes to the genome phylogeny based on the mitogenome-ML tree. Figure 1 The phylogenetic trees based on the PCG dataset using ML and BI methods. Branch lengths and topologies were obtained using Bayesian inference analyses. Among the Crossoptilon, C. auritum with C. mantchuricum and C. harmani with C. crossoptilon form two pairs of sister groups.
The relationships among the PBS value, length, singleton sites (S), parsimony informative sites (Pi), variable sites (V) and conserved sites (C) in different partitions are also studied. Ranking individual protein-coding genes by their respective contribution to the total PBS values shows that some genes, such as nad5 and nad4 provide a higher contribution compared to other markers, while atp8 and nad4L contribute less. But analysis also shows that PBS values are roughly correlated with gene lengths. The variable sites are closely related to parsimony informative  sites, which follow the same trend. The third codon has a high PBS value but with less conserved sites, while the tRNA gene cluster IQM has a lower PBS value with more conserved sites compared to other partitions.
The third codon has a high PBS value and therefore contributes most to the mitogenome-ML tree. The NADH genes have secondary high PBS, which might be the effect of the third codon in the relatively long sequence. To further understand the contribution of the 3rd codon, we reconstructed the phylogenetic tree using this single partition.
In the phylogenetic trees based on the third codon (Additional file 8), the structure of the ML tree is similar to the BI with the exception of Phasianus. The positions of Pucrasia macrolopha, Bonasa bonasia, Meleagris gallopavo and Polyplectron bicalcaratum in the 3rd-MP tree are different compared to the BI tree. The 3rd-BI tree is identical to the mitogenome-ML with relative high support except for the position of Polyplectron bicalcaratum.

Discussion
Phylogenetic relationship and divergence time of the Crossoptilon species In this study, the phylogenetic trees based on different datasets support the monophyly of Crossoptilon and the close genetic relationships of C. auritum and C. mantchuricum [4,12,54,55], C. harmani and C. crossoptilon [1,4,8,14,55]. Our analysis of the gene sequences of Crossoptilon based on different datasets confirms that Crossoptilon is the sister group to Lophura [55,56], indicating their closer relationship.
The genetic distance is rather small among some Crossoptilon species, for instance, between C. harmani and C. crossoptilon, C. auritum and C. mantchuricum. Morphological, ecological and behavioural studies have revealed the high similarity between C. harmani and C. crossoptilon, but with different plumage coloration [1]. Their distributions are overlapping, and C. harmani routinely hybridises with C. c. drouynii [1], which indicates their reproductive compatibility. Furthermore, the distribution of C. harmani is limited between the Himalayan and Nyenchen Tanglha ranges, which provided the geographical isolation required for the formation of a subspecies or species. Therefore, the taxonomic status and its relationship to C. crossoptilon are much questionable, which need to be re-evaluated by multiple markers and population genetics in the further studies.
Based on their significant differences in many aspects such as morphology and behaviour, C. auritum and C. mantchuricum are commonly considered to be independent species. However, the genetic distance of 0.002 indicates that they have a much closer genetic relationship. Consistent with this observation, Tsam et al. (2003) [12] observed that the interspecific differentiation between C. mantchuricum and C. auritum (0.18%) is less than the degree of differentiation among the C. crossoptilon subspecies.
Some Phasianidae species were found in the Pliocene epoch, and several modern taxa had already appeared in the Quaternary Pleistocene. The earliest fossils of Crossoptilon were found in the Cenozoic Zhoukoudian strata of Beijing, and in the Pleistocene strata at Yanjinggou in Wanxian of Sichuan Prov. [57,58]. The evolution of Crossoptilon mainly occurred corresponding to the Tibetan Plateau intensive uplift during the Tertiary Pliocene and the alternation of the glacial and interglacial periods in the Quaternary, which created profound and complex changes in the geographical and ecological environments [59,60] and was thought to have greatly affected its topography and avian species diversification [61,62]. Based on our results, the spliting time between C. harmani-C. crossoptilon and C. mantchuricum-C. auritum is consistent with the uplift during the Tertiary Pliocene and the fossil record in the Sichuan Wanxian Yanjinggou Pleistocene strata, which is approximate to the divergence times estimated by Jiang et al. (2014) [63] (3.78 Mya, 95% HPD = 1.17-6.56 Mya). Our study shows that C. mantchuricum and C. auritum are relatively ancient groups, while C. harmani and C. crossoptilon diverged later; these observations are not consistent with the Crossoptilon evolution hypothesis proposed by Lu et al. (1998) [1]. Lu et al. (1998) [1] suggested that the ancestor of C. crossoptilon and C. harmani originated in Sichuan Prov., Yunnan Prov., and the Tibetan border area. They noted that the clade containing the C. harmani plesiomorphy migrated to the hinterland plateau and subsequently differentiated into the ancestor of C. harmani and C. auritum. The ancestor of C. auritum dispersed in the northern plateau, and differentiated into C. auritum and C. mantchuricum. Based on the divergence time and distribution characteristics of Crossoptilon species, our results support the ancestor of Crossoptilon was initially distributed in Sichuan Prov., Yunnan Prov., and the Tibetan border area. The Crossoptilon first diverged into the ancestor of C. harmani-C. crossoptilon and the ancestor of C. auritum-C. mantchuricum. Considering the geographical distribution of Crossoptilon [1,[64][65][66], we propose that the ancestor of C. harmani-C. crossoptilon might migrated to platform of the plateau, and then splitted into two lineages during the last uplift of the Tibetan Plateau [67], one adapted the high altitude environment to form C. harmani in east Tibet, and the other evoluted into C. crossoptilon, while the ancestor of C. auritum-C. mantchuricum dispersed to the northern plateau and further differentiated nearly the same period with C. harmani-C. crossoptilon.

Selective pressure on protein-coding genes in Crossoptilon species
The Ka/Ks values in most protein-coding genes are less than 1 (Ka is lower than Ks) with P-values less than 0.001, which indicates that these genes in Crossoptilon are under purifying selection. The Ka/Ks values of each gene among C.har-C.man, C.man-C.cro, C.har-C.aur and C.cro-C.aur are similar. The Ka/Ks value of the atp8 gene is 50 in five groups except for C.har-C.cro, which reveals the strong positive selection. However, only the P-value for C.man-C. aur is far less than 0.01 (P = 0), which indicates significant differences. In contrast, the Ka/Ks values of the atp6 and cytb genes in C.har-C.cro are 50 with P-values less than 0.001, which also indicates strong positive selection. Previous studies have shown that atp6 gene is highly variable in the mitogenomes of humans living in extremely cold areas [68][69][70][71], the atp6 gene faces strong selection pressure with increasing altitude [72], and the variation of atp6 gene in Artemia tibetiana is the result of adaptation to the cold and hypoxic environment of the plateau [73]. The variation in atp6 may reflect the adaptation of C. harmani to the plateau environment; cytb might also play an important role in this process. The atp8 gene in C. mantchuricum and C. auritum may also have experienced a strong selection in adaption to the plateau environment.
However, the further study showed that only four amino acid sites of atp8 in the Crossoptilon are positively selected sites, and accordingly this gene may have experienced strong selective pressure for plateau adaptation. The corresponding TM-helix (IN-> OUT) or outside position might be important for the Crossoptilon adaptation. After discarding groups (C.har-C.cro and C.man-C. aur) and genes (cox1, atp8 and nad4l) with obviously different Ka/Ks values to others, the average value (0.0431) is similar with the strongly locomotive group (0.04) in the study of Shen et al. (2009) [74], which indicates the Crossoptilon accumulates fewer nonsynonymous mutations and is corresponding to their strong running ability.

Conclusions
In summary, the mitogenomes of the Crossoptilon species contain the same gene arrangements and similar compositions including base contents and secondary structures. According to the phylogenetic trees, C. auritum with C. mantchuricum, and C. harmani with C. crossoptilon form two pairs of sister groups. Crossoptilon have a closer relationship with Lophura. Based on the genetic distances, C. harmani is more closely related with C. crossoptilon, and is the most recent diverged descendent of Crossoptilon as a result of the plateau adaptation. According to the molecular dating results, the divergence time between C. harmani-C. crossoptilon and C. mantchuricum-C. auritum is consistent with the uplift of the Tibetan plateau and the subsequently climate change during the Tertiary Pliocene. The Ka/Ks analysis showed that atp8 gene in the Crossoptilon species may have experienced strong selection for the plateau environment adaptation.