Population genetics analysis of Tolai hares (Lepus tolai) in Xinjiang, China using genome-wide SNPs from SLAF-seq and mitochondrial markers

The main topic of population genetics and evolutionary biology is the influence of the ecological environment, geographical isolation, and climatic factors on population structure and history. Here, we estimated the genetic diversity, genetic structure, and population history of two subspecies of Tolai hares (Lepus tolai Pallas, 1778), L. t. lehmanni inhabiting Northern and Northwest Xinjiang and L. t. centrasiaticus inhabiting Central and Eastern Xinjiang using SNP of specific-length amplified fragment sequencing (SLAF-seq) and four mitochondrial DNA (mtDNA). Our results showed a relatively high degree of genetic diversity for Tolai hares, and the diversity of L. t. lehmanni was slightly higher than that of L. t. centrasiaticus, likely due to the more favorable ecological environment, such as woodlands and plains. Phylogenetic analysis from SNP and mtDNA indicated a rough phylogeographical distribution pattern among Tolai hares. Strong differentiation was found between the two subspecies and the two geographical groups in L. t. centrasiaticus, possibly due to the geographical isolation of mountains, basins, and deserts. However, gene flow was also detected between the two subspecies, which might be attributed to the Tianshan Corridor and the strong migration ability of hares. Tolai hare population differentiation occurred at approximately 1.2377 MYA. Population history analysis based on SNP and mtDNA showed that the Tolai hare population has a complex history and L. t. lehmanni was less affected by the glacial event, possibly because its geographic location and terrain conditions weaken the drastic climate fluctuations. In conclusion, our results indicated that the joint effect of ecological environment, geographic events, and climatic factors might play important roles in the evolutionary process of L. t. lehmanni and L. t. centrasiaticus, thus resulting in differentiation, gene exchange, and different population history.


