Three landmark schemes were employed: A 24–landmark scheme covering the complete occlusal outline and major morphological structures. A 12–landmark scheme focused on the anteroconid complex. A 6–landmark scheme targeting the most variable points identified in previous analyses. Morphometric measurements (taken from the photographs with aid of TPSDig software) were based on van der Meulen and Zagwijn (1974) and Nadachowski (1982) and adapted to integrate the landmarking approach of Wallace (2006). The positions of the landmarks, the descriptive scheme, and the linear variables are shown in Figure 2.

Morphotypes were classified using both (i) BRA4 angle–based system (after Smirnov et al. 1986; Ponomarev and Puzachenko 2017) measured as the degree of tilt of the fourth buccal re-entrant angle on the first lower molar, and (ii) Nadachowski’s (1982) structural classification, which distinguishes 13 morphotypes (A–M) based on enamel field arrangement and developmental complexity. See Figure 3 for details.

Geometric morphometric analyses were conducted using MorphoJ (Klingenberg 2011) and PAST5.3 (Hammer et al. 2001). Goodall’s F-values (Goodall 1991) and Pillai’s trace (Pillai 1955) were applied for basic between-sample comparisons. After Procrustes superimposition (Rohlf and Slice 1990), we conducted Principal Component Analysis (PCA) (Jolliffe 2002) to explore shape variance. Canonical Variate Analysis (CVA) (Fisher 1936; Klingenberg 2011) to assess stratigraphic and geographic clustering. Procrustes analysis of variance (ANOVA; Goodall 1991; Klingenberg and Monteiro 2005) to test the effects of geography (locality) and stratigraphy (MIS) on molar shape. Discriminant Function Analysis (DFA) (Klecka 1980) to evaluate classification success between groups. Procrustes and Mahalanobis distances were used to quantify divergence between populations. Regression analyses explored trends in size and shape across stratigraphical and spatial units. Morphotype frequency dynamics were evaluated using H' and coefficient of variation (CV), with additional consideration of skewness and kurtosis in key traits (L, W, AC1, AC2, B1–B3). Standard statistics and between-samples comparisons (using Pearson correlation, Bray-Curtis distances, significance test of differences, etc.) were computed in Microsoft Excel and PAST5.3 software and are presented in full in File S2.

Comparing overall variation in mean values of standard linear measurements (21 variables) across all populations, we found no significant differences between samples (ANOVA, F (df = 10) = 0.029; p = 1; Kruskal–Wallis test: H = 0.774; p = 0.999). All samples exhibited significant normal distributions (Shapiro–Wilk W = 0.808–0.835; p = 0.0035–0.0011) and high similarity in metric profile as indicated by high between-sample correlations (r: 0.991–0.999), high Bray–Curtis similarity values (0.93–0.97), and low Mahalanobis distances (0.581–1.271). These results suggest that all samples represent a single phenotypic unit consistent with the taxonomic identity of the mid-European Stenocranius anglicus.

However, the detailed analysis revealed notable variations in both metric and proportional characteristics across different stratigraphic and geographic contexts. In general, specimens from MIS 3 and older showed larger molar dimensions than those from Late Vistulian/Early Holocene localities. The largest molars were found in Bojnice (MIS 4) and Bišilu (MIS 2), with maximum recorded lengths of 3.900 mm and 3.088 mm, respectively. In contrast, the smallest were from Maštalná (MIS 1) and Stránska skála SSJ (MIS 12/14), with minimum lengths of 2.069 mm and 2.162 mm, respectively. The width of the molars followed a similar pattern, with the widest tooth recorded in Bojnice (1.333 mm) and the narrowest in Maštalná (0.660 mm). Among proportional characteristics, the anteroconid complex showed substantial variation, with the highest relative values observed in Zkamenělý Zámek and Bojnice, whereas the lowest values appeared in Holštejnská and Srnčí. The most extreme AC1/L ratio (relative anteroconid length) was observed in Stránska skála SSJ (0.562) and Bišilu (0.574), while the lowest values were recorded in Tučín (0.503) and Barová (0.502). Similarly, the AC2/AC1 ratio varied, with the highest values in Maštalná (0.832) and the lowest in Zkamenělý Zámek (0.775). The relative width B1/W was the largest in Bišilu (1.148) and the smallest in Zkamenělý Zámek (0.944). The variability in shape proportions, as measured by coefficients of variation, was higher in some Late Vistulian/Early Holocene samples, particularly in features of the anteroconid complex (AC1 and AC2). Detailed survey of biometric data is available in Files S1–S3.

