We examined 272 Blarinellini samples, including 127 individuals of Blarinella quadraticauda, 126 of B. wardi, 14 of Parablarinella griselda, and 5 of P. latimaxillata, representing both recognized genera and all four currently recognized species. Samples were obtained from 95 geographic localities, spanning most of the tribe’s known distribution (Fig. 1). All tissues (muscle or liver) were stored at −80 °C. Voucher specimens and associated tissues are deposited in national and provincial zoological collections. Animal sampling and handling followed the guidelines of the American Society of Mammalogists (Sikes et al. 2016) and were approved by the Animal Care and Ethics Committee of Anhui Normal University (Approval No. AHNU-ET2024005).

Total genomic DNA was extracted using a DNeasy blood and tissue kit (Tiangen, China). We amplified and sequenced two mitochondrial markers (cyt b, 1140 bp; 16S rRNA, ~524 bp) and three nuclear loci (ApoB, ~516 bp; BRCA1, ~770 bp; RAG2, ~750 bp) by PCR, using the primers and annealing temperatures described in Chen et al. (2023). Sequence availability varied among loci because of differences in amplification success and DNA quality. A total of 271 cyt b sequences were obtained from the 272 examined individuals, with one B. quadraticauda individual (voucher number ZIN96272; sample ID BQ95) lacking cyt b sequence data (Table SS1). In addition, we downloaded 179 gene sequences from GenBank, including two Blarina brevicauda and one Cryptotis parvus outgroup sequences (Table SS1). The 654 newly generated sequences were combined with downloaded sequences and aligned in MUSCLE (Edgar 2004), followed by manual inspection and correction in MEGA 11 (Tamura et al. 2021). Sample information and GenBank accession numbers are provided in Table SS1.

We inferred phylogenetic relationships within Blarinellini using Bayesian inference (BI) and maximum likelihood (ML) approaches. BI analyses were conducted in MrBayes v3.2 as implemented in PhyloSuite v1.2.2 (Zhang et al. 2020), while ML analyses were performed using RAxML via raxmlGUI 2.0 (Edler et al. 2020).

The phylogenetic reconstructions presented in the main text were based on a concatenated dataset comprising all five markers (cyt b, 16S rRNA, ApoB, BRCA1, and RAG2). This dataset was used to infer the primary BI and ML phylogenies and served as the basis for subsequent divergence time estimation. Optimal partitioning schemes and nucleotide substitution models were selected using PartitionFinder v2.0 (Lanfear et al. 2017) under the Bayesian Information Criterion (BIC) with a greedy search algorithm. Data blocks were defined by gene and, for protein-coding loci, by codon position, with branch lengths linked across partitions. The resulting partitioning scheme and substitution models were applied in all subsequent BI and ML analyses.

For BI analyses, two independent runs, each comprising four Markov chains, were conducted for 10 million generations, with trees and parameters sampled every 10,000 generations. Convergence was assessed using the average standard deviation of split frequencies and parameter traces, and the first 25% of samples were discarded as burn-in. ML analyses were conducted under the rapid bootstrap algorithm implemented in RAxML, with 1000 bootstrap replicates to assess node support.

To assess the consistency of phylogenetic relationships across marker types, we additionally reconstructed phylogenies from concatenated mitochondrial (cyt b + 16S rRNA) and nuclear (ApoB + BRCA1 + RAG2) datasets —using the same analytical framework. Single-gene trees were also inferred for each locus. These supplementary phylogenies are provided in the Figures S1–S8.

We estimated divergence times for major nodes within Blarinellini in BEAST v2.7.7 (Bouckaert et al. 2014) using the concatenated mitochondrial and nuclear dataset. The partitioning scheme and substitution models listed in Table SS2 were applied, with substitution models unlinked among partitions and a shared tree topology across all partitions. An uncorrelated lognormal relaxed molecular clock was applied to all partitions. Because the dataset included extensive intraspecific sampling and shallow population-level divergences, a coalescent constant-population prior was employed.

