Inter-glacial isolation caused divergence of cold-adapted species: the case of the snow partridge.

Deciphering the role of climatic oscillations in species divergence helps us understand the mechanisms that shape global biodiversity. The cold-adapted species may have expanded their distribution with the development of glaciers during glacial period. With the retreat of glaciers, these species were discontinuously distributed in the high-altitude mountains and isolated by geographical barriers. However, the study that focuses on the speciation process of cold-adapted species is scant. To fill this gap, we combined population genetic data and ecological niche models (ENMs) to explore divergence process of snow partridge (Lerwa lerwa). Lerwa lerwa is a cold-adapted bird that is distributed from 4,000 to 5,500 m. We found 2 genetic populations within L. lerwa, and they diverged from each other at about 0.40-0.44 million years ago (inter-glacial period after Zhongliangan glaciation). The ENMs suggested that L. lerwa expanded to the low elevations of the Himalayas and Hengduan mountains during glacial period, whereas it contracted to the high elevations, southern of Himalayas, and Hengduan mountains during inter-glacial periods. Effective population size trajectory also suggested that L. lerwa expanded its population size during the glacial period. Consistent with our expectation, the results support that inter-glacial isolation contributed to the divergence of cold-adapted L. lerwa on Qinghai-Tibetan Plateau. This study deepens our understanding of how climatic oscillations have driven divergence process of cold-adapted Phasianidae species distributed on mountains.

Unveiling the mechanisms driving species divergence is crucial to comprehend speciation (Coyne and Orr 2004; Edelman et al. 2019), which will contribute to understanding their current distribution (Boyer et al. 2007), adaptation (Lamichhaney et al. 2015), morphology (Chen et al. 2019), and their evolutionary direction (Sackton et al. 2019). The Quaternary, covering last 2.6 million years ago (mya), is a geological epoch characterized by glacial–interglacial cycles (Shi et al. 2006). The regularly repeated climatic changes have been recognized as one of the major forces driving divergence process with the shifts of the distribution ranges (Swart et al. 2009; Buainain et al. 2020; Tomasello et al. 2020). The temperate species that have conserved ecological niche may have a broad distribution during interglacial periods, but are distributed in isolated narrow refuges during glacial periods (Hewitt 2000). Whereas the cold-adapted species having conserved ecological niche may retreat to the arctic region or the higher mountains during the warm interglacial periods, but they can expand their distribution during the glacial periods (Hewitt 2011; Dapporto et al. 2019; Ehl et al. 2021). Facing ongoing climatic change, understanding how past climatic fluctuations influence divergence will benefit us in forecasting the future of the biota.

The glacial isolation or glacial refugia hypothesis holds the opinion that the different populations were fully or partly isolated by the glacier during the glacial period, and the genetic differences gradually accumulated in the isolated populations under genetic drift and/or ecological selection (Willis et al. 2004). However, this hypothesis can only explain the divergence of species that was unable to adapt to the cold environment. The cold-adapted species expanded their distribution during the glacial era and were isolated by geographical barriers during inter-glacial period (da Silva Ribeiro et al. 2020; Lu et al. 2020; Ehl et al. 2021). Although numerous studies have focused on how Quaternary climatic fluctuations drove divergence, few studies have explored the mechanisms of inter-glacial isolation driving this process (da Silva Ribeiro et al. 2020; Lu et al. 2020; Ehl et al. 2021); this happens, at least in part, because of the difficulty in collecting samples, especially animals, at high altitudes. The previous studies have found interglacial isolation drove divergence process of the cold-adapted butterfly (Boloria pales and Bolorianapaea; Ehl et al. 2021), and the cold-adapted Drab-Breasted Bamboo Tyrant (Hemitriccus diops) could expand its distribution during glacial period (da Silva Ribeiro et al. 2020). Hence, we hypothesized that the inter-glacial isolation may drive divergence process of other cold-adapted species.

