首页分享多种有蹄类的演化历史和种群动态揭示东喜马拉雅

多种有蹄类的演化历史和种群动态揭示东喜马拉雅

来源：花匠小妙招时间：2024-11-06 09:07

INTRODUCTION

The Quaternary represents one of the most significant periods of climatic change in Earth’s history, characterized by repeated glacial advances and retreats (Barnosky, 2005). These climatic shifts were largely driven by the periodic uplift of the Qinghai-Xizang Plateau, caused by regular changes in the Earth’s orbital parameters (Candy, 2023). The Quaternary is divided into two epochs: the Pleistocene, spanning approximately 2.6 million (Ma) to 11 700 years ago (Barnosky, 2005; Schaefer et al., 2016) and the Holocene, beginning approximately 11 700 years ago and continuing to present day (Gibbard, 2015; Head et al., 2015).

During the late Mesozoic and Quaternary periods, the Qinghai-Xizang Plateau underwent multiple uplift phases and experienced complex climatic transformations. The Pleistocene, in particular, was shaped by these repeated uplift events and the formation of extensive ice sheets, ice caps, and valley glaciers. According to the record, in addition to the Qinghai-Xizang Plateau and its adjacent region experienced tectonic rotation and strike-slip motion with minor uplift around 8 Ma (Zhang et al., 2006), it also underwent five major uplift phases, including the A phase (3.6 Ma), B phase (2.6 Ma), and C phase (1.7–1.66 Ma) of the Qinghai-Xizang Movement, followed by the Kunlun-Huanghe Movement (1.1–0.7 Ma) and Gonghe Movement (0.15 Ma) (Cui et al., 2011; Li et al., 2015). During the A phase of the Qinghai-Xizang Movement (3.6 Ma), the entire Qinghai-Xizang Plateau was elevated, leading to the disintegration of the primary planation plane (Raymo & Ruddiman, 1992). This uplift coincided with a critical period of global cooling, resulting in widespread glacier formation, although Asia remained largely unaffected due to insufficient mountain elevation for glacier development (Cui et al., 2011). During the B phase of the Qinghai-Xizang Movement (2.6 Ma), the plateau further uplifted to a critical height of 2 000 m, triggering the development of the Asian monsoon system (Li et al., 2001), aridification of the Asian interior, and cooling of the northern hemisphere (Cai et al., 2012; Zachos et al., 2001). During the C phase of the Qinghai-Xizang Movement, several Himalayan peaks surpassed 4 000 m in height, combined with the climate of the ice age at that time, the earliest Shishapangma ice age occurred (Zheng & Shi, 1976). The Kunlun-Huanghe Movement (1.1–0.7 Ma) further elevated the plateau beyond 3 000 m, creating conditions favorable for widespread glaciation. Glaciers expanded southward to the Himalayas, northward to the Kunlun and Tianshan mountains, eastward to the Hengduan Mountains, and westward to the Karakoram Mountains, marking the onset of the glaciation period across the region (Cui et al., 1998). Geological changes, such as mountain uplift, fragmented habitats and created barriers to gene flow, accelerating the genetic differentiation of local species (Antonelli et al., 2009; Xu et al., 2010). Fossil and moraine records indicate that the Qinghai-Xizang Plateau and surrounding regions experienced multiple Pleistocene glaciations during Pleistocene periods (Shi et al., 1992; Shi, 2002), including the Xixiabangma (1.17–0.8 Ma), Naynayxungla (0.7–0.5 Ma), Guxiang (0.3–0.13 Ma), and Baiyu (0.07–0.01 Ma) glaciations (Zheng et al., 2002).

In response to these spatiotemporal fluctuations in environmental conditions, refugia—areas that provide stable conditions for species survival—played critical roles in species persistence (Hewitt, 2000; Schmitt & Varga, 2012; Stewart et al., 2010). These refugia offered survival geographical opportunities for retreated species during adverse climatic periods, then they recolonizated during more favorable conditions (Morales-Barbero et al., 2018). Mountains, in particular, served as important refugia due to their high habitat heterogeneity, complex topography, climatic diversity, and low human disturbance (Zu & Wang, 2022). Fossil evidence from mainland China suggests that early Pliocene fauna adapted to warm, wet conditions retreated southeastward during glacial periods when northern regions became much cooler and drier (Li et al., 2004). Temporary mountain glaciers acted as physical barriers, isolating populations in refugia and promoting population differentiation and intraspecific diversification (Hewitt, 2000; Meng et al., 2015). Consequently, refugia tend to exhibit high species diversity (Gómez & Lunt, 2007; Whittington-Jones et al., 2008; Willner et al., 2009) and a high proportion of endemic species (Fjeldsaå & Lovett, 1997).

