A total of 79 E. miletus specimens were collected from September 2018 to January 2019 in Yunnan Province, including ALS (n = 15), JC (n = 22), LJ (n = 20), and DQ (n = 22). The geographic location, climatic characteristics, and sample number of each collection site were shown in Table 1 . We euthanize E. miletus using the anesthesia method. And then employed the phenol/chloroform method to extract genomic DNA from the liver (Arnason et al., 2018). We selected the reference genome of Microtus ochrogaster for electron enzyme digestion. We used SLAF-seq genomic data of 79 individuals collected from four regions and re-sequencing data of one individual from each region. 79 DNA samples were digested with Rsal, and sequenced with fragment lengths of 464–494 bp were defined as SLAFtag, and the obtained SLAFtag was purified by PCR amplification and then selected as a target fragment using 1.5% agarose gel electrophoresis and sequenced using pair-end sequencing on the Illumina HiSeq 2500 platform. A DNA small fragment library with an insert fragment size of 300–400 bp was constructed for resequencing in one randomly selected E. miletus from each population.

Reads containing more than 10% unrecognized nucleotides (N) were removed from the SLAF-seq genome, and reads with over 50% low-quality bases (Phred quality score ≤ 5) were discarded. For re-sequencing data, the raw short read sequences were filtered using SOAPnuke v.1.5.6 software (Li et al., 2009). Using M. ochrogaster as the reference genome, the re-sequencing data and SLAF-seq genomic data were mapped to the reference genome using the MEM module of BWA v.0.7.17-r1188 with default parameters (Bennett et al., 2022), respectively. The SAM files were subsequently generated using Samtools v.1.13 with default parameters (Li et al., 2009) for sorting and duplication in Picard v.3.1.1 (https://github.com/broadinstitute/picard/releases).

The Select Variants method in GATK v.4.4.0.0 software was used for SNP calling with default parameters across the all individuals (McKenna et al., 2010). We performed a filtering step with the following set of parameters: QD < 2.0, MQ < 40.0, FS > 60.0, SOR > 3.0, MQ Rank Sum < −12.5, Read Pos Rank Sum < −8.0 (Yuan et al., 2023). To obtain more reliable SNP loci, we further filtered using vcftools –max-missing with the parameter size set to 0.2 (Danecek et al., 2011).

We constructed the phylogenetic tree based on SLAF-seq genomic data and re-sequencing data using the maximum likelihood method implemented in RAxML v8.2.4, with the ascertainment bias correction and M. ochrogaster as the outgroup (Stamatakis, 2014). Before constructing the phylogenetic tree, we used the Python script vcf2pliphy.py for alignment and converted to a pliphy format recognizable by RAxML (Ortiz, 2019). We performed a Bayesian approach using re-sequencing data with M. ochrogaster as the outgroup in MrBayes v.3.2.7 with MCMC runs, iterated 1,000,000 times, and sampled every 1000 generations (Ronquist et al., 2012). The vcf file was converted to nexus format using the Python script vcf2pliphy.py before constructing the Bayesian tree, and then the nexus was converted to nexus format recognized by MyBayes v.3.2.7 using ALTER v.1.3.4 (Glez-Peña et al., 2010). The ML tree and the Bayesian tree obtained from the SLAF-seq genomic data and re-sequencing data were visualized in Figtree v.1.4.4 (Ngugi et al., 2023).

We employed several methods to test for admixture and genetic affinity among various populations. Outgroup f3-statistics were calculated using M. ochrogaster as the outgroup with the qp3pop program in Admixtools v.7.0.2, in the form of f3 (Pop1, ALS; M. ochrogaster), to assess the relative genetic similarity of ALS populations to other populations (Patterson et al., 2012; Raghavan et al., 2014). High f3 values indicate a strong degree of shared genetic history between populations (Patterson et al., 2012). Subsequently, we conducted an admixture f3-statistical analysis using the 3pop module in TreeMix v.1.13 (Alexander et al., 2009; Pickrell and Pritchard, 2012; Harris and DeGiorgio, 2017). The f3-statistic explicitly tests whether a taxon of interest (C) is the result of admixture between two others (A and B), considering the product of allelic differentials between (c − a)(c − b). Negative values suggest that the allelic frequencies of c are intermediate across many populations, which is consistent with admixture (Lopes et al., 2023). F4-statistics can visualize genetic drift shared between taxa (Harris and DeGiorgio, 2017). For an f4 (A, B; C, D) topology, without invoking admixture, the allele frequency difference between A and B (a − b) and between C and D (c − d) should be unrelated and hence results in f4 = [(a − b)(c − d)] = 0. Significantly positive f4 implies gene flow between A and C or B and D. Significantly negative f4 implies gene flow between A and D or B and C (Lopes et al., 2023). As used here, when D is an outgroup, it allows testing past introgression between A and either B or C (Lopes et al., 2023). A significantly positive f4 implies gene flow between A and C, and a significantly negative f4 implies gene flow between B and C (Patterson et al., 2012; Hu et al., 2024). The significance of the f4-statistics and f3-statistic is based on z-scores, significance positive (z > 3) and significance negative (z < −3) (Harris and DeGiorgio, 2017). Significant f4 and f3 values can also be interpreted as a rejection of a given topology (Peter, 2016; Zheng and Janke, 2018). Finally, we performed D-statistics analyses using qpDstat in the Admixtools v.7.0.2 package to detect gene flow events between the four populations in the form of D (Pop1, Pop2; Pop3, M. ochrogaster) (Patterson et al., 2012). We performed the following tests: (i) (ALS, DQ; LJ, M. ochrogaster); (ii) (ALS, DQ; JC, M. ochrogaster); (iii) (ALS, LJ, JC; M. ochrogaster); (iv) (DQ, LJ; ALS, M. ochrogaster); (v) (JC, LJ; ALS, M. ochrogaster). A significant positive value of D (z > 3) indicates high affinity between Pop1 and Pop3. If D is significantly negative (z < −3), it indicates high affinity between Pop2 and Pop3 (Zhang et al., 2020).

The deleteriousness of the missense mutations was also diagnosed using the Grantham Score (GS)—a measure of the physical/chemical properties of amino-acid changes (Grantham, 1974). Use the Python program “Grantham_score_calculator.py” to calculate the Grantham Score for each missense mutation site (Grantham, 1974). Additionally, count the sites with a score greater than 150, as these are defined as the genetic load of deleterious missense mutations (Li et al., 1984).

In general, positive selection gives rise to lower genetic diversity within populations and higher genetic differentiation between populations (Wu et al., 2014). The genetic differentiation index F_ST and the average proportion of pairwise mismatches over all compared sequences have been widely used to detect selection (Tajima, 1983; Weir and Cockerham, 1984). To detect selection signals possibly associated with local adaptation, we used a sliding-window method (10kb window, 1kb step) to calculate the genome-wide distribution of F_ST values and θ_π ratios for the two populations, implemented in vcftools v0.1.14 (Hu et al., 2020). We applied Z transformation for F_ST values and log2 transformation for θ_π ratios and considered the windows with the top 5% Z (F_ST) and log2 (θ_π ratio) values simultaneously as the candidate outliers under strong selection (Hu et al., 2020). All outlier windows were assigned to the corresponding snp and gene. We used Metascape to analyze the enrichment of selected genes located in specific regions (Zhou et al., 2019). Each significantly enriched category contained at least three genes, and a hypergeometric test was used to estimate significance (p < 0.05).

The SLAF-seq genomic data, with an average genome coverage of 10.21x, and the re-sequencing data, with an average genome coverage of 39.87x ( Supplementary Table S1 ), have been uploaded to NCBI. The SLAF-seq genomic data and re-sequencing data were mapped using BWA v.0.7.17-r1188. The raw SNPs underwent quality control using stringent filtering criteria, resulting in a total of 130,769,878 high-quality SNPs datasets for the re-sequencing data and 151,301 high-quality population SNPs datasets for the SLAF-seq genomes ( Supplementary Table S1 ).

ML tree constructed using SLAF-seq genomic data showed that E. miletus in ALS and JC basically cluster into separate clade, while E. miletus in DQ and LJ form sister clades. Among these four groups, ALS was the first to diverge and occupied a basal position, followed by JC, and the last two to diverge, DQ and LJ, were clustered into sister branches ( Figure 1A ). We then used the re-sequencing data to construct a ML and Bayesian tree ( Figures 1B, C ), and the phylogenetic tree constructed using the re-sequencing data also supported the existence of significant differences between the groups, and the results showed that the first to differentiate was still the ALS, followed by JC, and lastly DQ and LJ that were clustered into sister branches. Meanwhile, our results were also confirmed in the previous study on the genetic structure analysis from five regions in E. miletus, when k = 2, ALS was the first to be separated, and k = 3, JC were separated (Ren et al., 2023). And it was found that in the principal component analysis, using of PC1 and PC2 differentiated between ALS and JC, and it was the use of PC1 and PC3 that separated DQ and LJ (Ren et al., 2023). Therefore, based on the phylogenetic tree results, genetic structure and PCA analyses, we suggested that the ALS was the earliest to diverge from the common ancestor, whereas the divergence of the DQ and LJ was a notable recent event.

In order to assess the relative genetic similarity between ALS and other populations, we selected LJ, DQ, and JC as the reference populations, with M. ochrogaster serving as the outgroup. The f3 results indicated a significant positive value, with ALS exhibiting the highest f3 value in comparison to DQ, which suggests the greatest degree of shared genetic history between them. This was followed by ALS with JC, and finally ALS with LJ ( Figure 2A ). To explore the potential ancestral origin of the groups the admixture f3 (A, B; C) analyses were conducted, and the f3 were all significant positive values, indicating that there was no admixture event between the four regions ( Figure 2B ). The f4-statistics use four-taxa topologies in the form (A, B; C, D). Here, D represents M. ochrogaster. Taking turns in order, ALS, DQ, JC, and LJ serve as A, B, and C respectively. We found that significant gene flow was detected between all groups except between ALS and LJ ( Figure 2C ). Although these statistics cannot discern the direction of the introgressions nor distinguish between a scenario of continuous gene flow, the f4-statistics results strongly supported that ALS had been involved in past gene flow events in DQ and JC. In order to better understand the intergroup introgression signals, we performed D-statistics analysis using qpDstat in the Admixtools software package to test events of gene flow between the four groups. From the results of the D-statistics analysis, we found that the D values of all possible events detected when Pop1 was in ALS were significantly negative (z < −3), which suggested that the other three groups besides ALS had a high affinity. When Pop3 is ALS and Pop1 is DQ and JC, respectively, the D value showed a significant positive value, which indicated that ALS had a high affinity with DQ and JC ( Figure 2D ). The results of both f4-statistics and D-statistics analyses supported the existence of strong gene flow and high affinity between the ALS and both DQ and JC. E. miletus in Yunnan is mainly concentrated in the eastern part of the HDR, and its fossils had only been found in the Quaternary Pleistocene strata in China (Young, 1935). Some studies had shown that Eothenomys originated in the northern part of the subcontinental monsoon zone, and its ancestors entered the HDR from the northeastern part of the subcontinental monsoon zone to live mainly in the HDR with an altitude of 2,000–2,500 m sea level (Wang et al., 2022). Therefore, by integrating previous studies with data from our sampling sites, along with phylogenetic analyses and gene introgression between groups, as well as findings from earlier research on the genetic differentiation of E. miletus, we hypothesize that JC and DQ may have formed a group that continued to migrate to higher elevations after migrating to ALS from a common ancestor. But the question of whether the LJ was formed by the diffusion of ALS or evolved from an ancient ancestor needs to be followed up with a more detailed analysis.

We performed Grantham Score evaluation, we found that the maximum Grantham Score of all missense mutation sites in the four populations was 145, which was less than 150; there were no deleterious mutations in the missense mutation sites we detected ( Figure 3 ; Supplementary Table S2 ).

Considering that E. miletus lived in different geographical ranges and climatic environments and underwent a certain period of genetic differentiation, we focused mainly on the identification of selection and local adaptation between the two populations. Based on environmental differences, we selected the DQ population, which was relatively more rigorous compared to the other three populations, as our subject of study. Using F_ST and θ_π methods, between DQ and ALS, DQ and JC, as well as DQ and LJ, we selected 428, 452, and 429 genes in DQ, respectively ( Supplementary Table S3 ). The functional enrichment found that some genes were enriched in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of cancer pathway (ko05200, p < 0.01), renal cell carcinoma (ko05211, p < 0.01), choline metabolism in cancer (ko05231, p < 0.01), phospholipase D signaling pathway (ko04072, p < 0.01), chronic granulocyte leukemia (ko05220, p < 0.01), RAS signaling pathway (ko04014, p < 0.01), Rap1 signaling pathway (ko04015, p < 0.01), sphingolipid signaling pathway (ko04071, p < 0.01), Chagas disease (ko05142, p < 0.01), amoebiasis (ko05146, p < 0.01), toxoplasmosis (ko05146, p < 0.01), amoebiasis (ko05146, p < 0.01), and toxoplasmosis pathways (ko05145, p < 0.01) ( Figure 4 ; Supplementary Table S1 ) and the gene ontology (GO) term of olfactory placode formation (GO:0030910, p < 0.01), olfactory placode development (GO:0071698, p < 0.01), and olfactory placode morphogenesis pathways (GO:0071699, p < 0.01) ( Figure 5 ; Supplementary Table S1 ). We identified three selected genes associated with the olfactory placode in DQ: Six1, Six4, and Sox2.