Comparative phylogeography of two sister species of snowcock: impacts of species-specific altitude preference and life history

Abstract

Background

Phylogeographical patterns and population dynamics are usually interpreted by environmental disturbances and geographic barriers of the past. However, sister species may exhibit disparate patterns of genetic structures and population dynamics due to their habitat preference and altitude segregation. In this study, we tested how species-specific altitude habitat affected phylogeographical patterns in two sister snowcock species, Tibetan (Tetraogallus tibetanus) and Himalayan Snowcocks (T. himalayensis).

Methods

A panel of seven microsatellite loci and a fragment of Mitochondrial DNA Control Region were used to investigate genetic structures and population dynamics in hope of revealing the underlying evolutionary processes through the identification of possible past demographic events.

Results

Our results suggest that T. himalayensis showed a significant phylogeographical signal in mtDNA (F_ST = 0.66, p < 0.001) and microsatellite (F_ST = 0.11, p < 0.001) data and is stable during the glacial-interglacial cycles in the Pleistocene and followed demographic contraction until 0.003 million years (Mys) ago. The phylogeographical signal of T. tibetanus is lower than the level of genetic difference among populations in mtDNA (F_ST = 0.41, p < 0.001) and microsatellite (F_ST = 0.09, p < 0.001) data, likely benefiting from stable habitats over a long period of time. T. tibetanus has been experiencing expansion since 0.09 Mys ago. However, an abnormally haplotype H9 from T. himalayensis clustering with T. tibetanus was spotted.

Conclusion

Our results indicate that differences in habitat preference and altitude specialities were reflected in the genetic structure patterns and population dynamics of these two species. These dissimilarities in life history traits might have affected the dispersal and survival abilities of these two species differently during environmental fluctuations. The results of this study also enriched our knowledge on population differentiation and connectivity in high altitude mountain ecosystems.

Background

Montane glacial shifts and environmental changes during the Pleistocene are the main reasons to phylogeographical patterns and population dynamics (Fang et al. 2017). However, different phylogeographical patterns in closely related species might have been caused by different habitat preferences and altitude specificity (Qu et al. 2009). Especially, sister species, which shared recent ancestry, may show different patterns of population genetic structure due to biological, ecological, behavioural, or environmental factors (Bermingham and Moritz 1998). Hence, comparative analysis of the population genetic structures of sister species may provide the chance to disclose factors impacting evolutionary processes. Comparative phylogeography can reveal concordant phylogeographical patterns across sympatrically distributed species, demonstrating a relative historical stability of communities and reflections to resemble external (e.g. climatic) factors or incongruent patterns, showing a prevalence of particular, species-specific reflections (Correa Ribeiro et al. 2016).

The Qinghai–Tibetan plateau (QTP) is one of the most prominent topographical features on Earth. There is a growing understanding that avian evolution in the QTP region has been firmly associated with the orogenetic history of the flanking mountain chains (An et al. 2009, 2015; Qu et al. 2009; Lei et al. 2014). The Tibetan Plateau stretches from the Himalayan Mountains in the south to the Kunlun, Altun and Qilian Mountains in the north and from the Pamirs and the Hindu Kush in the west to the western part of the Qinling Mountain in the east (Zhang et al. 2018). However, the investigation on origin and evolutionary history of Tibetan endemic bird species has long been limited to the central plateau region, southern and eastern fringes (Prum et al. 2015). Although the biogeographical avian history of the western part of Tibetan Plateau has little information, research on the Snow Leopard (Panthera uncial) showed a clear biogeographical pattern in the distribution of three subspecies with Northern (the Altai region), Central (core Himalaya and Tibetan Plateau), and Western (Tien Shan, Pamir, Janecka et al. 2017). Previous research also found that differences among species in body size, shape as well as in foraging and feeding strategies had been involved in the initiative stages of SE Asian passerine radiations long before ancestors of extant species occupied their altitudinal niches in response to Sino-Himalayan vegetation belts (Kennedy and Price 2012; Price et al. 2014).