Climatic oscillations played a major role in driving biodiversity in mountain regions (Hoorn et al. 2018; Rahbek et al. 2019a, 2019b). The alternation of glacial and inter-glacial periods caused frequent contractions and expansions of the species distribution, leading to isolation as more likely to occur in mountain regions. These areas hold half of the biodiversity hotspots with high biodiversity and endemic species, which makes it an excellent region to study the influence of climatic oscillations on species divergence (Rahbek et al. 2019a, 2019b; Perrigo et al. 2020). The Qinghai–Tibetan Plateau (QTP) and its surrounding mountains are global biodiversity hotspots, which played a considerable role in promoting and maintaining biodiversity throughout time (Myers et al. 2000; Lei et al. 2014). The QTP is in the middle and low latitudes (26°N–39°N), having a complex landscape with a crisscross of mountains and canyons. The glacier only covered the high-altitude mountains on the QTP during glacial period (Shi et al. 2006); hence, the species distributed at high altitudes, compared with those at low altitudes, might have different genetic imprints of Quaternary climatic fluctuations. Most previous studies focused on the species distributed in the middle or low altitude of the QTP, and found their divergence caused, in great extent, by glaciation isolation (Qu et al. 2010; Lei et al. 2014; Liu et al. 2016; Jin and Brown 2019; Zhang et al. 2019). However, few studies have explored the divergence process of high-altitude species on QTP. These organisms can adapt to cold environments, which suggests that they may have expanded their distribution during glacial period and contracted their distribution during inter-glacial period. During inter-glacial periods, they may distribute in isolated high-altitude areas, leading to the divergence process happening at that time. The effective population size of the cold-adapted species may increase during glacial periods and decrease during interglacial periods.

Snow Partridge (Lerwa lerwa, Phasianidae, Galliformes) is distributed in the eastern and southern edges of the QTP (Figure 1), inhabiting alpine meadows, tundra, and rocky regions from 4,000 m to 5,500 m above the sea level, living around the glaciers (Zheng 2015). Three subspecies (L. l. major, L. l. lerwa, and L. l. callipygia) were recognized within this species based on morphological data, such as body size and plumage color. The L. l. major mainly distributed in western Sichuan and northwestern Yunnan of China; L. l. lerwa mainly distributed in southern and southeastern Tibet of China, Nepal, northern India, and eastern Pakistan; and L. l. callipygian mainly distributed in southern Gansu and northern Sichuan of China (Zheng et al. 1978). There is no apparent geographical barrier between the distributions of L. l. major and L. l. callipygian, however, there are several canyons, such as Lancangjiang, Jinshaiiang, and Yalongjiang, located in the adjacent region of L. l. major and L. l. lerwa. Because it is difficult to sample L. lerwa or do the field work at high latitudes, there have been limited reports about L. lerwa. Nevertheless, Yao et al. (2017) reported that L. lerwa ranged from 4,400 to 4,700 m above the sea level after breeding; and Wang et al. (2017) found multispecies coalescent ultra-conserved element analyses supported L. lerwa was a clade sister to a clade that included typical pheasants, turkey, grouse, and tragopans. To date, no study has reconstructed the divergence process within L. lerwa. Distribution in high altitudes and potential divergence within the species makes L. lerwa is an ideal organism to explore how climatic oscillations have driven the divergence process of cold-adapted species.

Figure 1