Introduction
Climate factors and habitat environment will profoundly impact the evolution of biological populations, thus leaving historical traces in the genetic diversity, population structure, distribution pattern, and other aspects of today's populations. For example, climate changes profoundly affected phylogeographic structure and the evolutionary history of brown hares (L. europaeus) by isolating populations in the distinct refugia where they adapted and differentiated in allopatry, leading to genome incompatibilities (Djan et al., 2017;Minoudi et al., 2018;Giannoulis et al., 2019). The southwest of Tarim Basin in Xinjiang, China, as the origin of rivers in the basin, was a glacial refugia for Yarkand hares (L. yarkandensis) during the Quaternary climate oscillations, providing a suitable environment for maintaining the relatively high genetic diversity of this species (Shan et al., 2011). In addition, genome-wide SNPs have confirmed differentiation between the southwest and northern populations (Ababaikeri et al., 2021) and higher diversity in the former population.
Lepus species have a wide distribution and can survive in various complex terrestrial habitats (Ben Slimen et al., 2018). Hence, they are indispensable for local ecosystems as part of the food chain. The Tolai hare (Lepus tolai Pallas, 1778) has sandy yellow, brownish gray, or gray dorsal pelage with a dark ripple (Smith et al., 2018). It occupies in various habitats, including desert, semi-desert, mountain steppe, foreststeppe, rocky habitats, and grasslands, and ranges from low to high elevations. Their strong ability to survive in diverse habitats makes them a good model for studying animal adaptation to the environment. However, since the low differences in morphology, and hybridization among hares, their classification has been challenging . Tolai hares were once classified as Cape hares (L. capensis) or Brown hares, but studies based on molecular biology (Wu et al., 2005;Wang and Yang, 2012;Shan et al., 2020b) and skull measurements (Cheng et al., 2012) did not support the classification of "L. capensis" in China as L. capensis. Therefore, the original Xinjiang's "L. capensis" population was divided into L. tibetanus and L. tolai, and the population distributed in the north of the Tianshan Mountains was classified as L. tolai (Shan et al., 2020a).
Xinjiang, China, is in the hinterland of the Eurasian continent, with a dry climate and unique geological structure called "three mountains nip two basins". Tolai hare distributes in the vast area of northern, central, and northwest in Xinjiang (Smith and Xie, 2008;Smith et al., 2018) including Altai Mountains and Junggar Basin in the north, Turpan-Hami Basin in the east, and the Tianshan Mountains crossing the central-eastern part of Xinjiang. Recent studies have shown two subspecies of Tolai hare in Xinjiang, with only slight different appearances (Smith et al., 2018;Shan et al., 2020a), among which L. t. lehmanni mainly inhabits the northern and northwestern regions. Except for Junggar Basin, these areas have many rivers, abundant water resources, and relatively humid climates. L. t. centrasiaticus mainly distributes in the central and eastern Tianshan Mountains. These areas are relatively dry with little rainfall and are dominated by arid and desert habitats.
Early molecular biological studies of Tolai hares were based on the misclassification of L. tolai into L. capensis. For example, Wu et al. (Wu et al., 2005) studied the phylogenetic relationships, biogeographic distribution, and species origin patterns between Chinese hare groups, including "L. capensis" in central Xinjiang. Wang et al. (Wang and Yang, 2012) sequenced the entire mitochondrial genome of Cape hares. They reconstructed phylogenetic relationships in genus Lepus based on CYTB, including so-called "L. capensis" distributed in Xinjiang. Liu et al. (Liu et al., 2011) used four mitochondrial DNA (mtDNA) fragments and nuclear gene to demonstrate that frequent introgression occurred through historical and recent interspecific hybridization among six Chinese hare species, including those in northern Xinjiang. Wu et al. (Wu et al., 2011) found extensive bidirectional mitochondrial DNA and SRY gene introgression in hybrids of Yarkand hare and Xinjiang "L. capensis". Recently, some studies have mapped the full mitochondrial genome of the Xinjiang Tolai hare . However, there has been no comprehensive evaluation of their genetic diversity and structure related to habitat and climate changes, and research on the biology of this species is scarce.
Specific-Locus Amplified Fragment sequencing (SLAFseq) is a high throughput, high-resolution, and low-cost marker development technology that has emerged recently (Sun et al., 2013). This technique focuses on finding single nucleotide polymorphisms (SNPs), an abundant form of genetic variation, in an economical way Wang et al., 2016;Ali et al., 2018;Qin et al., 2019). SNPs are found throughout the genome, and their distribution can reflect the population's genetic variation. SLAF-seq has been used to analyze several species' genetic diversity and phylogenetic structure Chang et al., 2019;Qin et al., 2019;Zhang J. et al., 2020a;Fang et al., 2020;Ababaikeri et al., 2021). In addition, mtDNA provides a different perspective of population genetic structure because it is maternally inherited and generally lacks intermolecular recombination (Allendorf, 2017).
This study used SLAF-seq to identify genome-wide SNP markers and combined them with four mtDNA markers Frontiers in Genetics frontiersin.org 02 (COI, ND4, CYTB, and D-LOOP) to evaluate the genetic diversity of two subspecies of Tolai hare in different habitats and reveal the effects of ecological environment, geographical isolation, and climate change on population structure characteristics and population demography history. This study will help to understand the evolutionary history of this species and provide essential data to maintain the biodiversity and stability of the Xinjiang ecosystem.