The Holocene, in contrast, has been marked by more variable climatic changes, influenced by both natural and human factors. Climatic fluctuations during the Quaternary are recognized as key drivers of evolutionary succession, shaping genetic and geographic structures, and influencing demographic dynamics of across various regions (Domingo et al., 2013; Lorenzen et al., 2011; Mapelli et al., 2012). Indeed, many species experienced population contractions and expansions in response to cold-warm climatic oscillations throughout the Quaternary, which significantly impacted their diversity and distribution (Fu & Wen, 2023). Climate change remains a major determinant of vegetation distribution and abundance (Bartlein et al., 1986; Ma et al., 2013), which, in turn, governs the population dynamics of herbivores by determining food availability and energy resources.

The East Himalaya-Hengduan Mountains (EHHM) region, a key extension of the Qinghai-Xizang Plateau, is one of the most well-known biodiversity hotspots in the world (Myers et al., 2000). This region not only serves as a center for the diversification and origin of many taxa (Hu, 1994; Jia et al., 2012) but is also considered as a critical glacial refuge for many species during the Quaternary glaciations (Frenzel et al., 2003; Wu, 1988; Yang et al., 2008). The EHHM harbors a rich diversity of mammals, including a notable number of endemic and rare species, such as Muntiacus vaginalis, Muntiacus gongshanensis, Muntiacus putaoensis, Naemorhedus cranbrooki, Naemorhedus evansi, Naemorhedus baileyi, Capricornis rubidus, Capricornis sumatraensis (Li et al., 2017; Lwin et al., 2021; Zhang et al., 2021b), Macaca munzala (Sinha et al., 2005), and Rhinopithecus strykeri (Geissmann et al., 2011; Yang et al., 2022). Given the complex topography of the region, the EHHM likely contained many refugia during glacial cycles, providing critical habitats for species survival and fostering biodiversity through isolation and subsequent diversification.

Mountain ungulates provide an ideal model for investigating the impacts of geological and climate changes on population history and dynamics due to their adaptation to alpine environments, where they endure harsh and highly variable conditions (Willisch et al., 2013). This is particularly relevant in the EHHM region, which has a complex geological and climatic history. Among mountain ungulates, serows are of significant interest, with four species listed in the IUCN (International Union for Conservation of Nature) Red List, including Japanese serow (Capricornis crispus), Formosan serow (Capricornis swinhoei), Burmese red serow (Capricornis rubidus), and mainland serow (Capricornis sumatraensis). The classification of Capricornis sumatraensis remains controversial, with ongoing debate about whether the Chinese serow (Capricornis milneedwardsii), Himalayanserow (Capricornis thar) and Sumatran serow (Capricornis sumatraensis) should be considered a species or subspecies as mainland serow (Capricornis sumatraensis) or three distinct species (Supplementary table S1 for details). Gorals, once classified within the same genus as serows, were reclassified into a separate genus in 1985 (Groves & Grubb, 1985). The IUCN Red List contains three goral species, including the long-tailed goral (Naemorhedus caudatus), Himalayan goral (Naemorhedus goral), and red goral (Naemorhedus baileyi). However, whether N. baileyi and N. cranbrooki represent the same or different species remains an unresolved taxonomic issue (Supplementary table S2 for details). Thirteen species of muntjacs are included in the IUCN Red List, with the main taxonomic dispute being whether Muntiacus gongshanensis and Muntiacus crinifrons are the same or different species (Supplementary table S3). Although Capricornis, Naemorhedus, and Muntiacus species share a common evolutionary history and similar life-history strategies, they differ in dispersal abilities and habitat preferences, which generally makes them more vulnerable to extinction compared to smaller mammals (Liow et al., 2008). At least three species or subspecies of each genus are found in the EHHM region. Among serows, C. rubidus, C. thar, C. sumatraensis, and C. milneedwardsii have been reported in the EHHM region (Dou et al., 2016; Lwin et al., 2021), and recent phylogenetic studies suggest the latter three should be classified as a single species—mainland serow (C. sumatraensis) (Mori et al., 2019). Similarly, four goral species have been reported in the EHHM region, including the Burmese goral (N. evansi) (Li et al., 2020), Cranbrook’s goral (N. cranbrooki) (Hayman, 1961; Li et al., 2020), N. baileyi (Velho et al., 2012), and N. goral (Liu & Jiang, 2017). Four muntjac species have been reported in the EHHM, including M. putaoensis (Li et al., 2017; Lwin et al., 2021), M. gongshanensis (Timmins & Duckworth, 2016; Zhang et al., 2021b), M. vaginalis (Lwin et al., 2021), and M. feae (Zhang et al., 1984).