We implemented two fossil-based calibration points using lognormal prior distributions. The first constrained the divergence between Blarinellini and its sister tribe, Blarinini. This calibration was based on the earliest fossil occurrence of Blarinini from the Barstovian of North America (~16.3–13.6 million years [Ma]; Repenning 1967) together with early Miocene records of Blarinellini from Europe and North America (Harris 1998; Rzebik-Kowalska 1998). The lognormal prior parameters were set to offset = 16.3, mean = 0.8, and standard deviation = 0.75, corresponding to a 95% prior probability interval of 26.0–16.8 Ma. The second calibration was based on the earliest known fossil record of Cryptotis (~9 Ma; Harris 1998). The lognormal prior parameters were specified as offset = 9.0, mean = 3.5, and standard deviation = 1.78, yielding a 95% prior probability interval of 16.0–9.0 Ma.

Markov chain Monte Carlo (MCMC) analyses were run for 100 million generations, sampling every 20,000 generations. We conducted two independent runs and discarded the first 10% of samples from each as burn-in. Convergence and mixing were evaluated in Tracer v1.7.2 (Rambaut et al. 2018), and all parameters achieved effective sample sizes (ESS) > 200. A maximum clade credibility (MCC) tree was generated in TreeAnnotator v2.7.7 using median node heights.

Because sampling localities for P. griselda and P. latimaxillata were limited, we conducted population genetic analyses only for B. quadraticauda and B. wardi. Due to missing cyt b data for one B. quadraticauda individual, cyt b-based analyses included 126 individuals of B. quadraticauda. Analyses were performed separately for the mitochondrial cyt b gene and the nuclear ApoB gene, which provided greater sample coverage than the remaining nuclear loci.

For each marker, haplotypes were identified and standard diversity indices, including the number of haplotypes (H), haplotype diversity (Hd), nucleotide diversity (π), number of segregating sites (S), and average number of nucleotide differences (k), were calculated in DnaSP v5.10 (Librado and Rozas 2009). Missing data were handled using pairwise deletion, and nucleotide diversity was calculated based on the available sites within the alignment. Population genetic structure was evaluated using molecular analysis of variance (AMOVA) and pairwise F_ST in Arlequin v3.5.2 (Excoffier et al. 2010). Historical demographic patterns were assessed using mismatch distribution analyses in DnaSP and neutrality tests (Tajima’s D and Fu’s Fs; Tajima 1989; Fu 1997) in Arlequin. Statistical significance was evaluated using 1000 permutations where applicable. Missing data were treated by pairwise deletion in both DnaSP and Arlequin. Sample sizes for each subclade are provided in Table SS1.

Haplotype networks for cyt b and ApoB were constructed using PopART v1.7 (Leigh and Bryant 2015) under the median-joining algorithm (Bandelt et al. 1999). Alignment gaps and ambiguous nucleotide sites were treated as missing data and excluded from network inference. Genetic distances (p distance) among lineages were calculated in MEGA 11 (Tamura et al. 2021).

To evaluate whether the distributions of Blarinellini are consistent with long-term climatic niche stability, ecological niche models (ENMs) were constructed for B. quadraticauda and B. wardi using MaxEnt. Potential distributions were reconstructed for the Last Interglacial (LIG; 120–140 ka), the Last Glacial Maximum (LGM; ~22 ka), and the present. Occurrence records were compiled separately for each species from verified localities. To reduce spatial sampling bias and autocorrelation, duplicate records were removed and occurrences were spatially rarefied by retaining a single record within each 5 × 5 km grid cell. After filtering, 121 occurrence records for B. quadraticauda and 122 occurrence records for B. wardi were retained.

The background extent for model calibration was defined to approximate the accessible area (M) of Blarinellini, encompassing the known fossil and extant distribution range of the group, with the Yellow River as the northern boundary, the Indochinese Peninsula as the southern boundary, the Brahmaputra River as the western boundary, and the East China Sea as the eastern boundary. Nineteen bioclimatic variables were obtained from WorldClim v2.1 (Fick and Hijmans 2017) for current climatic conditions at a resolution of 2.5 arc minutes. Corresponding LGM and LIG climatic layers were obtained at the same or higher resolution and resampled to 2.5 arc minutes prior to analysis.