Tibetan Snowcock (Tetraogallus tibetanus) and Himalayan Snowcock (T. himalayensis) inhabit the alpine and subalpine zone in the remote and rugged mountains on the Tibetan Plateau and its surrounding mountain ranges (Zheng 1978). They well adapt to the harsh climate, spending the hostile winter months in lower valleys and mostly moving in summer to high elevations for breeding as resident birds (An et al. 2015). They evolved in large mountains (Shen and wang 1963) and overlapped distribution in western Himalayan Mountains, Qilian Mountains, Qingha Nanshan and Kunlun Mountains (Liu 1999), with different ecological niches due to the body size (big vs small), habitat altitude (high vs moderately high), breeding time (late vs early) and diet spectrum (euryphagous vs stenophagous; Additional file 1: Table S1) (Liu 1999). T. tibetanus has a much less fragmented montane habitat and more cold tolerance than T. himalayensis. Based mainly on the head, abdomen, and adult plumage color, T. himalayensis and T. tibetanus were separately divided into dark-bellied and white-bellied groups (Bianki 1898). The voices of T. tibetanus and T. himalayensis are also distinct (Adeli et al. 1997). We hypothesize that the phylogeographical patterns and/or population dynamics of the two sister snowcock species may be different due to their different life history and habitat preferences, even though they share the same historical Pleistocene glacial fluctuations and environmental conditions. Habitat generalists with more ecological plasticity and broad geographic distributions have better chances to survive when in adverse conditions than specialists with specific habitats and narrow ranges (Zhang et al. 2012). Climate in high altitude mountains is usually cold and moist, and coniferous forests cover the mountains and higher elevation valleys, while low altitude foothills are warmer and drier (Zhang et al. 2013).

Snowcocks can survive only around the snow line or in low-temperature surroundings. Contrary to most other avian species, individuals spread only in glaciers and are isolated in interglacial locations (Liu 1999). Li et al. (2016) identified climate refugia of Snow Leopards due to the unique mountain environment, where snowcocks sympatrically distributed in High Asia, which maintains a relatively constant arid or semi-arid climate. With world climate changing, Snow Leopards’ suitable zones will move up towards mountaintops, causing the loss of Snow Leopard habitat and habitat fragmentation along (Li et al. 2016).

Although Ruan et al. (2005) and Liu (1999) documented that the two species never coexist on a hill in overlapped regions due to strong competition, an intriguing phenomenon was observed that the two species were mixed together based on recent years’ field observation in Qinghai–Tibetan Plateau (Ma et al. 2013; Zhang et al. 2015), and sympatrically distributed in Ladakh, Jammu and Kashmir, India (Namgail 2005). In this study, we examined the phylogeographical patterns and population structure and inferred the evolutionary history of T. himalayensis and T. tibetanus. The comparison of the two sister species, with different ecological habitat preferences, and altitude specificity, will offer more insights into the interaction between historical population dynamics and geological events in this region.

Methods

Sampling and DNA extraction

A total of 25 sample sites covered the distribution ranges of T. himalayensis and T. tibetanus on the QTP (Fig. 1). According to T. tibetanus geographical origin, the samples were divided into five regional groups: QLS (Qilian Mountain region), QDM (Qiadam Basin region), BKL (Baryan Har Mountains region), TGL (Tanggula Mountains region) and WKL (western Kunlun Mountains region) and according to T. himalayensis geographical origin, the samples were divided into six regional groups: HETS (eastern Tien Shan region), HWTS (western Tien Shan region), HWKL (western Kunlun Mountains region), HQDM (Qiadam Basin region), HQLS (Qilian Mountain region), HPME (Pamir Plateau) (Fig. 1 and Additional file 1: Table S2). The sample size of each population is shown in Table 1. Voucher specimens are held in the Lanzhou University Nature Museum. The research was conducted under permission, by the Forest Department of Gansu Province (China, Gansu), Xinjiang (China, Xingjiang Autonomous Region), Xizang (China, Xizang Autonomous Region). DNA was extracted using whole genome kit (Qiagen, Inc., Valencia, CA) except for some individuals of T. himalayensis specimens where DNA was already available (Wang et al. 2011). One hundred and nine T. himalayensis samples were genotyped at seven polymorphic microsatellite primers that were isolated originally from the Chicken (Gallus gallus) genome: (MCW135, MCW323, AB121114, GUJ0028, UBC0005, UBC0006 and GUQ0001).

Fig. 1