Sampling and DNA extraction
Muscle or skin tissue samples were collected from 106 Tolai hares from 12 geographic populations in northern, northwestern, central, and eastern Xinjiang from 2008 to 2019 (Table 1 and Figure 1). The 81 L. t. lehmanni samples included 59 samples from Altay Prefecture (24 from Altay, 26 from Fuhai, three from Burqin, four from Habahe, and two from Qinghe), three from Tarbagatay Prefecture, 13 from Bortala Mongol Autonomous Prefecture (including 12 from Jinghe and one from Wenquan), and six from Ili. Twenty-five L. t. centrasiaticus samples included 11 from Dabancheng, two from Tuokexun, and 12 from Kumul. The samples were divided into northern, northwest, central, and eastern groups, and the geographical details of the sampled populations are shown in Table 1. All Tolai hare samples were used for mitochondrial DNA analysis. For SLAF-seq analysis, samples that were stored for a long time, severely degraded, or whose DNA quality was too low to be sequenced were eliminated, and a total of 36 Tolai hare samples in two subspecies remained. These included samples from Altay, Tarbagatay Prefecture, Ili, Dabancheng, and Tuokexun (Table 1). Some samples in this study were confiscated from poachers and provided to us by local forestry bureaus, while others came from hares that died of natural causes. All experimental protocols involved in this study were approved by the Institutional Animal Care and Use Committee of the College of Life Science and Technology, Xinjiang University, Urumqi, China.
Muscle samples were preserved in sterile tubes with anhydrous alcohol at -80°C until total genomic DNA extraction using a DNA tissue extraction kit. Genomic DNA integrity was determined using 1.0% agarose gel electrophoresis.