In this study, we obtained complete mitochondrial genome (mitogenome) sequences from five C. sumatraensis from northern Myanmar, 10 N. cranbrooki (four from Motuo in Xizang, six from northern Myanmar), four N. evansi (two from Gongshan and Deqin in Yunnan, two from northern Myanmar), and four M. gongshanensis (three from northern Myanmar and one from Yunnan). We also obtained a complete mitogenome of N. baileyi from a wild individual for the first time, as previous samples were exclusively from zoo specimens. Using these mitogenomes, we examined the drivers of intrageneric and intraspecific differentiation in Capricornis, Naemorhedus, and Muntiacus species, with a particular focus on identifying potential evolutionary refugia in the EHHM region. Additionally, we investigated the historical demographic responses of four species (C. sumatraensis, N. cranbrooki, N. evansi, and M. gongshanensis) to Quaternary climate fluctuations in the EHHM region. This study provides critical insights into the relationships among Quaternary geological events, glacial cycles, climate change, species differentiation, and refugia formation, while also informing the development of effective conservation strategies for preserving species diversity.

MATERIALS AND METHODS

Sample sources

Five C. sumatraensis specimens, ten N. cranbrooki specimens, four N. evansi specimens, four M. gongshanensis specimens, and one N. baileyi specimen were obtained from the EHHM region (Figure 1), including three legs, 11 skulls, eight skins, and two dried meats (Supplementary table S4 for details). Samples were collected from ungulates that died naturally during multiple field expeditions and from specimens confiscated by forest police from hunters. All samples were initially identified based on morphological characters recorded in Groves & Grubb (2011). Only samples with a known locality were collected. All samples were dried and preserved in silica gel and deposited in the Southeast Asia Biodiversity Research Institute, Chinese Academy of Sciences, Myanmar Lab and the School of Life Sciences, Qinghai Normal University.

Figure 1. Geographic information on Capricornis sumatraensis, Naemorhedus cranbrooki, Naemorhedus evansi, and Muntiacus gongshanensis samples with known sampling locations

Samples included those obtained and those downloaded from the NCBI, with different species marked with different colored dots. Map in the upper right corner marks all species samples with known geographical locations. Red box marks East Himalaya-Hengduan Mountains (EHHM) region. MT and CY represent Motuo and Chayu in Xizang, China; GS and DQ represent Gongshan and Deqin in Yunnan, China; BX represents Baoxing in Sichuan, China; NM represents northern Myanmar.

DNA extraction, polymerase chain reaction (PCR) amplification, and sequencing

All molecular work was performed in a separate facility under in a sterile environment to avoid cross-contamination. DNA was extracted using a DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany), and its concentration and integrity were detected using 2% agarose gel electrophoresis. The volume of the master mixture for PCR was 50 µL, which contained 100 ng of template DNA, 1 μL (10 pmol) of each primer, 5.0 µL of 10×Taq buffer, 2 μL of dNTP (2.5 mmol/L), 2.0 U Taq DNA polymerase, and sterilized ultra-purified water. The PCR conditions were as follows: predenaturation at 95℃ for 8 min, 35 cycles of denaturation at 95℃ for 45 s, annealing at 52℃–64℃ for 45 s, extension at 72℃ for 1 min, and final extension for 7 min at 72℃. Reactions were performed on a Veriti Thermal Cycler (Applied Biosystems, USA) and always included a negative control. An ABI PRISM 3700 Genetic Analyzer (Applied Biosystems, USA) was used to sequence both ends of the PCR products. The primers refer to Supplementary table S4 of Li et al. (2020).