The location of the samples and the distribution of L. lerwa (blue, red, and purple circles indicate the sampling locations of L. l. callipygia, L. l. major, and L. l. lerwa, respectively. The corresponding shadow indicates the distribution of each sub-species, referring to the IUCN red list, field research experience of Nan Wang, and Zheng G-M. 2015. The stars mark the sampling locations of the 2 outgroups. The source of the map was from https://www.resdc.cn/).

In this study, we integrated population genetic data, ecological niche model (ENM), and effective population size trajectory information to explore the historical factors linked to the divergence process of L. lerwa. Specifically, this study focused on 4 questions: (1) How many genetic populations are within L. lerwa? (2) Did the splitting time within L. lerwa locate in inter-glacial period? (3) Did the effective population size of L. lerwa decrease during inter-glacial period? (4) Did L. lerwa expand its distribution during glacial periods? The answers to these questions are helpful to understand the divergence process of cold-adapted species. This study will help us to have deep understanding of the mechanism that underpins divergence or speciation occurring in high altitudes .

Materials and Methods

Sampling

The samples of L. lerwa were collected across its distribution area; and there were 11 (muscle), 15 (muscle), and 8 (7 muscle and 1 blood) samples were collected from L. l. lerwa, L. l. callipygian, and L. l. major, respectively (Figure 1 and Supplementary Table S1 ). The 34 samples covered most distributions of the 3 subspecies. To estimate phylogenetic relationships and divergence time, we collected blood from a Tibetan snowcock Tetraogallus tibetanus, a distantly related taxon, and a buff-throated partridge Tetraophasis szechenyii, a closely related taxon, as outgroups (Kimball et al. 2021). All samples were kept in a −80 °C ultra-low temperature freezer at the Beijing Forestry University.

Gene amplification and sequencing

Total DNA of each sample was extracted using a DNA extraction kit (TianGen Biotech, Beijing, China) according to the manufacturer’s protocol. There were 29 loci amplified using polymerase chain reaction (PCR), which included 3 mitochondrial and 26 unlinked autosomal exonic loci (Supplementary Table S2). The PCR protocol was as follows: reaction volume of 40 μL containing 200 ng of template DNA, 20 μL of premix (including DNA polymerase, deoxy-ribonucleoside triphosphate (dNTP), Mg2+; TaKaRa, Shiga, Japan), and 0.25 μM of each primer. The same PCR conditions were applied for all genes. There was an initial starting step of 94 °C for 5 min followed by 2 cycles of 94 °C for 30 s, 63 °C for 45 s, and 72 °C for 45 s; which was then followed by 20 cycles with lower annealing temperatures ranging from 63 °C to 52 °C with 0.5 °C steps for each cycle; and finally, 10 cycles of 94 °C for 30 s, 52 °C for 45 s, 72 °C for 45 s, and a final extension at 72 °C for 7 min. The PCR products were sequenced from both 5′- and 3′-directions using an Applied Biosystems 3730XL auto-sequencer (Foster City, CA, USA) at the sequencing center of Sangon Biotech Company.

Analyzing sequences, estimating population structure, and re-constructing phylogeny

The sequences were edited in Seqman version 11.2 (DNASTAR, Inc., Madison, WI, USA), including trimmed low-quality bases and manually mobilized the hybrid bases. The sequences of each locus were aligned in MEGA X (Kumar et al. 2018). Molecular clock tests were performed on each gene using a maximum likelihood (ML) ratio test in MEGA X with the following set. The bootstrap (bp) method was selected to estimate the phylogeny and the bp number was set to 1,000. The best-fit nucleotide substitution model was determined using the Akaike Information Criterion in jModelTest version 2.1.7 (Darriba et al. 2012). The unphased sequences were phased using DNASP, in which the haplotype probability threshold was set 0.9. The number of variable sites and parsimony information sites of each locus were calculated using DNASP (Rozas et al. 2017). We set the substitution rate of cytochrome b gene (cyt b) as 1.19 × 10−8 per site per year, which is an average substitution rate for Galliformes, calibrating by the fossils and bio-geographical events (Weir and Schluter 2008). We calculated the overall mean Tamura-Nei distance (d) for each locus, and then the substitution rate of each locus equaled the corresponding d divided by the d of cyt b and then multiple by 1.19 × 10−8. The nucleotide diversity (π and ϴw) of each locus within each subspecies was calculated using DNASP. The absolute differentiation (DXY) between pairs of subspecies on each locus was calculated in DNASP.

We employed Structure version 2.3.4 with admixture model (Hubisz et al. 2009) to analyze variable sites, inferring the population genetic structure of L. lerwa. The running parameters were set as follows: run length of burin period of 200,000, number of Markov Chain Monte Carlo (MCMC) replications after burin was 1,000,000, K value (potential genetic clusters) varied from 1 to 7, and the repeat number was 10. The results obtained from Structure were summarized in the CLUMPAK server (KopelmAn et al. 2015). Deltak was used to determine the optimal K value. The discriminant analysis of principal components (DAPCs) was also used to infer the population genetic structure using the adegenet package in the R software (Jombart 2008).

Phylogenetic reconstruction was performed using the following data sets: (1) nuclear gene sequences; (2) mitochondrial gene sequences; and (3) the combined nuclear and mitochondrial gene sequences. For each dataset, 2 algorithms (ML, and Bayesian inference BI) were used to infer the phylogenetic relationships. Using the RA×ML software (Stamatakis 2014), we conducted phylogenetic trees using the ML algorithm with 1,000 bp and GTR+GAMMA+I model. The BI trees were performed using MrBayes version 3.1.2 (Ronquist and Huelsenbeck 2003). During ran MrBayes, we performed 2 runs with 5,000,000 generations; 5 independent chain samplings every 1,000 generations with the burn-in was 25%.

Estimating divergence time

To explore whether the divergence occurred during glacial period or inter-glacial period, the *BEAST version 2.2.1 (Bouckaert et al. 2014) and Ima2p version 1.0 (Sethuraman and Hey 2016) were used to estimate divergence time between 2 genetic populations within L. lerwa. Both Structure and DAPC suggested there were 2 genetic populations within L. lerwa. The *BEAST hypothesized that there was no gene flow during the divergence process (Bouckaert et al. 2014); however, the Ima2p (Isolation with migration [IM] model) hypothesized that there were immigrants along the speciation (Sethuraman and Hey 2016). The combined dataset that included mitochondrial and nuclear loci was used for both analyses. *BEAST was performed in BEAST version 2.2.1. Because the lineages within L. lerwa have been completed sampling, the phylogenetic relationship was reconstructed under the Yule speciation process with a lognormal prior distribution (Steel and McKenzie 2001). The substitution model, clock model, and mutation rate of each allele are listed in Supplementary Table S2. The number of MCMC steps was set as 1,000,000,000, with every 5,000 steps sampled, and the first 10% of steps were discarded as burn-in. The program Tracer version 1.5 (Rambaut et al. 2018) was used to test the validity of the sampling results (effective sample size [ESS] >200). The final tree was summarized by the program TreeAnnotator version 1.7.4 (Bouckaert et al. 2014), using the common ancestor heights as node heights with the first 20% as burn-in and posterior probability limit was 0.5. Finally, we employed FigTree version 1.3.1 (Bouckaert et al. 2014) to visualize the tree. We performed the Ima2p program with 100 MCMC chains and a geometric increment model. The total number of MCMC steps was set to 10,000,000 with 100,000 steps, and the number of burin steps was 50,000. The substitution rate of each allele is shown in Supplementary Table S2. It is difficult to determine the average generation time of a population, so we assumed the generation time as 2 years, referring the study of other species in Galliformes (Hung et al. 2014; An et al. 2015). Both ESS value and distribution of the approximated posterior densities were employed to evaluate the convergence and stationarity of the result from Ima2p (Sethuraman and Hey 2016).

Estimating historical trajectory of effective population size

To explore whether the effective population size of L. lerwa decreased or increased during inter-glacial period, we employed high-throughput sequencing technology to sequence 1 individual of L. lerwa, and then used Pairwise Sequentially Markovian Coalescent (PSMC) model to estimate the historical trajectory of its effective population size (Li and Durbin 2011). Considering the different populations of L. lerwa may face similar climatic fluctuations and then they may have similar demographic history, we only selected one sample to construct PSMC model. The blood of one sample from Sichuan was used for sequencing at the Novogene Sequence Center (Beijing, China). The sequencing strategy was paired-end 150 bp on the Illumina platform after constructing a fragment library (insert size of 350 bp). After discarding the adapter sequences and 5 bp bases from both 5′- and 3′-ends for each read in the Trim Galore version 0.4.2 with other default parameters, the 298,888,504 clean reads were mapped to the autosomes of Gallus gallus (ftp://ftp.ensembl.org/pub/release-98) using the BWA-MEM algorithm (Li and Durbin 2009). The samtools version 1.12 program was used to calculate the mapping rate and average coverage of all positions (Li et al. 2009). We followed the best practice suggestions of the Genome Analysis Toolkit version 4.1.4 pipeline to sort the bam file and mark the duplication reads that were produced by PCR during constructed genome library (McKenna et al. 2010).

After obtaining the bam file, we first employed the mpileup function in the samtools version 1.12 (Li et al. 2009) to produce “pileup” textual format from the bam file. In this step, we set the “coefficient for downgrading mapping quality for reads containing excessive mismatches” as 50, the “minimum mapping quality for an alignment to be used” as 50, and the “minimum base quality for a base to be considered” as 20. Then, we employed the consensus-caller method of “call” function in the bcftools version 1.10.2 (Li et al. 2009) to identify the single nucleotide polymorphism (SNP). After that, we used the “vcf2fq” function in “vcfutils.pl” script from the bcftools program to produce consensus sequences as the input file for PSMC. For the high mean genome coverage (≥18×), setting 6× as the minimum depth has no bias for the effective population size trajectory (Nadachowska-Brzyska et al. 2016), so we set the minimum depth as 6. The maximum depth was set as 48 (twice average genome coverage), following the suggestion of Li (Li and Durbin 2011). During the PSMC model, we set the maximum number of iterations to 25, the maximum coalescent times to 15, the initial theta/rho ratio was 5, and the pattern of parameters was “4 + 25 * 2 + 4 + 6.” The genome-wide neutral substitution rate (4.02 × 10−9 per site per generation) and generation time (2 years per generation) were employed to plot the trajectory of the effective population size (Zhang et al. 2014; Zheng 2015). To evaluate the variability of the PSMC result, we performed 100 rounds of bootstrapping.

Constructing ENMs

To reveal the distribution of L. lerwa during glacial and inter-glacial periods, we employed an ENM of maximum entropy algorithm (Maxent version 3.4.1) to construct the current and historical distribution during 4 periods (Phillips and Dudík 2008), including Current (1970–2000 AD), Middle Holocene (∼6 thousand years ago (kya)), Last Glacial Maximum (LGM; ∼21 kya), and Last Inter Glacial (LIG; ∼120–140 kya). Past studies suggest Sichuan Basin is the refuge for some species during glacial periods (Wu et al. 2017), and the distribution of some organisms would experience the cycling of expansion and contraction across Sichuan Basin, Hengduan Mountains, QTP, and Himalayas (Khadka and James 2017; Wu et al. 2017; Wang et al. 2018), hence, the Sichuan Basin, Hengduan Mountains, QTP, Himalayas, and their surrounding regions were selected as the study region (Supplementary Figure S1). We downloaded 19 bioclimatic variables for each period from the WorldClim version 2.1 (www.worldclim.org) with resolutions: LIG with only resolution 30 s (∼1 km), and LGM and Mid Holocene with the same resolution 2.5 min (∼5 km). The 19 bioclimatic variables represent the monthly temperature and rainfall values and extreme weather. We used Data Management Tools of Resampl as 2.5 min resolution for LIG in Arcgis version 10.2. After removing the highly correlated variables (Spearman correlation was >0.7; Supplementary Table S3), 6 bioclimatic variables (annual mean temperature [°C], temperature seasonality, precipitation of driest month [mm], precipitation seasonality, precipitation of the warmest quarter [mm], and precipitation of the coldest quarter [mm]) were used to construct ENMs. The original occurrence points of L. lerwa were from the field survey by Hongyan Yao and Nan Wang (Supplementary Table S4). Because the spatial autocorrelation would bias the results, we randomly selected the occurrence sites that were >5 km from each other. There were 52 occurrence points that were retained after random selection from 107 points for constructing ENMs (Supplementary Table S5). During running Maxent for the 4 periods, 75% of all 52 points were selected as the training set and 25% were used to test the models randomly, with subsample as the replicated run type, the maximum iterations 500, and 10 replicates. The logistic probabilities of occurrence were set as the logistic format, with the suitability values ranging from 0 to 1. The model performances were evaluated as the area under the receiver operating characteristic (AUC). To construct ENMs of the 4 periods, we set the environmental layers of the current period as the input environmental layers for all the 4 ENMs, and we set the environmental layers of the 4 periods as the projection layers for the ENMs, respectively. For the logistic output of the ENMs of 4 periods, a binary transformation was used to obtain the potential suitable distribution based on the 10th percentile training presence logistic threshold for each species with Spatial Analysis Tools of Reclassify in Arcgis version 10.2.

Results

In order to test the hypothesis that the inter-glacial isolation drives the divergence process of cold-adapted species, this study obtained 1,951-bp mitochondrial (3 loci) sequences and 16,948 bp exon sequences (26 loci) for 36 samples, including the outgroups (Supplementary Table S2). An average of 28 loci was amplified per individual, the percentage of the missing data is 2.20% in the aimed 1,044 sequences. The average length of these loci was 648 bp, and the average number of variable sites and parsimony-informative sites was 8 and 7, respectively (Supplementary Table S2). The average substitution rate of all loci was 0.26 × 10−8 per site per year (Supplementary Table S2). Among the 3 subspecies, L. l. lerwa had the highest genetic diversity (ϴw = 1.02), the genetic diversity (ϴw) of L. l. callipygia and L. l. major was 0.77 and 0.84, respectively. The genetic difference between L. l. callipygia and L. l. major was the smallest (D = 0.14 × 10−2) among the 3 pair subspecies. The D between L. l. lerwa and L. l. callipygia was 0.46 × 10−2. The D between L. l. lerwa and L. l. major was 0.50 × 10−2. The detailed genetic diversity and genetic differences on each locus can be found in Supplementary Table S6.

Both Structure and DAPC analysis suggested that there were 2 genetic populations within L. lerwa (Figure 2). The Deltak value supported 2 was the best K value (Supplementary Figure S2). The samples from L. l. lerwa were from 1 genetic population, whereas the samples from both L. l. callipygia and L. l. major were from another genetic population (Figure 2). The BIC value from DAPC analysis also suggested there were 2 genetic clusters (Supplementary Figure S3). The components of genetic populations revealed by DAPC are same as those from structure (Figure 2). The phylogenetic trees support the samples from L. l. lerwa cluster into 1 clade, except for the phylogenetic tree constructed by BI algorithm with nuclear data set (Supplementary Figures S4 and S5). Based on either algorithm, the 3 datasets all do not support the samples from either L. l. callipygia or L. l. major cluster into 1 monophyletic clade, but the samples always cluster into a paraphyletic taxon (Supplementary Figures S4 and S5). In further analysis, the 2 clades (L. l. lerwa and L. l. major, L. l. major included L. l. callipygia) revealed by Structure, DAPC and phylogenetic were used for estimating divergence time and gene flow within L. lerwa.

Figure 2

The genetic populations revealed by Structure and DAPCs. (A) There were 2 genetic populations within L. lerwa; (B) The composition of the 2 genetic clusters revealed by DAPC; (C) The genetic compositions of the different genetic populations when K = 2 and 3 (the Deltak value suggested, K = 2, is the best cluster value).

The divergence time within L. lerwa populations was in the Mid-Pleistocene. The *BEAST analysis suggested that the divergence within L. lerwa occurred at ∼0.40 (95% highest posterior density [HPD] = 0.29–0.52) mya (Figure 3A). The divergence time estimated by the IM model was 0.44 (95% confidence interval [CI] = 0.30.–0.80) mya (Figure 3B and Table 1). The isolation with migration model indicates that the effective population size of the 2 clusters was similar (Table 1). There was an asymmetric gene flow (the P-value of log-likelihood ratio test was ˂0.01) from L. l. lerwa to L. l. major during divergence process, but the gene flow from L. l. major to L. l. lerwa was not significantly larger than zero (Table 1 and Supplementary Figure S6).

Figure 3

Divergence time within L. lerwa. (A) Divergence time estimated by the *BEAST; (B) Posterior probability of the divergence time estimated by the isolation with migration model (L. l. major includes L. l. callipygia. The gray bar means the 95% HPD of the divergence time.).

Table 1

The parameters of demographic history estimated by isolation with migration model

Parameters	t0 (mya)	q0 (×106)	q1 (× 106)	q2 (×106)	m (0 > 1)	m (1 > 0)
HiPt	0.44	0.90	0.84	1.61	0.00	0.23a
95%Lo	0.30	0.68	0.65	1.01	0.00	0.06
95%Hi	0.80	1.23	1.11	2.54	0.32	0.60
MPP	3.61	2.92	3.44	1.10	7.89	2.99

(t0: divergence time between L. l. lerwa and L. l. major; q0: effective population size of L. l. lerwa; q1: effective population size of L. l. major; q2: effective population size of the ancestor population of both L. l. lerwa and L. l. major; m (0 > 1): the migration rate from L. l. major to L. l. lerwa; m (1 > 0): the migration rate from L. l. lerwa to L. l. major; hipt: the value of the bin with the highest count; 95%Lo, the estimated point to which 2.5% of the total area lies to the left; 95%hi, the estimated point to which 2.5% of the total areas lies to the right.

The gene flow was significantly larger than the null expectation that the migrate rate was zero. Lerwa lerwa major includes the previous classified L. l. callipygia). MPP, marginal posterior probability.

The mapping rate and average coverage of the sequence data used for the PSMC model were 89.01% and 24.41, respectively. There were 6,152,602 SNPs were used to construct PSMC model. Moreover, the PSMC model revealed that the effective population size of L. lerwa continued to increase as the temperature decreased, from 38,428 (±standard deviation [SD]: 2,118) to 15,763 (±SD: 1,426) years ago (Figure 4), which almost happened during the LGM period. The effective population size began to decrease approximately15,763 (±SD: 1,426) years ago as the temperature increased (Figure 4).

Figure 4

The historical trajectories of the effective population size of L. lerwa and the climate (g, generation time; u, mutation rate. The red and blue lines represent the constructed and bootstrapped curves, respectively. The black line represents the historical climate condition. The data of the climate is from Song et al. (2007). The yellow and blue bars identify the inter-glacial and glacial period, respectively. Because the parameters estimated by PSMC have small variance between 2 × 104 years ago and 3 × 106 years ago, we only show the trajectory before 104 years ago following Li and Durbin (2011).

The AUC value of the ENMs for Current, Middle Holocene, LGM, and LIG was 0.98, 0.99, 0.98, and 0.99, respectively. The suitable distribution of L. lerwa during the 2 warmer periods (Middle Holocene and LIG) was 106,324 km2 and 2377 km2, respectively, which were smaller than the corresponding distribution in the LGM (544,325 km2) and in the present (361,919 km2) (Figure 5).

Figure 5

Suitable distributions of L. lerwa during the 4 periods. (A) Suitable current distribution, represented in violet; (B) Suitable distribution during Middle Holocene, represented in cyan; (C) Suitable distribution during the LGM, represented in blue; (D) Suitable distribution during the LIG period, represented in pink (the area surrounded by the dotted line is the study area. The arrows indicate the suitable distribution in LIG. The number under each figure indicates the size of the suitable distribution. The source of the map was from https://www.resdc.cn/).

Discussion

Although lots of studies focused on the relationship between climatic fluctuations and divergence, most of them identified glacial isolation driving as the divergence process (Petit et al. 2003; Wallis et al. 2016; Neiva et al. 2018). Two recent studies on divergence processes suggest a different scenario, where inter-glacial isolation contributed to the divergence process (Dong et al. 2017; Ehl et al. 2021). One explored the avian divergence in East Asia, which was not covered by glaciation during the Pleistocene (Dong et al. 2017); the other investigated the divergence process of butterflies in the high mountain areas of Europe (Ehl et al. 2021). Here, with comprehensive evidence included divergence time (Figure 3), the historical trajectory of the effective population size (Figure 4), and ENMs (Figure 5), we aimed to determine whether glacial isolation or inter-glacial isolation drove the divergence process on cold-adapted species on the QTP. The result will have implications in understanding divergence mechanisms.

Using the multi-locus sequence data, we found 2 genetic populations within L. lerwa (Figure 2) and the genetic difference between L. l. major and L. l. callipygian is smaller than other 2 paired subspecies (Supplementary Table S6). The results suggest that there may be 2 subspecies within L. lerwa, which is consistent with the opinions of some researchers (Zheng 2017; Clements et al. 2019). A specimen of L. l. lerwa was reported by Hodgson (1833), and then Meinertzhagen (1927) described the subspecies L. l. major. Later, in 1938, Stegmann (1938) reported L. l. callipygia after he compared the specimen with L. l. lerwa; however, he did not compare the specimen with L. l. major. After comparing the morphology of the specimens, Marien (1951) indicated that L. l. callipygia is a synonym of L. l. major. Zheng et al. (1978) reported that there are 3 subspecies within L. lerwa based on the stripe of tail feathers; however, the difference was not obvious. Here, we only provide the evidence that L. lerwa has 2 genetic populations. Further analysis needs to combine morphology and geography information to delimit the subspecies within L. lerwa.

One explanation for the species divergence process on mountains is the uplift-driven diversification (Xing and Ree 2017). This event can result in a geographical barrier to gene flow and complete speciation. For example, uplifts of the QTP and its surrounding regions contributed to the divergence between Primula nutans and Primulafasciculata (Ren et al. 2018). Nevertheless, by using both isolation model and isolation with migration model, we verified that the uplift of the QTP was not the driver of the divergence within L. lerwa. The divergence time estimated by *BEAST and Ima2p is ∼0.40 and 0.44 mya, respectively (Figure 3 and Table 1). Because the QTP and its surrounding regions had reached their current altitudes by the late Pliocene (2.58 mya) (Rowley and Currie 2006; Wang et al. 2014; Deng and Ding 2015), the recent divergence within L. lerwa was unlikely to be driven by this ancient geographical event. This gives us the opportunity to test 2 additional scenarios. The first is that glacial isolation was the driver of the recent divergence. The second is that inter-glacial isolation has driven recent divergence on the QTP.

The *Beast and Ima2p employed different models to estimate divergence time, 1 is “complete isolation” model (Bouckaert et al. 2014), the other is “isolation with migration (IM)” model (Sethuraman and Hey 2016). The natural ancestral populations are not instantaneously split into several subpopulations, the gene flow would exist at the various stage of speciation (Rundle and Nosil 2005). Under climate change, the cycles of distribution expansion and contraction provide opportunity for gene flow during divergence process (Miller et al. 2006). Compared with “complete isolation” model, the IM model may be close to realistic scenarios. However, the IM model assumes that the migration rate is constant because of the splitting time, which is also unrealistic (Wilkinson-Herbots 2008). Because either model is not flawless, we integrated the divergence time calculated by the 2 models. In the IM model, the divergence time is best estimated in the scenario of gene flow equals to zero (Quinzin et al. 2015), and in the gene flow case the coalescent time of the samples predate the divergence time (Wilkinson-Herbots 2008), so we put the divergence time estimated by IM as the right boundary. The “complete isolation” model does not consider gene flow, the divergence time estimated by *Beast may posterior to the divergence time (Wilkinson-Herbots 2008), so we put it as the left boundary. Based on the posterior probability, the divergence time within L. lerwa may be 0.40–0.44 mya (Figure 3). Even though the 95% HPD from *Beast and 95% CI from IM suggest the divergence time has a broad time scale, setting 0.40–0.44 mya as the divergence time has the highest posterior probability. Electron spin resonance dating of glacial deposits in the headwaters of Urumqi river in Tianshan Mountains suggests there was a glaciation during 0.45–0.48 mya, called as Zhongliangan glaciation (Zhou et al. 2001; Shi et al. 2006). The divergence within L. lerwa was ∼0.40–0.44 mya (Figure 3), happening in inter-glacial period after Zhongliangan glaciation.

We further estimated the historical trajectory of the effective population size of L. lerwa and found its contraction during inter-glacial period (Figure 4). The result suggests the climate fluctuation may have acted on the L. lerwa. During glacial period, ∼16–21 kya, the L. lerwa had maximum effective population size (Figure 4). As the temperature became warming, the effective population size of L. lerwa gradually declines . Consistent with our result, the organisms that distributed on the QTP and its surrounding regions, such as blue eared pheasant Crossoptilon auritum, also have a fluctuating effective population size (Wang et al. 2021). A cooling climate is thought to make the cold-adapted species have large population size, and if the climate becomes warming then the population size of cold-adapted species may decrease. Compared with the cold-adapted species, the temperate species may have opposite scenario. The temperate species may have small population size facing extensive glaciations and then expand their population size following the retreat of the glaciations.

After we re-constructed the trajectory of effective population size, we compared the potential distributions of L. lerwa among the 4 historical periods. The potential distribution of L. lerwa during the 4 periods revealed the L. lerwa expanded its distribution during glacial periods and contracted its distribution during the inter-glacial period (Figure 5). During LIG, the L. lerwa was distributed in the high elevations of southern Hengduan Mountains and the south of the Himalayas. The temperature during the LIG was ∼2 °C warmer than current, and the sea level was ∼6 m above the current sea surface (Otto-Bliesner et al. 2006). The warm climate may have pushed the cold-adapted species to the higher elevations, such as Blood pheasant Ithaginis cruentus, and Tibetan snowcock T.tibetanus (Zhan et al. 2011; An et al. 2015). The divergence time, 0.40–0.44 mya, within L. lerwa, located in the marine isotope stage (MIS) 11 period (Zhou et al. 2006). The temperature of MIS11 period was similar to that of LIG (MIS5; Zhou et al. 2006). The similar climate conditions between MIS11 and MIS5 suggest the L. lerwa distributed in the isolated restricted region at the initial divergence period. The inter-glacial isolation may induce divergence process within L. lerwa.

Combining the evidence of divergence time, effective population size, and potential distribution, we infer that the inter-glacial isolation has driven the divergence within L. lerwa. The cold-adapted species may have expanded its distribution with the development of glaciers but was discontinuously distributed on the high altitudes of the mountains after glacier melting at low altitudes. The inter-glacial isolation may drive the divergence process of cold-adapted species, and the divergence process of temperate species may be driven by glacial isolation.

Click here for additional data file.