To reduce multicollinearity among predictors, variable selection was conducted in three steps. First, an initial MaxEnt model including all 19 variables was used to assess preliminary variable importance. Second, highly correlated variables (|r| ≥ 0.9) were identified using Pearson correlation analysis and removed. Third, variance inflation factor (VIF) analysis was used to exclude variables with VIF values > 10. Variables with consistently low contributions in the preliminary model were subsequently excluded. The final predictor set for B. quadraticauda included bio2, bio4, bio10, bio13, bio14, bio18, and bio19, whereas that for B. wardi included bio2, bio3, bio5, bio14, bio15, bio18, and bio19.

Model complexity was optimized separately for each species by evaluating combinations of feature classes and regularization multipliers. Optimal settings were selected using the lowest Akaike Information Criterion corrected for small sample size (AICc) and subsequently used for projections under present, LGM, and LIG climatic conditions. Final models were run with 10 replicates, using 75% of occurrence records for training and 25% for testing. Model performance was evaluated using the area under the receiver operating characteristic curve (AUC). All analyses were conducted in R v4.3.3 (R Core Team 2024) with MaxEnt v3.4.3, and optimal model settings are provided in Table S3.

To assess distributional stability through time, habitat suitability predictions for the LIG, LGM, and present were converted to binary maps using a threshold of 0.5. Suitable areas (≥ 0.5) were overlaid across time periods to identify stable and unstable habitats. Cells predicted as suitable in all three periods were classified as long-term stable habitats (“Throughout”), those suitable in two periods as relatively stable habitats, those suitable in one period as unstable habitats, and those unsuitable in all periods as “Never suitable”. Following previous studies, resistance values of 1, 10, 100, and 1000 were assigned to these categories, respectively. Stability analyses were performed in ArcGIS v10.2. Variable contributions for each species are provided in Table S4.

Both ML and BI analyses of the combined datasets yielded highly congruent topologies, supporting the reciprocal monophyly of Blarinella and Parablarinella (posterior probability [PP]/ bootstrap values [BS] = 1.0/100; Figs 2, S1–S3) and their respective sister-species pairs. This interspecific framework was fully recovered by individual mitochondrial (Figs S4–S5) and ApoB gene trees (Fig. S6). In contrast, single-locus nuclear markers consistently lacked intraspecific resolution, and individual RAG2 (Fig. S7) and BRCA1 (Fig. S8) loci failed to support the monophyly of B. quadraticauda and P. griselda, likely due to slower evolutionary rates.

Within B. quadraticauda, five geographically structured lineages (Subclades I–V) were identified. The basal lineage (Subclade V) occurred in the eastern Yunnan Plateau and northern Vietnam, whereas the remaining lineages were distributed across the Yunnan–Guizhou Plateau, the margins of the Sichuan Basin, the Qinling Mountains, and the northern Yunnan–southwestern Sichuan highlands. Within B. wardi, four geographically structured lineages (Subclades A–D) were recovered. Subclade A occurred east of the Lancang River, Subclade B occupied the upper Nujiang drainage west of the river, Subclade C was distributed between the Lancang and Nujiang rivers, and Subclade D occurred in the lower Nujiang drainage.

Molecular dating estimated the stem and crown ages of Blarinellini at ~18.2 Ma (95% HPD: 21.6–16.5 Ma; Fig. 3) and ~14.0 Ma (95% HPD: 18.6–9.8 Ma), respectively. Interspecific divergences were dated to ~5.3 Ma (95% HPD: 7.6–3.5 Ma) for B. quadraticauda–B. wardi and ~3.0 Ma (95% HPD: 4.6–1.8 Ma) for P. latimaxillata–P. griselda. Following initial intraspecific splits at ~2.3 Ma in both species, internal lineages diversified rapidly through the Pleistocene (Fig. 3). Blarinella quadraticauda exhibited a sequential divergence of Subclades II, I, and III between ~2.3 and 0.9 Ma, whereas B. wardi split into Subclades D, C, and the sister pair A–B from ~2.3 down to 1.0 Ma, with terminal node radiation for both taxa concluding in the mid- to late-Pleistocene. The divergence-time analysis recovered the same major subclades as the phylogenetic analyses, although the branching order of the major B. quadraticauda subclades differed slightly between the two trees.