Assembly and alignment

Chromas software (http://www.technelysium.com.au) was used to analyze the obtained sequences, and waveform diagrams were used to correct the sequences. DNA Star v.5.0 was used to splice genes after correction (Burland, 2000). The sorted sequences were aligned in GENEIOUS v.2019.1.1 using the default parameters of MAFFT v.7.245 (Rozewicki et al., 2019), and manually checked and corrected.

Phylogenetic relationships based on mitogenomes

In addition to the obtained sequences, we downloaded all available complete mitogenomes from the NCBI database (https://www.ncbi.nlm.nih.gov). In total, 20 serows and 26 gorals were used as the ingroup and one Rupicapra rupicapra, one Oreamnos americanus, and one Pantholops hodgsonii (Bibi, 2013; Hassanin et al., 2012) were used as the outgroup to analyze the phylogenetic relationships of Capricornis and Naemorhedus species (Supplementary table S5 for details). Including the obtained and downloaded sequences, 20 complete Muntiacus mitogenomes were used as the ingroup and one Elaphodus cephalophus, one Cervus elaphus, and one Hydropotes inermis (Beller et al., 2006; Pang et al., 2008) were used as the outgroup (Supplementary table S5 for details). All sequences were subjected to pairwise alignment in ClustalW (Higgins & Sharp, 1988). Maximum-likelihood (ML) and Bayesian inference (BI) phylogenetic trees were constructed using IQ-tree and MrBayes, respectively, in PhyloSuite v.1.2.2 (Zhang et al., 2020). The most appropriate nucleotide substitution models were selected using Akaike information criteria (AIC) in ModelFinder (Kalyaanamoorthy et al., 2017). The optimal models for serows and gorals were TIM2+F+R8 (IQ-tree) and GTR+F+I+G4 (MrBayes) and the optimal models for muntjacs were TIM2+F+R3 (IQ-tree) and GTR+F+I+G4 (MrBayes). IQ-tree was run with 100 000 ultrafast bootstrap replicates. For Bayesian analysis, four Markov chain Monte Carlo (MCMC) chains were simultaneously run for 20 000 000 generations, with one tree randomly selected every 1 000 generations. Convergence diagnostics and discard values were assessed using Tracer v.1.7 (Rambaut et al., 2018), and the first 25% of samples were discarded as a conservative burn-in based on the diagnostic results. All trees were visualized using Figtree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree).

Divergence times

Bayesian analysis was conducted using the same dataset as the phylogenetic relationships, employing MCMC implemented in BEAST v.2.7.5 (Bouckaert et al., 2014) with a Yule prior and strict molecular clock. The nucleotide substitution models were all GTR+F+G4, with empirical base frequencies, determined based on AIC in ModelFinder within PhyloSuite v.1.2.2. Three calibration points were used to estimate the divergence time of Capricornis and Naemorhedus species. Two of these calibration points, with age ranges and prior probability distributions, were derived from Bibi (2013): internal calibration of the branching point between Caprini and Pantholops (95% highest posterior density interval (HPD)=9.3–11 Ma), and crown Caprini with a normal prior (mean=8.9 Ma, standard deviation=2 Ma) based on Aragoral mudejar (see Additional File 1 in Bibi, 2013). Another calibration point was based on the Sichuan (China) fossil record of Naemorhedus goral from the Mid-Pleistocene (0.126 –0.781 Ma; Mead, 1989). Three calibration points were used to estimate the divergence time of Muntiacus species: the crown Cervidae (12.4–9.0 Ma) given in Bibi (2013), the split between Cervini and Muntiacini at 9±1 Ma based on fossils from the Late Miocene of China (Dong, 2007; Petronio et al., 2007), and the split between the red and black muntjacs (M. crinifrons), estimated at around 3.29 Ma (95% HPD=3.36–5.88; Singh et al., 2021). Three independent MCMC runs were performed, with sampling every 1 000 generations for 100 000 000 generations. A burn-in of 20% was applied, and the convergence of all parameters was assessed using Tracer v.1.7. The maximum clade credibility criterion tree was obtained using TreeAnnotator and visualized with Figtree v.1.4.4.

Demographic history

The Extended Bayesian Skyline Plot (EBSP) was used to infer the historical population dynamics (Miller et al., 2018). EBSP analyses were implemented using the BEAST v.2.7.5 package. Complete mitogenomes of four species from the EHHM region (C. sumatraensis, N. cranbrooki, N. evansi, and M. gongshanensis) were selected for EBSP analysis with a strict molecular clock. The nucleotide substitution models were determined based on Bayesian Information Criterion (BIC) in ModelFinder within PhyloSuite v.1.2.2. The optimal model was HKY+F for all four species. The “Coalescent Extended Bayesian Skyline” model was used, with the “Population Model” factor set to 0.5 to account for female contribution to Ne (Heled & Drummond, 2008). Each population was run for 10 000 000 generations, with sampling every 1 000 steps and the first 10% of trees discarded as burn-in. Each dataset was run twice to confirm repeatability. Results were analyzed using Tracer v.1.7, ensuring that the effective sample sizes (ESS) for all relevant parameters exceeded 200. Demographic reconstructions were plotted in R v.3.2.3 using the R script plot EBSP.R (Heled & Drummond, 2008).

RESULTS

Phylogenetic relationships

The ML and BI analyses, based on complete mitogenomes, yielded similar tree topologies. Bootstrap values indicated significant support for all nodes. Naemorhedus and Capricornis were monophyletic sister genera (Figure 2). Three major phylogenetic clades were identified within Capricornis: 1) C. crispus; 2) C. swinhoei and C. rubidus; and 3) C. milneedwardsii, C. thar, and C. sumatraensis. The five Capricornis sequences collected in this study clustered within the third branch, along with C. milneedwardsii, C. thar, and C. sumatraensis. Two major phylogenetic clades were identified within the Naemorhedus: 1) N. cranbrooki, N. baileyi, and N. evansi—with N. cranbrooki and N. baileyi forming monophyletic taxa; and 2) N. griseus, N. goral, and N. caudatus. In addition, N. caudatus was divided into two distinct subclades, with our N. cranbrooki sequences from Xizang clustering in one subclade and those from northern Myanmar clustering in another. Similarly, N. evansi was divided into two different subclades, with our sequences from Yunnan and northern Myanmar grouped together, and N. baileyi was divided into two different subclades. The Muntiacus species formed two major clades (Figure 2): 1) M. reevesi, M. vuquangensis, and M. putaoensis; and 2) M. feae, M. gongshanensis, M. crinifrons, red muntjac (M. vaginalis, M. muntjak, and M. malabaricus), and M. atherodes. Both M. crinifrons and M. gongshanensis were identified as monophyletic, with M. gongshanensis further divided into two subclades. The sequences of M. gongshanensis from northern Myanmar were dispersed across both subclades.