Distributions and sample locations for the two snowcock species: a Tetraogallus himalayensis, marked with blue rectangular solid, b T. tibetanus, marked with pink triangular solid. Sample sites are abbreviated as follows: Datong (DT), Sunan (SN), Tianzhu (TZ), Atushi (AT), Wuqia (WQ), Yecheng (YC), Haixi (THX), Delingha (TDL), Zhiduo (ZD), Biru (BR), Baqin (BQ), Anduo (AD), Suoxian (SX), Qumalai (QM) for T. tibetanus sampling sites, Subei (SB), Pishan (PS), Hetian (HT), Kashi (KS), Germ (GM) Delinha (HDL), Qitai (QT), Houxia (HHX), Akesu (AK), Zhaosu (ZS), Taxkorgan (TS) for T. himalayensis. Group definitions see Additional file 1: Table S2. QLS (Qilian Mountains region), QDM (Qiadam Basin region), BKL (Baryan Har Mountains region), TGL (Tanggula Mountains region) and WKL (western Kunlun Mountains region) for Tetraogallus tibetanus, and for T. himalayensis geographical origin, the samples were divided into six regional groups: HETS (eastern Tien Shan region), HWTS (western Tien Shan region), HWKL (western Kunlun Mountains region), HQDM (Qiadam Basin region), HQLS (Qilian Mountain region), HPME (Pamir Plateau)

Table 1 Genetic diversity of populations of T. tibetanus and T. himalayensis, with mitochondria DNA control sequences (mtDNA) and 7 nuclear microsatellite markers (STR)