Construction of SLAF-seq library and high-throughput sequencing
The domestic rabbit (Oryctolagus cuniculus) OryCun 2.0 genome (Lindblad-Toh et al., 2011) from the National Center for Biotechnology Information (NCBI: ftp://ftp.ncbi. nlm.nih.gov/genomes/all/GCF/000/003/625/GCF_000003625.3_ OryCun2.0/) served as the reference genome for genome electronic enzyme digestion prediction. The selection principle for the enzyme digestion scheme was as follows: the proportion of restriction fragments located in the repetitive sequence was as low as possible, the restriction fragments were distributed as evenly as possible across the genome, the length of restriction fragments was highly consistent with the experimental study system, and the number of restriction fragments obtained matched the expected number of tags. A RsaI-EcoRV-HF ® restriction enzyme was used to digest the genomic DNA. From the constructed SLAF library (Sun et al., 2013;Zhang et al., 2013) DNA fragments 314-344 base pairs (bp) long were selected for sequencing on an Illumina HiSeq 2,500 system (Illumina, Inc., San Diego, CA, United States) at Beijing Biomarker Technologies Corporation (Beijing, China). The Oryzasativa ssp. japonica genome (http://rapdb.dna.affrc.go.jp/ ) was selected as the positive control for sequencing, and SOAP v2 software (Li R. et al., 2009b) was used to calculate the enzyme digestion efficiency, fragment insertion distribution, and alignment efficiency to evaluate the reliability of the enzyme digestion experiment and the accuracy of library construction.

Quality control and SNP calling
Reads were obtained for each sample, and after filtering the connectors of the sequencing reads, the sequencing and data volumes were evaluated. We calculated the ratio of highquality reads based on two key quality indicators, Q30 (read quality score of 30, indicating a base sequencing error probability of 0.1%) and GC content. All sample reads were mapped to the OryCun 2.0 genome sequence using BWA v0.7.5a-r405 (Li and Durbin, 2010). Then, we mined the SLAF tags according to the restriction fragment size defined by the enzyme digestion scheme. SNP calling was performed using GATK v3.3.2 (McKenna et al., 2010) and SAM tools v0.1.18 (Li H. et al., 2009a), and the intersection of SNP markers called by the two packages was chosen to build the original SNP dataset. Plink v1.07 (Purcell et al., 2007) was used to filter high-quality SNPs according to the criteria of site information integrity (INT) ≥ 0.5 and minor allele frequency (MAF) ≥ 0.05. Finally, the selected high-quality SNPs were used for further analysis.
Frontiers in Genetics frontiersin.org
Phylogenetic analysis based on SNPs was reconstructed using the Neighbor-Joining (NJ) method and Maximum-Likelihood (ML) method. The NJ phylogenetic tree was performed in Mega X (Kumar et al., 2018) with 1,000 bootstraps, and the ML phylogenetic tree was performed in IQ-TREE (Nguyen et al., 2015) with 3,000 bootstraps on the Cipres (https://www.phylo. org/index.php/). mtDNA-based phylogenetic analysis was conducted on the Cipres using the Bayesian Inference (BI) method by MrBayes 3.2.7a (Ronquist et al., 2012) with five million Bayesian Markov chain Monte Carlo (MCMC) generations and ML tree by RAxML-NG (Kozlov et al., 2019) with 1,000 bootstraps. The rabbit (O. cuniculus) genome was used as the outgroup (GenBank accession: GCF_000003625.3 for SNP and AJ001588.1 for mtDNA). ModelFinder was used to identify the best fit base-pair substitution model according to the Bayesian information criterion (BIC) using Phylosuite v1.2.2 (Zhang et al., 2018), and the general time-reversible (GTR) model was selected as the optimal model for phylogenetic analysis of the SNPs and mtDNA concatenated dataset. Principal component analysis (PCA) was conducted using Plink v1.07 (Purcell et al., 2007). The population structure among the Tolai hare populations was inferred using ADMIXTURE v1.3.0 (Alexander et al., 2009). Treemix v.1.13 (Pickrell and Pritchard, 2012) was used to infer multiple population splitting and mixing events using genome-wide allele frequency data. The median-joining network was conducted by Popart (Leigh and Bryant, 2015). We performed a hierarchical analysis of molecular variation (AMOVA) using Arlequin ver3.5 (Excoffier and Lischer, 2010).
Divergence times were estimated based on mtDNA using Beast v1.10.4 . Four points were calibrated to build the time tree: the split of Lepus (approximately 11.57 million years ago, MYA), the divergence time of L. americanus (approximately 8.6 MYA), the divergence time of L. europaeus (approximately 1.84 MYA) (Ge et al., 2013), and Tolai hare fossil calibration (0.78 MYA) (Erbajeva and Alexeeva, 2000). Best substitution models according to BIC (HKY model for COI, ND4, and CYTB and GTR model for D-LOOP) were found using modelFinder in Phylosuite v1.2.2 (Zhang et al., 2018), and uncorrelated relaxed lognormal clock with a prior coalescent tree of constant size was used. The Markov chain Monte Carlo analysis was run thrice for 1 × 10 8 generations, sampling every 1,000 generations. Tracer v1.7.2     was used to summarize the tree data, and the first 10% of trees were discarded as a burn-in. The tree and divergence times were displayed and edited in Figtree v1.4.3 (http://tree.bio.ed.ac. uk/software/figtree/). SMC++ v1.15.5 (Terhorst et al., 2017) was used to check the population history of the two subspecies based on SNPs. Tajima's D and Fu's Fs analyses were performed to test neutrality based on mtDNA in Arlequin ver3.5 (Excoffier and Lischer, 2010). The mismatch distributions of pairwise sequence differences for the Tolai hare populations were estimated using DnaSP6 (Rozas et al., 2017). To estimate changes in effective population size through evolutionary time, we explored population history by constructing Extended Bayesian Skyline Plots (EBSP) in Beast. v2.4.7 (Bouckaert et al., 2014) with mtDNA. The prior tree was set as the Coalescent Extended Bayesian Skyline and sampled every 1,000 steps for 1 × 10 8 steps. The ESS values for all the parameters were assessed using Tracer v1.7.2 .

SLAF-seq and SNP discovery
The rabbit genome served as the reference genome for predicting electron enzymes and identified a restriction fragment length of 314-344 bp, ultimately defined as a SLAF tag. Sequencing results for the positive control indicated that the efficiency of paired-end comparison was 96.84%, and the enzyme digestion efficiency of the control was 94.74%, indicating that the process was normal and reliable.
We reconstructed phylogeny to explore genome-wide relationships among Tolai hare populations. The topology of ML ( Figure 2A) and NJ trees (Additional file 3: Supplementary Figure S1) was consistent. The phylogenetic tree results showed that the Tolai hares analyzed in this study based on SNPs were divided into two main clusters with high confidence, including three clades. The first branch was located at the tree's root, containing two samples from the TKX population. The second branch consisted of samples from the DBC population and three individuals from the ALT population. These two clades were grouped into one cluster because it included most samples from the central group of L. t. centrasiaticus, except for three individuals from the ALT population of L. t. lehmanni. The third branch consisted mainly of samples from the northern and northwestern groups, including samples from the ALT, TC, and YL populations of L. t. lehmanni, and one sample from the DBC population. The clustering relationships among populations were also evident in the PCA, in which L. t. centrasiaticus individuals were clustered, and L. t. lehmanni individuals were clustered together except for individuals from the TC population ( Figure 2B). The population relationships were largely consistent with the geographical distribution of the samples, consistent with the phylogenetic tree.
To evaluate population structure, the ADMIXTURE was performed ( Figure 2C, Additional file 4: Supplementary Figure  S2). Since the cross-validation errors at K = 2 is close to the lowest cross-validation errors at K = 1 (Additional file 5: Supplementary Figure S3), the assessment of the population structure indicated a clear subdivision of two main ancestral lineages when K = 2. Except for DBC3, most samples in L. t. lehmanni formed an ancestral cluster (shown in green in Figure). Most samples in L. t. centrasiaticus, and three ALT individuals formed another ancestral cluster (shown in red in the image). In addition, fifteen samples have mixed lineages.
Genetic differences among the geographical groups were examined using AMOVA (Table 3), and the results indicated moderate genetic differences among the geographical groups (7.65%, F ST = 0.0765, p < 0.01). When individuals were pooled into the northern group and northwest group, only 6.57% (p < 0.01) of the variability was partitioned among groups, and differentiation among groups was 0.0657 (p < 0.01). When individuals were pooled into the L. t. lehmanni and L. t. centrasiaticus, only 6.57% (p < 0.01) of the variability was partitioned among populations, and differentiation among subspecies was 0.0657 (p < 0.01).
To further infer population gene flow, we performed a Treemix analysis based on SNPs ( Figure 2D). Three migrations were detected. One event began from the TKX to the ALT population (migration weight 0.2775). The second began from the ALT population to the DBC population (migration weight 0.184,875). The third event is from the DBC to the YL population (migration weight 0.0999,906).

Genetic diversity and population structure of mtDNA
Sequencing of the four mitochondrial DNA segments (COI, ND4, CYTB, and D-LOOP) with an alignment length of 2,885 bp in 106 Tolai hare individuals resulted in 99 haplotypes. Nucleotide diversity ranged from 0.0295 (eastern group) to 0.0449 (northwest group) ( Figure 3B). The π value of L. t. lehmanni (0.0410) was slightly higher than that of L. t. centrasiaticus (0.0355) ( Figure 3B). The total π value of the Tolai hare in Xinjiang was 0.0442. To explore the phylogenetic structure, ML ( Figure 3A) and BI evolutionary trees (Additional file 6: Supplementary Figure S4) were constructed, and the topological structures of the two trees were consistent. The phylogenetic tree based on combined haplotype sequences of four mitochondrial genes showed that the Tolai hare samples analyzed in this study contained four branches and could be further divided into two main clusters ( Figure 3A). The first cluster was located at the tree's root and predominantly comprised of L. t. centrasiaticus individuals from the central and eastern groups. Each individual from the QH, ALT, and FH populations were included. The second cluster was comprised of L. t. lehmanni individuals from the northern and Frontiers in Genetics frontiersin.org 08 northwest groups, and two individuals from the DBC population. In addition to the evolutionary tree, a median-joining network was constructed to further elucidate the phylogenetic relationships among haplotypes ( Figure 3C). The haplotype grouping in the median-joining network was consistent with the clustering in the evolutionary tree. The network verified the existence of two distinct clades separated by several mutational steps.
Genetic differentiation among the groups was next examined using AMOVA. The results showed that when pooling individuals into northern, northwest, central, and eastern groups, the genetic variation among groups was 26.00% (p < 0.01), and differentiation was 0.2600 (p < 0.01) ( Table 4). Genetic variation and differentiation among groups was lower when individuals were pooled into the northern and northwest group (15.58%, F ST = 0.1558, p < 0.01) than when individuals were pooled into the central and eastern group (27.24%, F ST = 0.2724, p < 0.01) ( Table 4). In addition, when individuals pooled into L. t. lehmanni and L. t. centrasiaticus, genetic variation among populations was 23.76% (p < 0.01), and F ST was 0.2376 (p < 0.01).

Divergence time estimation and population history
We used mtDNA data to construct a Tolai hare differentiation time-merging tree. The tree showed that Tolai hare population differentiation occurred at approximately 1.2377 MYA, and the eastern HM population showed the earliest differentiation. The divergence between L. t. lehmanni and L. t. centrasiaticus occurred at approximately 1.1898 MYA ( Figure 4A).
To estimate the historical population dynamics of Tolai hare in Xinjiang, we performed SMC++ analysis to track changes in effective population sizes (Ne) over time based on SNPs. The dynamics of historic Ne for L. t. lehmanni and L. t. centrasiaticus are shown in Figure 4B. The effective population size of L. t. centrasiaticus decreased in the last glacial maximum (LGM), while L. t. lehmanni population remained relatively stable. Both subspecies expanded during the interglacial period after the last glacial period (LGP), especially L. t. lehmanni. About 0.001-0.0015 MYA, however, populations of both subspecies declined. In addition, the Tajima's D, Fu's Fs, mismatch distribution, and EBSP based on mtDNA were used to analyze the population history of Tolai hares. The neutrality test results showed that Tajima's D and Fu's Fs values for L. t. lehmanni were 2.06761 (p = 0.9800 > 0.05) and -7.73931 (p = 0.0480 < 0.05), and -0.00167 (p = 0.5720 > 0.05) and -3.95338 (p = 0.0310 < 0.05) for L. t. centrasiaticus. The mismatch distribution analysis revealed that L. t. lehmanni and L. t. centrasiaticus showed multi-peak curves (Additional file 7: Supplementary Figure S5). The demographic scenario for Tolai hare populations determined through EBSP analysis suggested a more pronounced population expansion of L. t. lehmanni beginning from 0.02 MYA, and the L. t. centrasiaticus population showed a less distinct tendency to expand ( Figure 4C).

Discussion
The relatively high genetic diversity of Tolai hare populations The amount of genetic diversity reflects the evolutionary potential of a species, and populations with more genetic diversity are expected to adapt better to environmental changes such as climate change, habitat loss, overharvesting, invasive species, and disease than populations with low genetic diversity (Kardos, 2021). In this study, we used a set of SLAF-seq genome-wide SNP and mtDNA markers to estimate the genetic diversity of Tolai hares in Xinjiang. At the genome-wide level, the genetic diversity of Tolai hares based on SNPs was higher (Table 2) (π = 0.05926, He = 0.37412, Ho = 0.40226, PIC = 0.2974) than that of other reported species, including rabbits (PIC = 0.2-0.2281, He = 0.2511-0.2857, Ho = 0.3072-0.3418, (Ren et al., 2019), and Yarkand hares (π = 0.0655, He = 0.3130, Ho = 0.2852, PIC = 0.2543 (Ababaikeri et al., 2021). The estimated nucleotide diversity of Tolai hares (π = 0.0442) based on mtDNA was also relatively high compared to other Lepus taxa, such as brown hares (π D-LOOP = 0.030, (Minoudi et al., 2018), Yarkand hares (π D-LOOP = 0.033, π CYTB = 0.008, (Shan et al., 2011), and Italian hares (L. corsicanus, π D-LOOP = 0.018, (Pierpaoli et al., (Pironon et al., 2017). Larger effective population sizes may have also contributed to this (Kimura, 1983;Hague and Routman, 2016). Frontiers in Genetics frontiersin.org 10 SNPs and mtDNA-based genetic diversity analysis of the subspecies showed that L. t. lehmanni polymorphism level was higher than that of L. t. centrasiaticus, likely due to a more favorable living environment and distribution (Table 2 and Figure 3B). In Xinjiang, L. t. lehmanni is mainly distributed in the northern and northwestern regions, which are more humid and include vast forests, grasslands, and croplands, providing more water and food for survival. While L. t. centrasiaticus is mainly distributed in the central and eastern Xinjiang, including the Turpan-Hami basin, which is primarily the Gobi desert; thus, survival there is likely more challenging.
The coexistence of genetic differentiation and gene flow in Tolai hare populations Geological formation events can cause physical barriers to disperse species, such as mountains and rivers, leading to divergence and speciation (Chaves et al., 2011). Thus, geographical isolation is an important factor in genetic differentiation (Duncan et al., 2015;Gao et al., 2021). The effects of geographic isolation on genetic structure have been reported in several species, such as Sarcophanops (Campillo et al., 2020), Orestias ascotanensis (Cruz-Jofré et al., 2016), Weigela coraeensis (Yamada and Maki, 2012), L. capensis, and L. europaeus (Ben Slimen et al., 2008). In this study, phylogenetic trees based on SNP and mtDNA markers (Figures 2A, 3A) showed that Tolai hares were divided into two main clusters (L. t. lehmanni and L. t. centrasiaticus). Most samples from the same or adjacent regions clustered together, and their clustering relationships were consistent with their geographic distribution, indicating a rough phylogeographic distribution pattern. This was also supported by the network analysis ( Figure 3C) and the PCA ( Figure 2B). In addition, combined with the F ST and significant p-values in AMOVA based on SNPs and mtDNA (Tables 3 and 4), differentiation occurred between subspecies L. t. lehmanni and L. t. centrasiaticus. These two subspecies distribution areas span across the Tianshan Mountains in Xinjiang. We speculate that the high mountains, the barren basins and deserts, and long geographical distances serve as geographical barriers that limit the dispersal and communication between two subspecies, leading to genetic divergence. In addition, the different populations may have evolved different ecotypes with different feeding habits or habitats, hindering gene exchange and leading to differentiation (Sarabia et al., 2021).
A strong differentiation was also found between the Central and Eastern groups in L. t. centrasiaticus according to the phylogenetic analysis and AMOVA results from mtDNA (Table 4). Although the Dabancheng (DBC) and Tuokexun (TKX) of the Central Group are adjacent to the Kumul (HM) of the Eastern Group, they are separated by the huge Turpan-Hami Basin and the Gobi Desert and are relatively distant geographically. Moreover, the food shortage caused by drought, little rain, and hot habitat further prevent the migration, spread, and exchange of hares, thus promoting differentiation.
In this study, the genetic structure analysis showed that geographical isolation caused by geographical barriers, such as distance, mountains, basins, and deserts, as well as the ecological environment of the habitat, affected the genetic structure of hare populations. This also occurs in Yarkand hare in Xinjiang. The population genetics of Yarkand hares based on SNPs (Ababaikeri et al., 2021) and mtDNA (Shan et al., 2011) showed a systematic geographical distribution pattern and genetic differentiation between the north and southwest populations, which may have been caused by geographical isolation and different living environment. It is concluded that geographical isolation and complex habitats can affect hare populations' genetic structure and differentiation.
On the other hand, parapatric and sympatric subspecies or populations are accessible to gene flow, thus affecting the genetic structure and evolutionary history. Despite clear population differentiation and significant variation among subspecies and geographical groups, our phylogenetic analyses, PCA plot, ADMIXTURE, and Treemix results revealed a certain degree of lineage admixture and gene flow between subspecies. This is likely due to the strong adaptability to environmental changes, large effective population sizes, and relatively strong migration capability, promoting gene exchange among populations (Ababaikeri et al., 2021). Moreover, since ancient times, the Tianshan Corridor in the Silk Road has connected the northern and southern Xinjiang (Chen, 2014). Although Tianshan Mountains lie across Xinjiang, the biological corridors may contribute to admixture. Gene flow also has an important impact on population differentiation. Generally speaking, differentiation occurs when there is limited gene flow between populations or strong natural selection is sufficient to overcome the dilution effect of gene flow between populations on population differentiation (Nosil, 2008;Li et al., 2014). Many studies reported gene flow coexisting in species or population differentiation (Nadachowska and Babik, 2009;Ababaikeri et al., 2021). In short, the geographical factors and the properties of hares might contribute to the coexistence of genetic differentiation and exchange in Tolai hares.

The complex population history of Tolai hare populations
Climate-driven environmental changes during the Pleistocene influenced the evolution of many terrestrial species (Meiri et al., 2020). However, the influence of climate fluctuation on the Tolai hares population history has not been studied. In this study, neutrality tests, nucleotide mismatch distribution, EBSP analysis based on mtDNA, and SMC++ analysis based on SNP were used to explore the population history of Tolai hares in Xinjiang. The integrated results indicated that the Tolai hare population has a complex history. SMC++ shows that L. t. centrasiaticus effective population size decreased during the LGM while L. t. lehmanni population remained stable. L. t. lehmanni was less affected by LGM might because it is distributed between the mountains, which cushion drastic climate fluctuations (Yuan et al., 2019). In the interglacial period after the LGP, the population of both subspecies had a small increase with the growth of L. t. lehmanni being more evident, consistent with the Frontiers in Genetics frontiersin.org EBSP analysis based on mtDNA. The population of the two subspecies decreased from about 0.001 to 0.0015 MYA. Previous studies have shown that early humans in Xinjiang expanded into the present region during this period (Tan et al., 2022), and that the early diet of humans in northwest China included wild animals (Ren et al., 2017). Meanwhile, previous researchers believed that Lepus could be influenced by human activities (Ge et al., 2013). Therefore, We hypothesize that human activity during this period led to population decline in both subspecies. Unlike other species in Europe and eastern China, such as the Brown hare (Minoudi et al., 2018), Panda (Ailuropoda melanoleuca) (Zhao et al., 2013), and Lizards (Shinisaurus crocodilurus) (Xie et al., 2022), our comprehensive indexes of two markers showed that L. t. lehmanni population remained stable, even in the LGM (0.0265-0.019 MYA) (Clark et al., 2009;Liang et al., 2017). Other studies have shown that among the Tibetan Plateau species, population size did not decline in the LGM (Liang et al., 2017), implying that the geographical distribution of species and climate oscillation indeed play an important role in population history.
Tolai hare differentiated at 1.2377MYA, and the divergence between subspecies occurred between 1.1183 and 1.1898MYA, later than the formation time of the Tianshan Mountains and the Turpan-Hami Basin (Li et al., 2006). Therefore, mountains and basins act as geographical barriers for gene exchange between hare populations, resulting in apparent differentiation between the two subspecies and within L. t. centrasiaticus.
However, the limitations of the SLAF-seq and the mtDNA fragments used in population history analysis should also be considered. In the future, whole genome resequencing and complete mitogenome sequencing will be performed to obtain more reliable results on the Tolai hare population history.

Conclusion
This paper used SNP and mitochondrial markers to explore the effects of geological formation and environmental and climate factors on the Tolai hare. We found that Tolai hares have a high genetic diversity due to their strong adaption and migration ability. In addition, geological events such as the formation of mountains, basins, and other geographical factors led to differentiation and gene flow between subspecies L. t. lehmanni and L. t. centrasiaticus. Thus, it is due to the different geographical and geological conditions resulting in a relatively steady climate, making L. t. lehmanni stable and less affected by LGM compared to L. t. centrasiaticus.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/bioproject/PRJNA850843. The Accession number of mtDNA data for the Xinjiang Tolai hare in this study see the Additional file 2: Supplementary Table S2.

Ethics statement
All experimental protocols involved in this study were approved by the Institutional Animal Care and Use Committee of the College of Life Science and Technology, Xinjiang University, Urumqi, China.