Figure 2. Unified maximum-likelihood (ML) and Bayesian inference (BI) tree and estimation of divergence time based on complete mitogenomes

A: Unified tree of Capricornis and Naemorhedus species. B: Unified tree of Muntiacus species. Black numbers at nodes indicate maximum-likelihood and Bayesian posterior probability, respectively. Stars indicate values of 100 (ML) and 1.00 (BI). Blue numbers at nodes represent average time (Ma), blue bar represents 95% confidence interval. Sequences in bold are our samples. ☆ indicates sample is an ancient individual (sub-fossils).

Divergence times

The evolutionary tree topologies for Capricornis, Naemorhedus, and Muntiacus species were consistent with the phylogenetic tree topologies (Figure 2), thus the two trees were combined into a single unified tree.

The split between Muntiacus species and Elaphodus cephalophus occurred at 7.82 Ma (95% HPD=7.27–8.39). The split between Capricornis and Naemorhedus species occurred at 5.19 Ma (95% HPD=4.66–5.74), with two major branches of Naemorhedus separating 4.83 Ma (95% HPD=4.33–5.34). Muntiacus atherodes split from the common ancestor of red muntjacs+M. feae+M. gongshanensis+M. crinifrons at 3.11 Ma (95% HPD=2.85–3.38), while M. reevesi split from the common ancestor of M. vuquangensis+M. putaoensis at 3.04 Ma (95% HPD=2.76–3.11) and M. vuquangensis split from M. putaoensis at 2.26 Ma (95% HPD=2.03–2.5). The red muntjacs split from the common ancestor of M. feae+M. gongshanensis+M. crinifrons at 2.69 Ma (95% HPD=2.46–2.92). Capricornis crispus diverged from the common ancestor of Capricornis species at 2.64 Ma (95% HPD=2.33–2.95). The split between N. evansi and the common ancestor of N. cranbrooki and N. baileyi occurred at 2.52 Ma (95% HPD=2.23–2.82). The common ancestor of C. swinhoei and C. rubidus diverged at 1.51 Ma (95% HPD=1.31–1.7), while the split between N. cranbrooki and N. baileyi occurred at 1.5 Ma (95% HPD=1.3–1.68). Muntiacus feae split with the common ancestor of M. gongshanensis+M. crinifrons at 1.24 Ma (95% HPD=1.1–1.37), while the split between M. vaginalis and M. muntjak occurred at 0.91 Ma (95% HPD=0.8–1.02).