Analysis of 271 cyt b sequences identified 109 haplotypes (H1–H109; Table SS1). The median-joining network recovered four well-separated haplotype groups corresponding to B. quadraticauda, B. wardi, P. griselda, and P. latimaxillata (Fig. 4). No haplotypes were shared among species.

Within B. quadraticauda, five distinct mitochondrial haplogroups corresponding to Subclades I–V were recovered. These haplogroups were separated by multiple mutational steps and numerous inferred intermediate haplotypes, indicating substantial genetic differentiation among geographic regions. Similarly, four well-defined mitochondrial haplogroups corresponding to Subclades A–D were identified within B. wardi. The distribution of haplotypes within both species showed strong geographic structure, with no haplotypes shared among major subclades.

Parablarinella exhibited markedly lower levels of mitochondrial variation. Parablarinella griselda comprised a compact cluster of closely related haplotypes, whereas P. latimaxillata contained only two closely related haplotypes. Despite their low internal diversity, both species remained strongly differentiated from Blarinella.

The nuclear ApoB dataset exhibited markedly reduced variation compared to cyt b and produced a shallow haplotype network (Fig. S9). Two primary haplotype groups corresponding to B. quadraticauda and B. wardi were identified, but several haplotypes were shared among mitochondrial subclades. Although overall geographic structure was weaker than in the mitochondrial dataset, a number of haplotypes were restricted to particular regions, especially within B. quadraticauda Subclades II and IV and B. wardi Subclade A.

AMOVA based on cyt b indicated pronounced population subdivision in both widespread species (Table 1). Most genetic variation was partitioned among subclades, accounting for 66.77% of total variation in B. quadraticauda (F_ST = 0.6677, p < 0.001) and 87.52% in B. wardi (F_ST = 0.8753, p < 0.001). Pairwise F_ST values and genetic distances further demonstrated substantial divergence among mitochondrial lineages (Fig. 5), whereas differentiation inferred from ApoB was considerably weaker (Table S5–S6).

Haplotype diversity was high in most subclades (Table 2). Significant negative Tajima’s D and Fu’s Fs, together with unimodal mismatch distributions, were observed only in B. quadraticauda Subclade III and B. wardi Subclade A (Fig. 6), indicating recent demographic expansion. Other subclades showed nonsignificant or inconsistent signals. ApoB diversity was much lower, with significant negative neutrality values limited to B. quadraticauda Subclades III and IV (Fig. S10; Table S5).

Ecological niche models showed high predictive performance for both species (B. quadraticauda: AUC = 0.919–0.956; B. wardi: AUC = 0.987–0.989), and the predicted present-day suitable areas were generally consistent with the observed occurrence records (Fig. 7).

The analysis predicted extensive suitable habitat for B. quadraticauda across southwestern and central China under all climatic scenarios (Fig. 7A–D). Suitable habitat was most extensive during the LGM and became more fragmented during the LIG. Areas predicted to remain suitable throughout all three periods were concentrated in the Hengduan Mountains, eastern Tibetan Plateau, and adjacent montane regions, indicating long-term climatic stability across these areas.

In contrast, suitable habitat for B. wardi was consistently more restricted and concentrated within the southern Hengduan Mountains and southeastern Tibetan Plateau (Fig. 7E–H). Although moderate shifts in distribution were observed among climatic periods, areas of long-term stability were largely confined to the southern Hengduan Mountains. Stable habitats across all three climatic periods were identified within the Hengduan Mountains for both species.