PCR was performed in 50 μL reaction mixtures containing 1.5 mM of MgCl₂, 0.6 μM of primer pair, 0.4 mM dNTP, 100 ng of the DNA template, 0.4 U of Taq DNA polymerase (Sigma) and 1× PCR buffer (Sigma). The PCR protocol used was: (1) initial DNA denaturation at 94 °C for 4 min; (2) 35 successive cycles of strand denaturation at 94 °C for 40 s; (3) primer annealing at 45‒60 °C for 40 s and (4) DNA extension at 72 °C for 45 s. A prolonged extension step at 72 °C for 6 min was added after the final cycle. Amplified products, which ran on 8% polyacrylamide gel by electrophoresis, could be visualized by silver staining. PUC19 DNA/MspI (HpaII), as a size standard, ran on each gel to determine fragment size using Bandscan 4.30 software (http://moleco.sjtu.edu.cn).

Genetic diversity and phylogeographic structure

All 67 T. himalayensis mitochondrial control region sequences from Wang et al. (2011) were aligned with 80 sequences of T. tibetanus in An et al. (2015). We identified haplotypes and calculated standard genetic indices (haplotype diversity (h) and nucleotide diversity (π)) for each population and combined diversity for all populations using DNASP 5.0 (Librado 2009). The sequence of Alectoris magna was used as outgroup. Evolutionary relationships among all haplotypes were performed by constructing phylogenetic tree using Bayesian (BA) Markov chain Monte Carlo (MCMC) phylogenetic analyses (Drummond and Rambaut 2007). The computer program jMODELTEST 2.1.1 (Darriba et al. 2012) was used to select the best-fit model of evolution. For Bayesian analyses, four independent MCMCs were initiated with random starting trees and each run for 5 × 10⁶ generations, sampling every 100 generations. A 50% strict consensus tree was computed for the 1000 bootstrap trees. We further constructed unrooted haplotype networks using the median-joining algorithm (Bandelt et al. 1999) as implemented in Network v4.6.1.1 (http://www.fluxus-engineering.com). This method allows the visualization relationships and frequencies among mtDNA haplotype.

All T. himalayensis microsatellite data and fragment lengths were analysed with the internal size marker GENESCAN-500 ROX (Applied Biosystems), and scored with GENEMARKER 3.7 (SoftGenetics). The published microsatellite data (An et al. 2009) were involved in this study. For microsatellite data, we estimated the standard genetic diversity for each group using the following indices: the number of alleles at each locus (N_A), average allelic richness (A_R), and observed heterozygosity (H_O); heterozygosity expected from the Hardy–Weinberg assumption (H_E) and exact tests of linkage disequilibrium (LD) between pairs of loci for each population were calculated with ARLEQUIN 3.5 (Excoffier and Lischer 2010) and FSTAT 2.9.4 (Goudet 2001). STRUCTURE 2.3.4 (Pritchard et al. 2000) was used to determine the most probable number of genetic clusters (K). Analyses were performed using an admixture model, correlated allele frequencies, a burn-in period of 50,000 cycles, and 500,000 additional cycles (determined from test runs to be enough for parameter stabilization). The admixture model with correlated allele frequencies was applied and with sampling locality priors for magnifying the signals as suggested by Hubisz et al. (2009). Analyses were repeated 20 times for each species with K = 1 to 5 or 6, and the two species together with K = 1 to 11, where K is the number of genetic populations, and posterior probability, ln[P(D)], was used to infer the most likely number of genetic populations as described in Pritchard et al. (2000). We used STRUCTURE HARVESTER to process the results from multiple runs of STRUCTURE (Earl and vonHoldt 2012).

Genetic differentiation between regional groups was evaluated based on pairwise values of F_ST. The statistical significance of the estimates was assessed after 10,000 permutations. Gene flow (N_m) was calculated to ascertain the conditions of gene communication among populations and was estimated as follows: N_m = (1 − F_ST)/4 F_ST (Rogers and Harpending 1992). F_ST and N_m were calculated using the software ARLEQUIN.

Historical demography and divergence time estimate

The dynamics of population size fluctuations and divergence time were estimated using the Bayesian skyline plot (BSP) method implemented in BEAST 1.7.4 (Drummond and Rambaut 2007). This approach incorporated uncertainty in the genealogy by using Markov chain Monte Carlo (MCMC) integration under a coalescent model, in which the timing of dates provided information about effective population sizes through time. Chains were run for 100 million lengths and the first 10% discarded as ‘burn-in’. The substitution model was selected according the result of jMODELTEST 2.1.1 (Darriba et al. 2012). We applied 10 grouped coalescent intervals and constant growth for the skyline model. Because no fossil data were available to calibrate the mutation rate, we assumed a conventional molecular clock for the avian mitochondrial DNA gene (2 × 10⁻⁸ per site per million years; Weir and Schluter 2008). Demographic history through time was reconstructed using Tracer 1.5 (Drummond and Rambaut 2007).

We conducted two types of analyses for mtDNA in T. tibetanus and T. himalayensis as implemented in Arlequin to examine possible demographic expansions with each population in the species. First, Tajima’s D (Tajima 1989) and Fu’s F_S (Fu 1997) were calculated to assess possible expansions; large negative values of D and F_S statistics are in accordance with the expansion hypothesis. Significance was calculated using 10,000 replicates. Second, pairwise mismatch distributions (Rogers and Harpending 1992) were calculated to examine the demographic expansions of three sympatric locations of the two snowcock species.

Results

Genetic diversity

Genetic diversity parameters, revealed with microsatellites and mitochondria DNA in the two snowcock species, are presented in Table 1. In addition, the compilation of mitochondria genes from GenBank yielded 30 haplotypes of T. tibetanus (Table 1, accession number: JX136799–JX1368336) and 27 haplotypes of T. himalayensis (Table 1, accession number: GQ343513–GQ343549). T. tibetanus in mtDNA showed considerably less haplotype and nucleotide diversity, and fewer polymorphic sites than T. himalayensis (h: 0.83 vs 0.97; π%: 0.43 vs 1.23; S: 34 vs 48). Microsatellites were polymorphic in all populations of T. tibetanus and T. himalayensis. Both average allele number and expected heterozygosity of T. tibetanus were lower than T. himalayensis (Table 1).

Phylogeographical reconstruction

In T. himalayensis, a significant phylogeographical structure and high level of genetic differentiation across the whole range were supported by mtDNA (F_ST = 0.66, p < 0.001) and microsatellites (F_ST = 0.11, p < 0.001). By comparison, the level of genetic differentiation across the whole range of T. tibetanus was lower (mtDNA: F_ST = 0.41, p < 0.001; microsatellite: F_ST = 0.09, p < 0.001). The values of F_ST in T. tibetanus between KLS group and the other groups are much higher than the other pairwise F_ST values based on microsatellite and mitochondria data (Tables 2 and 3). However, the values of F_ST in T. himalayensis between HQLS group and the other groups are much higher than the other pairwise F_ST values based on microsatellite and mitochondria data (Tables 2 and 3).

Table 2 Group pairwise F_ST estimates based on mitochondria DNA in Tetraogallus

Table 3 Group pairwise F_ST estimates based on 7 microsatellite data in Tetraogallus

Despite the broad geographic sampling (Fig. 1), there was little phylogenetical structure in T. tibetanus (Fig. 2). The haplotype network in T. tibetanus revealed no major branching events (Fig. 2), although two groups of haplotypes could be differentiated. The first cluster in T. tibetanus was mostly connected with haplotype T4, which was observed in 25% of individuals and distributed in all groups except WKL. All haplotypes from WKL group were unique and clustered together in other groups with low branch. Interestingly, haplotype H9, which was shared by one individual of T. himalayensis from HQLS and two individuals of T. himalayensis from HWKL group, was clustered in T. tibetanus in the second cluster. For T. himalayensis, HWKL group has some haplotypes, which individually shared with HQLS group or HWTS group.

Fig. 2

Phylogenetic relationships among mitochondrial DNA haplotypes. a Mitochondrial DNA Bayesian phylogeny for T. himalayensis and T. tibetanus. The two species structured clades are separated. Each haplotype of T. himalayensis was named Capital H before numbers and T. tibetanus was named capital T before numbers. The color and the source population of haplotypes are as follows: light blue from QLS, dark blue from TGL, dark pink from BKL, dark blue from QDM, green from KLS, brown from WTS, light pink from ETS, dark green from PME. Abbreviations see details in Fig. 1. Numbers above branches indicate BI posterior probabilities. Only bipartitions with bootstrap or posterior probability values above 50 and 0.5, respectively, are shown. b Median-joining network among mitochondrial DNA haplotypes. The haplotypes are indicated by circles, the size of each circle being proportional to the observed frequency of each haplotype. Lines drawn between haplotypes represent mutation events and small red circles represent missing alleles that were not observed

In the two snowcock species, the FCA plotting of individual microsatellite genotypes (Fig. 3) and the results of STRUCTURE (Fig. 4) detected similar patterns as mtDNA without clean geographic clusters. Although the best-fit K for T. tibetanus was five, and for T. himalayensis it was three, we compared the results from K = 2 to 6. For T. tibetanus, at K = 2 to 5, with WKL groups remained as a single cluster. At K higher than three, the QDM split out, while the QLS fell into the rested clusters with no clear relation to geographic locations (Fig. 4). For T. himalayensis, the bar plots showed three-clade structure (HETS, HQDM and HPME) when K was set to 3. At K = 4, HQLS was separated from the rest clusters, while the other clades showed no significant structure (Fig. 4).

Fig. 3

Principal co-ordinates analyses (PCoA) and genetic structure of the two Tetraogallus species. Factorial correspondence analysis (Benzecri 1973), computed using GENETIX 4.02, shows relationships among the multilocus genotypes of different T. himalayensis and T. tibetanus populations. Axis 1, Axis 2 and Axis 3 are the first, second and third principal factors of variability

Fig. 4

Bar plot derived from Bayesian-based cluster analysis in STRUCTURE. Posterior probabilities (q) for T. himalayensis and T. tibetanus groups are analyzed by STRUCTURE. The proportion of color in each bar represents the assignment probability of an individual, corresponding to different groups of T. tibetanus and T. himalayensis

Historical demography and divergence time estimate

The effective population sizes and demographic trends estimated by the Bayesian skyline plot (BSP) analysis indicated recent population size increases in T. tibetanus (Fig. 5a). Bayesian skyline plots indicated that T. himalayensis maintained a relatively stable population size from 0.06 to 0.003 million years (Mys) ago, after which its population size decreased rapidly (Fig. 5b). The total effective population size of T. tibetanus was always larger than T. himalayensis (Fig. 5). The divergence time between T. tibetanus and T. himalayensis was estimated to be 2.0 Mys ago (95% HPD: 0.9‒3.7 Mys ago).

Fig. 5

Bayesian skyline plot representing historical demographic trends in sampled Tetraogallus tibetanus (a) and T. himalayensis (b). Time is reported on the X-axis as years before present. Estimations were based on a mutation rate of 2 × 10⁻⁸ substitutions per site per year. Ne, the product of the effective female population size and the generation time (in years, log-transformed), is reported on the Y-axis. Estimates are joined by a solid line, whereas dashed lines mark the 95% highest probability density limits

The pairwise comparisons of two snowcocks with their respective outgroup species found no significant deviation from neutrality (p > 0.05) in the McDonald–Kreitman neutrality test (McDonald and Kreitman 1991). The significantly negative values of Tajima’s D and Fu’s F_S support a recent demographic expansion of T. tibetanus (Additional file 1: Figure S1) but stable population size of T. himalayensis with insignificantly positive values of Tajima’s D and Fu’s F_S (see Table 1 and Additional file 1: Figure S1).

Discussion

Genetic variation and the influence of species-specific biological traits

In this study, we investigated phylogeographical patterns for two sister avian species of the Qinghai–Tibetan Plateau and adjacent mountains, which showed a distinct difference in genetic diversity, despite some similar life history traits. Thirty haplotypes with 80 individuals were found in mitochondria DNA regions in T. tibetanus whereas 27 haplotypes with only 67 individuals were detected in the same mitochondria DNA regions in T. himalayensis. The lack of consistency in genetic patterns of sympatric species could reflect historical differences between the species possibly attached historical population isolation and sizes of refugial populations, response to biogeographical barriers, rates of molecular evolution or selective pressures. In 3 of the 25 sampled areas where the two species co-occur, T. himalayensis is more abundant in the western QTP and Pamir Plateau than T. tibetanus (Fig. 1), where the dry season is longer and more severe than in the center of QTP. T. himalayensis was documented to be more tolerant to drought than T. tibetanus (Liu 1999). It is probable that the two sister species might have evolved different tolerances to stresses for cold, hypoxia and drought.

We conducted a comparative assessment of phylogeographical structure in the two snowcocks with biological differences in the Qinghai–Tibetan Plateau and adjacent regions. A striking difference was observed in the genetic diversity, phylogeographic and population dynamics of the two species, which highlights the importance of species-special biological traits, such as breeding behavior and habitat preference on genetic variation and distribution patterns (Fang et al. 2017). It supported previous comparisons of the phylogeography in sympatrically distributed birds, which have revealed the influence of species-specific features on the genetic variation of populations (Qu et al. 2009; Zhang et al. 2012).

Effect of the quaternary glaciations on the phylogeographic pattern and historical demography of two species

Phylogeny and network both point to the importance of intrinsic biological features in structuring genetic diversity and population dynamics. T. tibetanus and T. himalayensis haplotypes are clearly distinct except for H9 from T. himalayensis, which presented in T. tibetanus clusters. It might indicate that at least 4.5% (3/67) of wild T. himalayensis individuals were offspring of female T. tibetanus (Fig. 2). In general, pairwise F_ST was relatively high, with few localities sharing haplotypes in T. himalayensis. T. himalayensis exhibits divergent and structured lineages. Species on the western low-elevation region, with its north extending mountains, appear to have experienced colonization via dispersal followed by isolation and divergence, which has many discontinuous mountain ranges. They appear to have experienced fragmentation, sometimes staged, of wide-ranging ancestral populations. The most extensive divergence occurred within Gansu, Qinghai and Xinjiang, separating eastern and western populations (Wang et al. 2011). For T. himalayensis, a vicariant event, as a consequence of geological isolation and desert expansion, might have produced the significant divergence between the groups from Gansu and Xinjiang (Janecka et al. 2017). We also found KLS group had the biggest pairwise values among other groups in T. tibetanus, which is consistent with the fragmented distribution of the Snow Leopard in this area (Li et al. 2016). As T. himalayensis inhabits middle to high elevation, human activities including mining and infrastructure are more harmful to T. himalayensis genetic diversity. Smaller effective population sizes in T. himalayensis than in T. tibetanus throughout history could also explain the differences between these two species.

The species T. tibetanus did not display a clear phylogeographic pattern. The lack of distinct mtDNA lineages is consistent with previous studies in the QTP region showing weaker Pleistocene refugia effects (Lei et al. 2014; Janecka et al. 2017). In the Central Asian mountains, the Quaternary glaciations have been well developed and intensively influenced the birds of this mountain region (Lei et al. 2014). During cold periods, many warm-adapted species were locally extirpated, due to the advance of major ice sheets. With respect to the mountain cold-adapted animals, the glacier during the cold glacial periods has created connectivity and interglacial refugia that might have served as a gene pool and maintained high levels of genetic diversity during the warming periods (An et al. 2015; Li et al. 2016). One of the haplotypes T4 existed in each group but WKL group (Fig. 2). It means the founder effect could be prevailing during the glacial cold in the QTP. Refugia usually have a relatively stable climate and complex landscape topography, and thus offer the best chance for survival of many taxa when climate changes (Ashcroft et al. 2012), which might explain no genetic structure in T. tibetanus. The absence of T. tibetanus genetic structure is also probably due to recent expansion, which can be seen in Bayesian skyline analysis and mismatch distribution of genetic variations. The two species are sympatric in the Kunlun Mountains and Qilian Mountains, although they have different habitat preferences. T. tibetanus prefers habitat with dwarf shrub and may not extend to areas with elevations lower than 3000 m. In contrast, T. himalayensis in the warmer and drier habitats in low-elevation environments appears to have been less affected by the relatively recent, milder glaciations, and more so by harsher and extensive glaciations (Liu 1999). Hence, the T. himalayensis is more influenced by geographical barriers to gene flow than T. tibetanus in higher altitude (more than 3000 m) and human activities, which are likely to decrease habitat connectivity and cause population decline. Notably, T. himalayensis populations experienced decline since 0.003 Mys ago due to habitat loss, fragmentation, and increasing human activities in the relatively low altitude. We consider these are the main factors shaping the genetic structure of T. himalayensis. The ridges and valleys of the mountains create physical barriers that limit animal dispersal and cause deterministic local variation in rainfall during periods of disturbances of its natural habitats. In contrast, with more plastic ecological requirements, better adaptability and less complex breeding behavior, T. tibetanus has better chances to expand its distribution, resulting in less geographically structured genetic patterns. T. tibetanus as a bird which inhabits the alpine zone between snow line and tree line, is less affected by human activities, which is consistent with the sympatric species Snow Leopards (Li et al. 2016). These results indicate that species with different life histories may respond differently when facing environmental fluctuations during ice ages as well as present-day habitat fragmentation (Qu et al. 2012).

With various behavioral and ecological characteristics, T. tibetanus is more of a habitat generalist than T. himalayensis, with less specific requirements on resources and with more plasticity in food selection (Liu 1999). Interaction between vicariance, dispersal and habitat fragmentation produced the current distribution and genetic diversity. Thus, elevation, topography and cold tolerance appear to drive evolutionary patterns of diversification and demography even among closely related taxa. The comparison of multiple species in genealogical analyses can lead to an understanding of the evolutionary drivers (Qu et al. 2012).

Our results showed low levels of species-level nuclear (F_ST = 0.098) and high mitochondria genetic differentiation (F_ST = 0.601), with no haplotype sharing between the two species. Phylogenetic tree and haplotype networks reaffirmed the sympatric population of H9 of physical barriers in structuring neutral genetic diversity. The two species might have experienced different evolutionary history throughout their current distribution. In the present study only a limited number of samples from Western Kunlun Mountain were analyzed, as the region is difficult to access. H9 shared with the three T. himalayensis individuals was clustered in T. tibetanus, but the fact that the unique mitochondria haplotype representative of these two species was not found elsewhere in the populations analyzed suggests at best a limited role of their coexistence. This could be further clarified by more extensive sampling in Western Kunlun Mountains.

In addition, we discuss the main probable factors maintaining species integrity. How do two snowcock sister species T. tibetanus and T. himalayensis coexist? In the last decades, natural breeding in combination with artificial rearing and massive relocations either for harvesting purposes (hunting or fishing) or for sustaining populations of non-game species mostly dealt with vertebrates. This impressive pace of wildlife relocation is raising increasing conservation concern, as the consequential reshuffling may affect taxa in the wild to such an extent that it results into the loss of biological distinctiveness among regions following the replacement of native biotas by locally expanding non-natives.

Conclusion

By comparing genetic structures estimated from mtDNA and nuclear microsatellite DNA, we show that T. himalayensis and T. tibetanus display different spatial genetic structure and demographic history in response to past climate changes. This study supports the hypothesis that the difference in the snowcocks’ habitat preference would drive their different phylogeographical patterns, even after experiencing several similar geological events, which might have contributed to the different genetic patterns.