Regarding the EHHM region, the split between M. gongshanensis and M. crinifrons occurred at 0.74 Ma (95% HPD=0.64–0.84), while the two subclades of M. gongshanensis split at 0.18 Ma (95% HPD=0.14–0.22). The mainland serow (Capricornis sumatraensis) predominantly divided into two main subclades 0.42 Ma (95% HPD=0.36–0.48). The two subclades of N. cranbrooki split at 0.42 Ma (95% HPD=0.36–0.5), while the two subclades of N. evansi split at 0.35 Ma (95% HPD=0.29–0.42). A series of rapid divergences occurred within mainland serow between 0.38 Ma and 0.14 Ma.

Population dynamics

The EBSP analysis of C. sumatraensis, N. cranbrooki, N. evansi, and M. gongshanensis from the EHHM region is presented in Figure 3. The population of C. sumatraensis remained stable before 1.5 Ka (kilo annum), after which it experienced an expansion. Similarly, the population of N. cranbrooki remained stable until approximately 1 Ma, followed by expansion. The populations of both N. evansi and M. gongshanensis were stable until approximately 8 Ka, with population expansion occurring after that time.

Figure 3. Extended Bayesian skyline plot (EBSP) analysis of Capricornis sumatraensis, Naemorhedus cranbrooki, Naemorhedus evansi, and Muntiacus gongshanensis in the EHHM based on complete mitogenomes

A: EBSP analysis of C. sumatraensis. B: EBSP analysis of N. cranbrooki. C: EBSP analysis of N. evansi. D: EBSP analysis of M. gongshanensis. X axis represents time in millions of years (Ma). Y axis represents effective population size. Black dashed line represents median of population size, and green shaded area represents 95% confidence interval.

DISCUSSION

Phylogenetic relationships, species differentiation, and refugia