Geographically, specimens from Moravian localities (e.g., Stránska skála SSJ, Balcarka, Zkamenělý Zámek) exhibited slight but consistent differences from Slovak samples (e.g., Šarkanica, Bojnice, Muráň 3), particularly in the width and proportions of the anteroconid and basin structures. The highest coefficients of variation (CV) were found in Šarkanica, especially for AC2 and AC1/L, suggesting higher morphological diversity, while the lowest variability was observed in Holštejnská and Barová. Extreme skewness and kurtosis values for specific characteristics, such as the B2 and B3 proportions in Šarkanica and Tučín, the sites with the highest dominance of S. anglicus, suggest strong tendencies toward rearrangements of specific phenotypic traits.

First, we compared the strength of three different landmark schemes against the geographical and stratigraphical positions of particular samples using ANOVA and Canonical Variates Analysis (CVA). Both analyses revealed that the geographical position of the localities, as well as their stratigraphy (recorded by the marine isotope stage (MIS) of the respective locality/stratigraphic layers), significantly affect shape variation across all landmark schemes.

Procrustes ANOVA for the 24–landmark scheme showed a strong effect of locality on shape (F = 28.28, p < 0.0001; Pillai’s trace = 2.51) and MIS (F = 51.12, p < 0.0001; Pillai’s trace = 1.67). When the number of landmarks was reduced to 12, Goodall’s F-values slightly increased for both locality (F = 31.63) and MIS (F = 57.86), but Pillai’s trace values decreased (1.86 for locality; 1.33 for MIS), indicating a reduction in the proportion of variance explained. The 6–landmark scheme showed the weakest effects, with locality (F = 29.19; Pillai’s trace = 0.85) and MIS (F = 56.47; Pillai’s trace = 0.63), confirming that shape differentiation remains significant but is notably weaker with fewer landmarks (Table 2). Overall, the effect of geographical position consistently shows more substantial shape variation effects than stratigraphy (MIS) across all landmark schemes. Locality has consistently higher Pillai’s trace values than MIS in the 24–landmark scheme (2.51 vs. 1.67), meaning that geographic differences explain more of the shape variation than stratigraphic differences. Goodall’s F-values are higher for MIS than for locality in the 12– and 6–landmark schemes, but this likely reflects statistical sensitivity rather than a true biological pattern (Table 2).

Mahalanobis and Procrustes distances are generally larger for geographical comparisons, suggesting that shape differences among localities are greater than those among stratigraphic stages. Geographical position could therefore be a more dominant factor influencing shape variation. However, MIS still has a significant effect, particularly in the full landmark dataset, although it explains less overall shape variation. This pattern holds across all landmark schemes, but when reducing the number of landmarks, the stratigraphic effect appears less affected than the locality effect. Pairwise Mahalanobis distances indicated that shape differentiation among localities and MIS groups was strongest in the 24–landmark scheme. For locality, the highest Mahalanobis distance was observed between Stránska skála SSJ and Barová (5.82), while the lowest (1.74) was between Maštalná and Bojnice. In the 12–landmark scheme, the highest distance decreased to 5.14, and in the 6–landmark scheme, it dropped to 2.34, indicating a progressive reduction in shape distinctiveness. For stratigraphy, the 24–landmark scheme showed the highest Mahalanobis distance spanned between MIS 12 and MIS 4 (5.10), i.e., Stránska skála SSJ and Bojnice, while the 6–landmark scheme reduced the span to 2.93. The smallest distances were also reduced, suggesting that fewer landmarks lead to less separation among stratigraphic groups. Procrustes distances followed a similar trend. While group separation remained significant in all schemes (all p < 0.0001), Procrustes distances were consistently lower in the 6–landmark scheme, meaning the amount of shape variation captured was reduced (Table 3).

Overall, the 24–landmark scheme provides the most comprehensive representation of shape variation, while the 12–landmark scheme retains significant, though slightly reduced, explanatory power. The 6–landmark scheme still detects locality and MIS effects but at a much weaker level, suggesting that landmark reduction may lead to a loss of biologically meaningful shape variation. These findings indicate that landmark reduction should be carefully considered, particularly when studying complex morphological differences across space and time.