Geological events, particularly the uplift of the Qinghai-Xizang Plateau, are considered key drivers of species divergence and population differentiation (Fu & Wen, 2023). The uplift of mountains creates physical barriers that promote species differentiation. The differentiation within Capricornis, Naemorhedus, and Muntiacus species is closely associated with multiple phases of uplift of the Qinghai-Xizang Plateau. The divergence between Muntiacus species and Elaphodus cephalophus, which occurred 7.82 Ma, may have been triggered by a slight uplift of Qinghai-Xizang Plateau at 8 Ma. The Naemorhedus and Capricornis species split around 5.19 Ma, with two major subclades of Naemorhedus diverging at approximately 4.83 Ma, corresponding to slight uplift of the Qinghai-Xizang Plateau and Hengduan Mountains. The geographic barriers formed by the A phase of the Qinghai-Xizang Movement may have contributed to the divergences of M. atherodes and M. reevesi, which occurred at 3.11 Ma and 3.04 Ma, respectively. The split between the three red muntjac species and the common ancestor of M. feae+M. gongshanensis+M. crinifrons occurred at 2.69 Ma, while the divergence between N. evansi and the common ancestor of N. cranbrooki+N. baileyi occurred around 2.52 Ma. Both of these events may have been triggered by mountain uplift during the B phase of the Qinghai-Xizang Movement. The basal species C. crispus, which inhabits the islands of Japan, diverged from the common ancestor of other serows around 2.64 Ma, coinciding with the appearance of the Tsushima Strait land bridge during the late Pliocene, which connected Japan and mainland China (Kimura, 1996, 1977). After disappearance of the land bridge at 3 Ma, species may be isolatied (Kimura, 2002; Ogasawara, 1994). Additional divergences occurred between C. sumatraensis and the common ancestor of C. swinhoei+C. rubidus approximately around 1.51 Ma, between N. cranbrooki and N. baileyi around 1.5 Ma, and between M. feae and the common ancestor of M. gongshanensis+M. crinifrons approximately at 1.24 Ma. These divergences were likely driven by habitat fragmentation resulting from the C phase of the Qinghai-Xizang Movement.

Pleistocene glaciers may have isolated populations into refugia and promoted population differentiation and intra-specific diversification without apparent geological breaks (Hewitt, 2000). Our divergence analysis supported the interspecific and intraspecific differentiation of Capricornis, Naemorhedus, and Muntiacus species were closely linked to a series of Pleistocene glaciers. The divergence between M. vaginalis and M. muntjak occurred approximately 0.91 Ma, coinciding with the Xixiabangma glaciation (0.8–1.17 Ma; Zheng et al., 2002). This glaciation likely contributed to speciation, as arid regions of central India were unsuitable for many species, forcing populations to seek refuge in more favorable forested habitats in the northern and southern regions of the Indian subcontinent (Martins et al., 2017). The divergence between M. gongshanensis and M. crinifrons, two independent species, occurred around 0.74 Ma, likely due to the significant mountain uplift during the Kunlun-Huanghe Movement (1.1–0.7 Ma). The development of glaciers possibly forced these species to retreat to different refugia within the EHHM region. Subsequently, the distribution range of M. crinifrons likely expanded during the interglacial period. Consistent with the findings of Mori et al. (2019), our results showed that C. thar, C. sumatraensis, and C. milneedwardsii clustered into a single clade, but were not monophyletic, suggesting they may represent a single species (C. sumatraensis). The Naynayxungla glaciation, the most extensive glaciation event impacting the EHHM region (Wu et al., 2022; Zheng et al., 2002), likely contributed to the intraspecific differentiation of C. sumatraensis and N. cranbrooki around 0.42 Ma, and of N. evansi around 0.35 Ma (95% HPD=0.29–0.42). In the case of C. sumatraensis, the population was divided into two subclades, with samples from northern Myanmar dispersed across both clades. Similarly, N. cranbrooki samples from northern Myanmar and Xizang were divided into two subclades. We speculate that glacial expansion forced these populations into different refugia within the EHHM region, where they survived. A series of rapid intraspecific differentiations occurred in C. sumatraensis between 0.39 Ma and 0.07 Ma, potentially driven by the glacial and interglacial cycles, suggesting repeated dispersal occured centering from different refugia during the Quaternary period. According to Song et al. (2023), the regional extirpation of serows in northern China may have occurred without significant impact on genetic diversity, as individuals from Beijing and Baoxing (Sichuan) are closely related. A possible regional replacement in Guizhou, driven by maternal lineage, is indicated by sub-fossils from Guizhou and all modern serows from the same area clustering into different subclades. The intraspecific differentiation of N. cranbrooki was likely triggered by the Naynayxungla glaciation, which forced populations to retreat into refugia in Xizang and northern Myanmar. Naemorhedus evansi was divided into two subclades: one consisting of samples from Thailand and Guizhou (China), and the other from Myanmar and Yunnan (China). This divergence potentially occurred on both sides of the Nu-Salween River, which formed in the middle Neogene period (4.2–5.0 Ma) (Zhao et al., 2011). Glaciers and seasonal snowpacks, essential for regional water cycles, influence the flow of downstream water resources (Hugonnet et al., 2021; Zemp et al., 2019). Thus, we speculate that the Nu-Salween River may have experienced several low flow periods under the influence of the Xixia Bangma Glacier, similar to those of the Yarlung Zangbo River (Wu et al., 2022), potentially enabling N. evansi to cross the river at shallow points. In addition, the valley glaciers formed during the glacial period may also have served as natural bridges for N. evansi dispersal. The intraspecific differentiation of M. gongshanensis around 0.18 Ma also suggests the presence of different refugia within the EHHM region during the Guxiang glacial period.

Species population expansion

The EBSP analyses revealed that the population expansions of N. cranbrooki, N. evansi, and M. gongshanensis in the EHHM region began between 10 and 8 Ka, coinciding with the early to mid-Holocene period (11–7 Ka). During this time, the Northern Hemisphere experienced rapid warming and increased humidity, primarily due to the retreat of continental ice sheets and orbital-induced changes in seasonal insolation (Zhang et al., 2022). These climatic shifts during the early Holocene promoted vegetation growth, particularly temperate deciduous broadleaf forests. In the middle Holocene (7–3 Ka), although the climate became slightly drier, it still provided a relatively warm environment (Zhang et al., 2021a). Global warming likely had a pronounced effect on vegetation phenology in high-altitude mountain regions, extending the vegetative period (Kullman, 2004) and increasing food availability and energy for herbivores, thus facilitating population growth. Consequently, for herbivorous species inhabiting extreme mountain environments, where food availability is limited by cold temperatures during the growing season (Ciach & Pęksa, 2018), this climatic improvement likely had a positive effect on their survival and reproduction. Winter is a particularly critical period for offspring survival (Willisch et al., 2013), and the warmer winter climate of the Holocene, compared to the Pleistocene, would have been conducive to survival. As herbivorous ungulates adapted to high mountain regions, the population expansions of N. cranbrooki, N. evansi, and M. gongshanensis were highly consistent with the warmer Holocene climate. EBSP analysis indicated that the population expansion of C. sumatraensis in the EHHM region began approximately 1.5 Ka, during the late Holocene. At that time, environmental conditions on the Qinghai-Xizang Plateau deteriorated, with lower temperatures and reduced precipitation leading to forest retreat toward the southeastern margin (Zheng & Li, 1999). This environmental shift may have driven the serow to migrate toward the southeastern regions, contributing to its population expansion in the EHHM. Furthermore, after 1.6 Ka, intensified human settlement and grazing activities on the northeastern Qinghai-Xizang Plateau (Huang et al., 2017) may have further displaced C. sumatraensis, promoting its population expansion in the EHHM region.

In this study, we report on the collection of four samples of N. cranbrooki obtained in Motuo (Xizang) and two samples of N. evansi obtained in Yunnan (one from Gongshan and one from Deqin) for the first time, extending the distribution range of both species and suggesting a previous underestimation of their ranges. Naemorhedus cranbrooki was traditionally thought to be restricted to northern Myanmar (Li et al., 2020); however, our discovery confirms its presence in Xizang, China. Similarly, while N. evansi was previously thought to be widespread from Thailand to central China (Li et al., 2020), no molecular evidence had confirmed its presence in Yunnan until our samples. It's worth noting that the sequence MH155202 of serow ought to be C. milneedwardsii which belongs to a subspecies of C. sumatraensis in the original documentary record (Zhao et al., 2019), but it is listed as C. sumatraensis in NCBI. Accurate species identification is crucial for subsequent studies, and we recommend that this sequence be correctly classified as C. milneedwardsii in NCBI to ensure taxonomic precision.

ACKNOWLEDGMENTS: This work was supported by the Forest Research Institute, Hponkanrazi Wildlife Sanctuary, and Hkakaborazi National Park Offices, Myanmar. Laboratory work was completed in the Central Laboratory, Public Technology Service Center, Xishuangbanna Tropical Botanical Garden (XTBG), CAS and the laboratory in the School of Life Sciences, Qinghai Normal University.