The correspondence between stratigraphic classification based on CVA (Canonical Variates Analysis) and the actual stratigraphic position of particular samples revealed varying levels of classification accuracy. The highest classification success was observed in MIS 6 (Tučín), with 88% of specimens correctly classified, indicating distinct molar morphology during that glacial period. Similarly, the Last Glacial Maximum (MIS 2, Šarkanica) achieved 80% accuracy, supporting the notion that glacial periods produced morphologically cohesive vole populations. MIS 12 (Stránska skála SSJ) showed moderate accuracy (58%), consistent with its early Middle Pleistocene position and potential ancestral morphology. In contrast, Late Vistulian/Early Holocene MIS 1 (multiple localities) had 69.5% accuracy, with notable misclassifications into pre–LGM periods, suggesting morphological convergence or stabilization. MIS 3 (Balcarka, Zkamenělý Zámek) had 71% accuracy, while MIS 4 (Bojnice) was lower at 38.6%, with substantial misclassification into MIS 3 and MIS 1, implying morphological overlap during this period. Regarding the centroid positions of particular MIS units, those of the Middle Pleistocene age (i.e., MIS 12 and MIS 6) show considerable distance from the Late Pleistocene and Holocene samples. At the same time, MIS 4 and MIS 2 (both extreme pleniglacials) appeared close to each other, similarly to MIS 3 and MIS 1 (Fig. 4).

Across localities, the classification accuracy from the CVA varied notably. Tučín displayed the highest correct classification rate at 82%, suggesting that its narrow-headed vole molars possessed distinct morphological traits, likely shaped by the MIS 6 glacial environment. Other localities with relatively high accuracy included Muráň 3 (68.1%), Býčí Skála (66.5%), Šarkanica (64.2%), and Balcarka (63.2%), all of which indicate strong morphological differentiation. In contrast, Bojnice had the lowest accuracy (24.1%) and was frequently misclassified as Šarkanica, suggesting morphological convergence or shared traits. Similarly, Maštalná, Zkamenělý Zámek, and Srnčí showed moderate to low accuracy (around 41%), with misclassifications spread across several groups, indicating either internal variability or morphological overlap with other localities. Worth mentioning is a resemblance of the sites Barová (57.3%), Dzeravá (50.0%), and Býčí Skála, which displayed moderate classification success, corresponding to their similarities in other phenotype comparisons, and the distant position of Stránska skála SSJ (with 58% of correct classification) and Tučín, both representing a Middle Pleistocene context (compare centroid positions in Fig. 5).

Both classification schemes demonstrated the highest accuracy for cold-adapted populations during glacial stages- Tučín (MIS 6) and Šarkanica (MIS 2) stood out in both analyses with high classification rates (~82% and ~80%, respectively). Localities tied to transitional periods, such as Bojnice (MIS 4), showed poor classification accuracy in both schemes (24% locality-based, 39% stratigraphic), underscoring morphological variability or convergence during these times. Holocene localities had moderate success in both systems (~50–57%), suggesting some morphological stabilization before extinction. Overall, both classification matrices show consistent patterns that highlight the narrow-headed vole’s evolutionary responses to climatic fluctuations. A relatively low number of correctly classified specimens (locality ~55%; stratigraphy ~70%) indicates that the phenotypic characteristics of the European narrow-headed vole are stable but exhibit plasticity, with a high potential for phenotypic rearrangements during periods of unfavorable environmental conditions (e.g., Holocene localities).

Morphotype scheme I – BRA4 angle classification (Ponomarev and Puzachenko 2017)

During the Middle to Late Pleistocene (pre–LGM), diversity values ranged widely across both Czech and Slovak localities (Fig. 6). In the oldest locality, Stránska skála SSJ, the diversity index was moderate (H' = 1.637). This pattern continued into Tučín (H' = 1.585) and Bojnice (H' = 1.639), suggesting relatively stable diversity during early glacial cycles. A notable increase in diversity occurred during MIS 3, particularly at Zkamenělý Zámek (Czech Republic), which recorded the highest pre–LGM diversity (H' = 1.980) across all localities, exceeding values in all other samples. Similarly, Balcarka showed high diversity (H' = 1.787). Contrasting patterns of diversity across Slovakia and the Czech Republic were evident during the Last Glacial Maximum (LGM). In Šarkanica (Slovakia), an apparent decline in diversity was observed (H' = 1.520). While the Muráň 3 sequence (Slovakia), Muráň 3/6 (24.4ky BP) exhibited one of the highest diversity values recorded during the LGM (H' = 1.923), Muráň 3/4 (16.3ky BP) and Muráň 3/3 declined to H' = 1.230–1.332, foreshadowing population contractions toward the beginning of the Holocene. In the Czech Republic, diversity during the LGM suggests high morphological variability. Following the LGM, a general increase in diversity was evident across both countries, with most post–glacial localities exhibiting H' values above 1.6. The highest diversity post–LGM was recorded in Maštalná/13 (H' = 1.955) and Barová 14–10ky (H' = 1.979), while Dzeravá/3 (Slovakia) reached H' = 1.889, showing high variability across regions. In Holštejnská (Czech Republic), layer 5 (H' = 1.759) and layer 6 (H' = 1.864) showed elevated diversity, supporting a trend of increasing variability toward the Holocene. In Slovakia, the Dzeravá sequence showed stable high diversity (H' = 1.796–1.889), consistent with long-term population stability or refugial continuity. Similarly, Býčí (Czech Republic) showed high diversity in layers dated to ~11.2–12.1ky BP (H' = 1.686–1.845), aligning with early Holocene warming. Overall, Czech localities (e.g., Balcarka, Barová) tended to exhibit slightly higher and more variable diversity compared to Slovak sites, which, while also diverse, showed more consistent values (Dzeravá and Muráň 3 sequences).

Morphotype scheme II – morphotype classification after Nadachowski (1982)

Morphological diversity assessed using Nadachowski’s (1982) classification proposal revealed trends that broadly align with those derived from the BRA4 angle-based classification, still with notable local and chronological distinctions (Fig. 7). Stránska skála SSJ (MIS 12) exhibited lower diversity (H' = 1.501), consistent with a more primitive morphotype spectrum and less morphological variation. Diversity then peaked during the MIS 5–3 interval, with the highest values at Bojnice (H' = 1.93) and Zkamenělý Zámek (H' = 1.843). Tučín (MIS 6) also showed relatively high diversity (H' = 1.616), possibly reflecting morphotype turnover during the pre–LGM cold period. During the LGM (MIS 2), diversity remained moderate to high, as seen at Šarkanica (H' = 1.527), suggesting the continued presence of multiple morphotypes and potentially adaptive responses to glacial conditions. However, some late LGM localities, such as Muráň 3/6 (H' = 1.484) and Holštejnská/4 (H' = 0.9503), began to exhibit decreasing diversity, indicating incipient morphological homogenization preceding the species’ decline. In the post–LGM (Late Vistulian/Early Holocene) period, diversity patterns varied geographically: Býčí Skála and Dzeravá localities generally retained moderate to high diversity (e.g., Býčí/8a: H' = 1.831, Dzeravá/7: H' = 1.713), suggesting refugial stability or persistent morphological variability until the species’ final disappearance. Srnčí and Maštalná populations showed lower and fluctuating diversity (e.g., Srnčí/8: H' = 0.6365, Maštalná/10: H' = 0.4506), likely reflecting small, isolated populations with restricted morphotype spectra at the edge of extinction. Notably, Barová layers and Maštalná/14–15 exhibited high Shannon diversity (H' > 1.95), comparable to peak glacial values. In contrast, Holštejnská/4 and Maštalná/10 had the lowest diversity values (H' = 0.9503 and 0.4506), underscoring their role as terminal occurrences of the species in the region.

Overall patterns of these two approaches align (Fig. 8). Both classification systems identify high diversity during MIS 3–LGM and a post–LGM decline, but Nadachowski’s scheme (II) highlights fine-scale morphotype turnover and rare morphotype occurrences (J–M) in late periods.

Geographic differences are more pronounced in Slovak localities, consistently showing higher diversity than the Bohemian and Moravian sites. It is markedly illustrated by a comparison of morphotype diversity (scheme II) in Czech Republic, Slovakia and Poland arranged in stratigraphical units (Fig. 9) proposed by Nadachowski (1982). It demonstrates similarity in the phenotype diversity trends among the Czech Republic (mostly Moravia) and Poland during the glacial stages, contrasting to the specific situation in the Slovak Carpathians (characterized by significantly higher diversity during MIS 4 – MIS 2). The characteristics of the Holocene samples, both in Slovakia and the Czech Republic, are situated quite apart from the cluster of glacial samples. Worth mentioning are also the diversity values respective to all studied samples (Czech and Slovak sites together – Fig. 9A), which (particularly in MIS 4–2 samples) considerably exceed values of individual regions (Fig. 9B), indicating thus a greatly pronounced effect of beta diversity in phenotype pattern of the species.

Phenotype variation along the dataset means

To assess the degree of morphological divergence of the first lower molars across localities, each locality’s mean shape was compared to the dataset mean using Procrustes and Mahalanobis distances. Significant shape differences were observed in most localities, as confirmed by permutation tests (1000 runs; p < 0.001 unless noted otherwise; Fig. 10).

Regarding the Procrustes distances, the highest shape divergence was observed at Stránska skála SSJ (0.0510, p < 0.0001), likely representing a primitive Middle Pleistocene (MIS 12) morphotypes. Other notable deviations were found at Tučín (0.0369) and Muráň 3 (0.0354), indicating distinctive morphologies associated with pre–LGM and glacial periods. Šarkanica (0.0327) also exhibited significant divergence, potentially reflecting cold-adapted forms. In contrast, Holocene sites such as Maštalná (0.0133) and Srnčí (0.0156) showed the lowest divergence, suggesting shape convergence toward the dataset mean in post–LGM populations. Srnčí was the only site with non-significant shape differentiation (Procrustes p = 0.236), indicating high similarity to the global mean, however, the number of specimens from this locality was very low (22 individuals), which explains the statistical insignificance. Regionally, Slovak sites (Muráň 3, Šarkanica, Bojnice, Dzeravá skala) tended to show greater shape divergence than Bohemian and Moravian sites, except for Maštalná and Srnčí, which were closest to the mean (Fig. 10). The Mahalanobis distances indicate clear morphological differentiation among the examined localities, with all pairwise comparisons being highly significant (p < 0.0001), except, again, for a case of Srnčí, which shows relatively low differentiation from Bišilu (p = 0.0461), Holštejnská (p = 0.0383), and Maštalná (p = 0.1575). Stránska skála SSJ shows the highest Mahalanobis distances across comparisons, reinforcing its distinctiveness and alignment with early Pleistocene morphotypes. Conversely, the smallest distances are observed among Late Vistulian/Early Holocene sites such as Maštalná, Srnčí, and Zkamenělý Zámek, indicating greater morphological homogeneity. The distinctiveness of Tučín and Šarkanica further highlights morphological shifts associated with the extreme glacial conditions and eudominant position of S. anglicus in mammal communities (84% and 62%).

Between-population phenotype variation: common patterns and divergences

As demonstrated above, a particular population might differ significantly in both metric and non-metric traits. While there is extensive overlap in metric variables, the proportions, the frequency of specific morphotypes, morphotype diversity, and ratios of metric variables exhibit clearly pronounced between-population differences, particularly when confronted with the actual dominance of the species in small ground-mammal communities (Figs 6, 7). Despite specific differences between the two morphotype sets, the overall pattern of morphotype dynamics is the same: as population abundance declines, morphotype diversity decreases dramatically (Figs 8, 9). The other phenomenon worth mentioning here is a significant between-population variation in the sites representing continuous records of local history – Dzeravá skala (53–20 ky) and Maštalná (16–8 ky). The extensive range of variation recorded in these two sites encompasses the states of the respective variables in most other populations. This is also demonstrated by the results of PCA (Figs 6B, 7B), where a distinct outlier position shows either (i) the populations from a deeper Pleistocene past (SSJ, Tučín); (ii) peak glacial stages in the Carpathians (MIS 4 Bojnice, MIS 2 Šarkanica); and (iii) the populations of a terminal stage of the species extinction (Srnčí/4, Holštejnská/4, Maštalná/10, Muráň 3/3–4).

A comparison of metric variables (Figs 11, 12) revealed an unexpected pattern of between-site divergence. Both in mean differences in metric variation from overall mean values (Fig. 11) and PCA first component of all linear variables (Fig. 12), the total set of samples splits into two distinct clusters, separating closely related samples of Býčí Skála, Barová, and Dzeravá skala from all the remaining samples, as revealed also by a comparison of non-linear variables (compare Figs 5 and 10).

The PCA results of a high-resolution biometric record (155 linear distances among individual landmarks) from particular sites (all layers included), demonstrated in Figure 13, split the set of sites into two distant clusters: A – sites of the Carpathian region and those from deeper Quaternary past; C – sites from the Bohemian Massif (Bohemia + Moravia) with a dense cluster B of Býčí, Barová and Dzeravá situated in between them. The respective PCA results indicate that, for populations in cluster C, the essential source of variation is PC1 (roughly size characteristics); in cluster A, the role of additional variation (shape of the anteroconid complex, proportions) is distinctly more pronounced. In contrast, cluster B seems to show a balanced response to both drivers.

Even more distinct differences between the B cluster sites and all other sites appear in the proportion ratios of metric variables (Figs 14, 15), which exhibit values larger than the overall mean of the respective ratio variables and show consistent relations with variation in metric variables (R² = 0.680).

Worth mentioning is that the results of our study concerning the phenotypic divergences among the clades from Czech Massif and Carpathians and MIS 4 sample are robustly supported also by results of aDNA analyses (Fig. 16). The early haplotypes (EF) present in the Carpathian MIS 4 sites (Bojnice, Peskö/12) were replaced during MIS 3 by modern haplotypes EB (in the Carpathians) and even EA distributed mainly in Bohemian Massif populations.