This study is based on 339 adult skulls belonging to the Japanese badger (M. anakuma) and the Asian badger (M. leucurus sensu lato; Table 1; File S1; Table S4).

The additional material included 536 adult skulls of the European badger (M. meles) and the Southwest Asian badger (M. canescens). Age classes were defined by scoring morphological skull features, such as the development of crests, the obliteration of sutures and tooth wear (Hysing-Dahl 1959; Wiig 1986).

Most of the material was collected in the 20^th century. A small part of our sample was collected in the 19^th century and at the beginning of the 21^st century. We examined the skulls of M. anakuma and M. leucurus deposited in the collections of 16 museums and institutions throughout the world (File S1). A significant number of M. leucurus skulls from the eastern part of species range was collected by the second author, V. Yudin. Additional material on the European and Southwest Asian badgers was described in details earlier (Abramov and Puzachenko 2013). Specimens of M. leucurus were provisionally divided into the geographical samples associated with the two subspecies – M. m. leucurus and M. l. amurensis. This expert split (a priori classification) was used as an external hypothesis for comparison with the results of sample clustering based on morphometric data.

Thirty measurements were made using a digital sliding calliper to the nearest 0.1 mm. A list and a scheme of cranial measurements are presented in Figure 1.

Figure 1. 

Measurements taken of badger skulls according to Abramov and Puzachenko (2005): p1 – condylobasal length, p2 – neurocranium length, p3 – viscerocranium length, p5 – palatal length, p6 – maxillary tooth-row length, p7 – length of upper carnassial tooth Pm4, p8 – greatest length between oral border of the auditory bulla and aboral border of the occipital condyle, p9 – length of the auditory bulla, p10 – zygomatic width, p11 – mastoid width of skull, p12 – postorbital width, p13 – interorbital width, p14 – width of rostrum, p15 – greatest palatal width, p16 – width of the auditory bulla, p17 – width of upper molar M1, p18 – cranial height, p19 – total length of the mandible, p20 – length between the angular process and infradentale, p21 – mandibular tooth-row length, p22 – length of lower carnassial tooth m1, p23 – height of the vertical mandibular ramus, p24 – minimum palatal length, p25 – length of upper molar M1, p26 – length of lower premolar pm2, p27 – length of lower molar m2, p28 – width of lower molar m2, p29 – talonid length of lower carnassial tooth m1, p30 – length of upper canine, p31 – width of upper canine.

Sexual dimorphism in skull size (SSD) is weakly expressed in the genus Meles (Abramov and Puzachenko 2005). Males and females can be reliably distinguished by the size of their canines, which is typical for mustelids and other carnivorous mammals. Because Asian badger subspecies differ in skull sizes (Abramov and Puzachenko 2013), the presumed sex of specimens was assessed separately within “leucurus” and “amurensis” samples. Our sample also included 109 (32%) Asian badger skulls where the sex was unknown. The discriminant analysis was used to find out whether these specimens were males or females. For the assessment of SSD, standard methods of univariate analyses were used. Two indices were calculated for each measurement: SSD1 = 100*[(Mm – Mf)/(Mm + Mf)], SSD2 = 100* [(Mm – Mf)/Mm]. We used SSD1 in (Abramov and Puzachenko 2005) to account for the contribution of both sexes to SSD. Smith (1999) provides some criticism of SSD2, which was widely used. We used it here, both traditionally and because this ratio shows how much larger (in percentage) males are on average than females. The statistical significance of differences between males and females was assessed using the Welch’ F test in Past (Hammer et al. 2001). When interpreting the test results, sample size must be taken into account, as statistical significance is positively correlated with sample size. In other words, low p-values (<< 0.05) do not guarantee the high morphological differentiation that is actually observed. The male and female samples were then examined separately.

First, we formally defined a “morphological system” as the logical intersection of a set of skulls (elements), a set of measurements (variables, which can assume different states) and a set of metrics (methods of measuring similarities or dissimilarities between elements) (Abramov and Puzachenko 2005; Abramov et al. 2009; Puzachenko 2023). The two types of morphospaces were modelled using a non-metric multidimensional scaling method (Davison 1983). The first model describes size variability (SZM, based on a matrix of Euclidean distances) and the second model describes shape variability of the skull (SHM, based on Kendall’s tau-b associations matrix (Kendall 1957)). Some more details see in File S1 in the section “The concepts and applications of morphosystem and morphospace”. In the text, SZM and SHM coordinates of the models are denoted by the letters E (from Euclidean distance) and K (from Kendall’s tau-b metric), respectively. The quality of the description of the raw measurement was assessed by the square of the multiple correlation coefficient of multivariate regression models using the coordinates of the SZM and SHM models as predictors.

(For more details see File S1). The assumption (null hypothesis) of homogeneity of the Asian badger sample was tested using various agglomerative and divisive cluster analysis methods (Rencher 2002). Textbooks usually recommend comparing the classification results obtained by different methods in order to assess the “stability” of the partitioning (Aldenderfer and Blashfield 1984).

Three sets of variables – the coordinates of the SZM or SHM morphospaces individually and in combination – were used to partition the combined M. leucurus and M. anakuma sample. Five different partitioning methods were used. Three of them were non-hierarchical clustering methods: K-means (KM) (Rencher 2002), partitioning around medoids (PAM) (Kaufman and Rousseeuw 1990), and Gaussian finite mixture model with expectation-maximization algorithm (GFM) (Quiles et al. 2005; Scrucca et al. 2016). The next two methods, UPGMA (Sokal and Michener 1958) and McQuitty (MQ) methods (McQuitty 1966), belong to the hierarchical clustering methods (Charrad et al. 2014). We used the following R packages: NbClust 3.0.1 “Determining the Best Number of Clusters in a Data Set” (Charrad et al. 2014), mclust 6.1.1 “Gaussian Mixture Modelling for Model-Based Clustering, Classification, and Density Estimation” (Scrucca et al. 2016), and fpc 2.2–12 “Flexible Procedures for Clustering” (Hennig 2024). All packages implement criteria to find the best/optimal/relevant number of clusters for different methods of sample partitioning. Based on the results of the tests, we were able to estimate the range for the possible optimal number of clusters and the number of clusters that occurred most frequently. The latter was accepted as the presumed most likely optimal partitioning of the sample, least dependent on the clustering method. Note that if the most likely number of clusters was 1, this would be evidence in favor of the hypothesis of sample homogeneity. This hypothesis can also be supported by the random variation in the number of clusters over a wide range. The result of calculations is given in File S1.

We calculated cross-tabulations between formal clusters and the expert division of the sample into subsamples of M. l. leucurus, M. l. amurensis and M. anakuma. Maximum likelihood chi-square (Özdemir and Eyduran 2005; Howell 2011) was used to test the hypothesis of independence of the two compared distributions of specimens across clusters (File S1). In addition, we calculated the correctly classified rate (CCR) in the M. l. leucurus – M. l. amurensis pair and separately in the M. anakuma – M. l. amurensis pair. The CCR is the ratio of correctly classified specimens to the total number of specimens in the groups of badgers being compared. The term “correctly classified” is directly applicable only to the M. anakuma. In all other cases, it means only the agreement between the expert and the formal classifications.

To test the stability of the partitioning of the Asian and Japanese badger sample, we calculated the frequencies of occurrence of each specimen in different clusters under different combinations of coordinates of the SZM or SHM and clustering methods. The set of elements that repeatedly (with a frequency greater than 0.9) fell into the same cluster, regardless of the clustering method or combination of coordinates used, was termed the “kernel”. The kernels were combined in training sample in linear stepwise discriminant analysis with forward and backward variants. The raw skull measurements (p1–p31) were used as predictors. Finally, for every element that was not included in the kernels, the average of posterior probabilities obtained in forward and backward variants of analysis was used as the probability of their assignment to one or another cluster. The purpose of this analysis is to identify skulls with “intermediate” or “outlier” skull measurements.

To assess the position of the Far Eastern badger relative to M. anakuma or M. l. leucurus sensu stricto, we evaluated the positions of all of them in the general context of the morphological variability of all Eurasian badgers Meles. For this purpose, SZM and SHM models were constructed as described above. Both models were then combined into a single model using the principal component method (PCA). Thus the latter model described both size and shape variability in Meles skull.

The PCA of the variance-covariance matrix of log-transformed measurements were used to calculate the multivariate allometric coefficients (MAC) (Jolicoeur 1963, Klingenberg 1996). In this study, the MAC was used to compare males of M. leucurus or M. l. amurensis (File S1).

In this study, we follow a simplified ‘rough’ chronological scheme of Pleistocene Eurasian fossil badger divergence, focusing primarily on the likely ancestors of recent Asian badgers (Koepfli et al. 2008; Jiangzuo et al. 2018). In particular, the data were pooled for the following chronotaxa: 1) “thorali”-like badgers – M. thorali, M. iberica, M. dimitrius, and M. hollitzerii; 2) “chiai”-like badgers – M. chiai, M. minor, M. magnus and M. teilhardi. We collected 723 sites/localities across Eurasia with the remains of Meles dated from the Early Pleistocene to the Early-Middle Holocene. All of them were dated using different qualitative or quantitative methods (File S3). A series of four maps were produced showing badger sites in the Early Pleistocene (~2.5–0.9 Ma; 53 sites), Middle Pleistocene (~0.8–0.13 Ma; 102 sites), Late Pleistocene (~0.13–0.0117 Ma; 337 sites) and Early to Middle Holocene (~11.7–3.0 kyr BP; 233 sites). The choice of such intervals on the Pleistocene timescale was mainly determined by the distribution of palaeontological sites, which increases almost exponentially as we approach modern times. In order to assess the representativeness of the palaeontological data, the maps illustrating badger distribution in the Early and Middle Pleistocene show the points of 1174 localities where no remains of Meles were found (File S3). The generalised boundaries of the ice cover were shown on the maps of the different marine isotope stages (MIS) according to (Batchelor et al. 2019): MIS16, MIS12, MIS6, and MIS2 glaciations.

In addition to the R packages listed above, we used the following software NCSS 12 Statistical Software (2018) (Professional License (Perpetual)) to perform nonmetric multidimensional scaling, Kendall’s tb calculation, cross tabulation procedures, and maximum likelihood chi-squared tests. Free software PAST v. 4.14 (Hammer et al. 2001) used to perform PCA, calculation of Euclidean distance, MAC, descriptive and univariate statistics.

For males and females, the optimal dimensionality of the SZM and SHM morphospace models was two. The low dimensionality indicates highly coordinated variability of individual skull parts in Asian badgers. Despite the low dimensionality of the SZM, its coordinates reproduce well the variation in morphological distances and cranial measurements. The correlation coefficient between morphological distances (Euclid or Kendall’s tb) and the distance between individuals in the SZM morphospace was 0.97 for males and 0.98 for females, and in the SHM model – 0.78 and 0.86, respectively. The mean r² values in multivariate regression models for males and females were approximately the same and equal to ~0.64 (0.12–0.94) (File S1 [table S1]).

In males and females, the E1 coordinate describes the variability due to variation in general cranial dimensions. Indicators of this variability (r > 0.9) in both sexes were condylobasal length, palatal length, maxillary tooth row length, total mandibular length, mandibular tooth row length, and the length between the angular process and the infradentary. The coordinate E2 contains information on variability independent of general size of skull. In females, it was clearly present in the length of Pm2, the length and width of M1 and m2. In males, this variability was generally less pronounced and was observed maximally in width/length of m2, length of talonid of m1, length of m2, width of M1. Thus, in both sexes, dental variability appears partially independent of variations in cranial size. The talonid of the lower carnassial tooth plays the greatest role in this, as its variations do not correlate with the general dimensions of the skull (File S1 [table S1]). Some measurements did not correlate with morphospace coordinates and their variability is mainly interpreted as “stochastic” (postorbital width, interorbital width (in females) and minimum palatal length).

The interpretation of the SHM morphospace is less trivial because its coordinates take into account the variance in skull proportions. Since the correlation coefficients of skull measurements with K1 strongly correlated with similar coefficients for E1 (for males, r = 0.94, for females, r = 0.86), it can be argued that K1 mainly describes allometric variability associated with variation in skull size. The coordinate K2 in females and males describes the nonlinear variability of the relative sizes (a ratio of measurement value to condylobasal length) of the postorbital width, interorbital width (females only), minimum palatal length (females only) and mastoid width. Non-linearity in this case means that the variability is not described by a linear regression but can be reproduced by a polynomial regression. Further discussion of this topic is outside the scope of the present study. From the above results, we concluded that the SZM and SHM models contain most of the information about the variation in initial measurements in the sample of skulls of the Asian badger studied. Therefore, the morphospace coordinates can be used to test the hypothesis of homogeneity of the Asian badger sample.

For the joint sample of M. l. leucurus, M. l. amurensis and M. anakuma, a negative correlation was found between E1 and K1 with geographical longitude (r² are 0.50 and 0.49; p < 0.01). Formally, E1 is weakly but statistically significantly correlated with latitude (r² = 0.3, p < 0.01). However, there is no reason to take these correlations as evidence of cline (gradual) geographic variability. Between 45ºE–130ºE and 38ºN–57ºN, E1 and K1 do not correlate with either longitude or latitude. The above correlation is only due to the abrupt changes in the pattern around 130ºE and 45ºN. These changes are interpreted as follows: to the west of 130ºE and to the north of 45ºN, badgers with larger skulls are distributed, whereas to the east and south of this longitude/latitude, badgers with smaller skulls are found. These geographical forms also differ in skull proportions.

The SZM and SHM coordinates (E1, E2, K1, K2) were used separately and in combination to partition samples of male and female Asian badgers. Detailed results described in the File S1. The optimal/relevant number of clusters obtained by combining all coordinates (E1, E2, K1 and K2) was two, independently of the cluster analysis methods. When the SZM or SHM coordinates were used separately, the optimal number of clusters was two (63%), three (22%) or four (15%), depending on the partitioning method and the sex of the badgers. This instability in classification results generally reflects the fact that within Asian badgers, including the Japanese badger, there is no hiatus in skull measurements between geographic forms.

Comparison of the a priori classification of M. anakuma, M. l. leucurus and M. l. amurensis with the formal partitioning of the sample by the coordinates of the SZM and SHM models (File S1 [crosstabulation results sections]) showed that the similarity between the two partitions was not random in all cases (p < 0.0001; Table 2). This result allows us to reject the null hypothesis of homogeneity of the Asian badger sample.

The mean CCR (Table 2) for the pair of M. l. amurensis– M. l. leucurus was 0.84 (males) and 0.86 (females). In the M. anakuma – M. l. amurensis pair, this value was low – 0.25 (males) and 0.26 (females). When the SZM or SHM model coordinates were used separately, the general tendency to group Far Eastern and Japanese badgers into one cluster was generally maintained, as was the tendency to discriminate between M. l. amurensis and M. l. leucurus (File S1). The obtained sample partitioning demonstrates the high morphometric specificity of Far Eastern badgers with respect to other M. leucurus and, in parallel, brings them closer to the island species M. anakuma.

The division of Asian badgers into two clusters is the most probable (in terms of frequency of occurrence) and optimal according to most criteria. In this case, Japanese badgers are included in one cluster with the Far East badgers and contrasted with the rest of the continental Asian badgers. The proportion of specimens included in the “kernels” of two clusters was about 67% in the male sample and 88% in the female sample (File S2 [tables S4, S5]). Stepwise discriminant analysis was used to classify specimens not included in the cluster’s “kernels” (Fig. 2). In the male sample, 89.7% of the specimens had a posterior probability of membership in one of the two clusters greater than 0.90. This proportion was higher for females – 0.98.5%. As a result, we refined the a priori classification of continental Asian badgers into M. l. leucurus and M. l. amurensis. About 10% of the skulls had morphometric characteristics intermediate between M. l. leucurus and M. l. amurensis. These specimens are of particular interest because their presence in the sample leads to the observed instability of the classification.

Figure 2. 

Scatter diagram of canonical axes (Axis 1) obtained by stepwise discriminant analysis using forward and backward methods. Cluster 1 corresponds to Meles leucurus amurensis and M. anakuma, cluster 2 – to M. l. leucurus. A males, B females.

Identification of specimens characterised by intermediate cranial features between M. l. amurensis and M. l. leucurus

Males and females for which the posterior probability was below 0.9 were divided into two groups. They included specimens from a region where an overlap between the ranges of M. l. leucurus and M. l. amurensis cannot be excluded at the present time, i.e. is somewhere ~113–134ºE (Fig. 3). The group “amurensis–leucurus” included 6 skulls (all males) that were similar to the representatives of M. l. amurensis, but in some features they were close to the M. l. leucurus. The group “leucurus–amurensis” included 9 specimens (8 males and 1 female) that were more similar to the skulls of M. l. leucurus, but also in some features resembled the skulls of M. l. amurensis. We labelled the next group as “deviant”, which included skulls of 12 badgers (10 males, 2 females) from the western part of the range where the presence of M. l. amurensis is highly unlikely. They were either consistently classified as “amurensis” with a high probability, or had a low probability of being assigned to the morphological cluster associated with “leucurus”.

Figure 3. 

Distribution of Asian badgers classified as “amurensis” or “leucurus” and locations of the “intermediate” and “deviant” individuals. The approximate range of the Asian badgers from Abramov (2016).

Morphometric description of M. l. leucurus and M. l. amurensis

In terms of general skull size, all Asian badgers form a series (from smaller to larger): Meles anakuma < M. l. amurensis < M. l. leucurus (File S1 [tables S2, S3]). Differences between the males of M. l. leucurus and M. l. amurensis have a high level of statistical significance and cover the entire skull and mandible with the exception of cheek teeth but not canines: width/length of M1, width/length of m2, talonid length of m1, length of lower premolar Pm2 and upper Pm4. Thus, M. l. amurensis has relatively larger premolars and molars compared to M. l. leucurus. The variability described above can be interpreted as geographical, but not as cline. It should be emphasised that no morphological hiatus between M. l. amurensis and M. l. leucurus was found for any of the measurements (example in Fig. 4).

Figure 4. 

Bivariate plot of condylobasal length and length between the angular process and infradentale of mandible: A males, B females. Ellipses show 95% prediction interval for regression of variables.

“Intermediate” and “deviate” specimens

There are several specimens in our collection that occupy an intermediate position between M. l. amurensis and M. l. leucurus (Fig. 5; Supplementary Data SD1, Table S2 and Table S3). The proportion of “intermediate” specimens was 6.66% in the male sample of M. l. amurensis and 7.54% in the male sample of M. l. leucurus. Only one “intermediate” skull was found among M. l. amurensis females. Such specimens were absent in a relatively small sample of the females of M. l. leucurus.

Figure 5. 

Cauterised box and whiskers plots of selected skull measurements (males) in the M. l. amurensis (M. l.a.), M. l. leucurus (M. l.l.), “intermediate” groups “amurensis-leucurus” (“a-l”), “leucurus–amurensis” (“l-a”) and “deviate” group of M. l. leucurus (“d”).

On average, these skulls were slightly larger than those of the typical “amurensis”, but slightly smaller than those of the typical “leucurus”. The average condylobasal skull length of males from the “amurensis–leucurus” group is large significantly (p = 0.0003, according to Welch F’ test) than that of M. l. amurensis and its minimum and maximum values are close to the upper limit of this measurement in M. l. amurensis. A similar variation was described for cranial mastoid width, interorbital width, length between angular process and infradentale, and several other measurements. The maximally pronounced (p < 0.000001) difference from the typical “amurensis” was found for the mandibular length. In addition, it is worth noting that the interorbital width in this “intermediate” group was the same as in the typical “leucurus” (Fig. 5).

Specimens from the “leucurus–amurensis” group were smaller than typical “leucurus” in a number of measurements, most notably, in the condylobasal length (p < 0.000001) and the mastoid width of the skull (p < 0.00005). The measurements in this group were slightly larger or almost equal to those in the “amurensis–leucurus” group, but in some cases they were slightly smaller (e.g. the width of auditory bulla; Fig. 5).

The interesting group of “deviate” skulls of M. l. leucurus differed from the typical “leucurus22.5.2025 г. 15:10:24” in smaller cranial dimensions, mainly because of a shortening of the viscerocranial part of the skull. The most significant differences (p < 0.000001) were observed for the length of the mandible. It is also worth emphasising that the small distance between orbits was very close to that of M. l. amurensis (Fig. 5). The characteristics of tooth rows were relatively stable, i.e. their variability was little or practically not affected by the variability of overall cranial dimensions. Such size independence is most clearly expressed in the length variability of the second lower molar (Fig. 5).

Comparison of multivariate allometric patterns in typical M. l. amurensis and M. l. leucurus

MAC values were between 0.5 and 1, indicating weak negative allometry or isometry), but this does not apply to dental measurements. In M. l. amurensis, the length of the upper carnassial tooth is characterised by isometry or a weak positive allometry. In M. l. leucurus, this measurement shows moderately high positive allometry, i.e. the tooth length increases slightly faster than the general cranial dimensions. In M. l. leucurus, the same positive allometry is found for the width of the upper molar M1 and the length of the lower carnassial tooth m1. At the same time, the length of the upper molar M1 does not increase at the same rate, so that this tooth is relatively shorter and wider in larger skulls. The talonid of the m1 is characterised by a strong positive allometry, which is higher than that of its length. Because of this allometry, larger animals have a relatively larger talonid. The length of the lower premolar pm2 shows a very strong positive allometry, especially in M. l. amurensis.

The high individual variability of the morphometric parameters of pm2, M1, m1 and m2 were reflected in the very wide 95% confidence intervals of the corresponding MACs. The correlation coefficient between the MACs of M. l. amurensis and M. l. leucurus was relatively low (0.73, p < 0.00001, r² = 0.53).

We calculated MACs for different geographical samples of the European badger M. meles (Scandinavia: Norway and Sweden; Continental Europe) and the Southwest Asian badger M. canescens (Abramov et al. 2009; Abramov and Puzachenko 2013). Then we developed a PC model based on a variance-covariance matrix of the MACs including MACs of M. l. amurensis and M. l. leucurus (Fig. 6). The variability of the allometric pattern in badgers depends mainly on the allometry of the length of pm2 (p26), the length of m2 (p27), the talonid of the lower carnassial tooth (p29) and the length of Pm4. Within this variability, Asian badgers (M. l. leucurus, M. l. amurensis) differ from the other Meles.

Figure 6. 

Biplot of the first two PCs calculated from the variance-covariance matrix of the MACs of several Meles taxa. The descriptions of the skull measurements (p1–p31) are shown in Figure 1.

Position of M. l. amurensis relative to M. anakuma in the context of morphological variability in Eurasian badgers

To assess the degree of divergence in the M. l. amurensis – M. anakuma pair and the M. l. amurensis – M. l. leucurus pair, SZM (three coordinates) and SHM (two coordinates) morphospaces were constructed for all Meles. The PC model combining SZM and SHM (Fig. 7) is interpreted as a metamorphospace of cranial size and shape variations. PC1 correlates with overall skull size (coordinate E1) and explains 78.4% of the variance. All Meles taxa are arranged along PC1, from the smallest Japanese badger to the largest European badger. PC2 correlates mainly with the K2 coordinate. The K2 correlates (r > 0.5) with the width of the mastoid process, the width of auditory bulla, and somewhat less (0.5–0.4) with the length of the auditory bulla and the height of skull. PC3 correlates with the coordinates E3 and K1. The E3 contains information on the variability of the talonid length of m1 and the length of m1 (r > 0.5). The K1 describes the allometry due to the variation of cranial size and is strongly correlated (r > 0.7) with many cranial measurements – viscerocranial length, palatal length, skull height, total length of the mandible, maxillary row length, length of the lower premolar pm2, etc.

Figure 7. 

Biplots of PC1×PC2 (A) and PC1×PC3 (B) illustrating size/shape variability of the skull of Eurasian badgers, Meles ssp. E1-E3, K1, K2 - coordinates of SZM and SHM morphospaces.

The separation of M. anakuma from M. l. amurensis was observed along only the PC3. Consequently, the degree of morphological difference between M. l. amurensis and M. anakuma should be recognised as small one. Therefore, the Japanese badger is more closely related to M. l. amurensis than to M. l. leucurus.

Pleistocene and Holocene paleogeography of Eurasian badgers

(For details see File S3 [tables S6, S7].) The western (“thorali”) and eastern (“chiai”) lineages of the badger were geographically isolated at least since the Early Pleistocene (Fig. 8) (~2.59 Ma). The absence of badger remains at this time in Eastern Europe, Middle and Central Asia, and West Siberia does not seem to be a consequence of insufficient study. The Prorva site (layer 6, Kazakhstan) is the only known find of Meles sp. in southern West Siberia (Gaiduchenko 1986, cited by Shpansky 2018). The fauna of this stratum is correlated with the Podpusk-Lebyazhe faunal complex, dated to the Matuyama chronozone, below the Reunion subchronozone, i.e. > 2.2 Ma (Vislobokova 1996). During the Middle Pleistocene, the two badger lineages continued to evolve independently (Fig. 8). Badgers of the eastern lineage lived within the presumed modern range of M. m. amurensis. There are known finds of Meles from the Kurtak (= Berezhekovo) site in the upper Yenisei River (Russia), which are restricted to strata roughly corresponding to the MIS9 (0.38– 0.30 Ma) or MIS7 (0.24–0.19 Ma) interglacials (Chlachula 2001; Haesaerts et al. 2005; Malikov 2018). The badger colonisation of Japanese islands most probably occurred during the MIS12 glaciation (0.48–0.42 Ma). This is supported by a find in the Matsugae Cave on Kyushu Island (Ogino et al. 2009; Handa et al. 2021). The Matsugae fauna was correlated with the Quaternary Mammal Zone 4 of the Japanese Islands (0.5–0.3 Ma) and the Asian Choukoutien fauna (Loc. 1) (Ogino et al. 2009).

Figure 8. 

The Early and Middle Pleistocene sites with remains of Meles in Eurasia (Supplementary Data SD3).

During the Late Pleistocene (0.13 – 0.01 Ma), the range of the Asian badger expanded. This is supported by numerous finds in southern Siberia (Fig. 8). In particular, badgers appear to have lived permanently in the area of the Sayano-Altai Mountain Region and in Transbaikalia. These areas can be considered the zonal refugia in the Last Glacial Maximum (LGM, MIS2). We assume that at least during the LGM badgers occurring in the ranges of M. l. amurensis and M. l. leucurus were isolated from each other. At that time, the geographical boundary between them could be the Great Khingan Ridge.

Radical changes in the distribution of M. l. leucurus occurred during the Holocene (<0.01 Ma) when badgers dispersed to the south of Western Siberia (Devjashin et al. 2017) and Central Asia (Fig. 8). The species reached the Ural Mountains and the Volga River basins in the late Holocene (Kinoshita et al. 2020). In the latter case, they displaced the M. meles (Gasilin and Kosintsev 2010), which had lived there since at least the Last Interglacial (MIS5, ~0.125 Ma) (Fadeeva et al. 2010; Danukalova et al. 2020; Gimranov and Kosintsev 2020). Note, that this expansion is still in progress today (Oparin et al. 2024). It can be assumed that the range expansion of M. l. leucurus occurred simultaneously in the eastward direction, in the Amur River basin. As they moved in this direction, they could have reached the range of the Far Eastern M. l. amurensis.

As it is evident from the brief review above, in Eurasia the Meles range became continuous only in the Holocene. It is highly likely that the separation of the ancestors of M. anakuma from those of M. l. amurensis occurred no later than in the Middle Pleistocene (~MIS12). Thus, the onset of divergence between M. l. leucurus and M. l. amurensis was probably not earlier than in the MIS8 (0.30–0.24 Ma).

In the present study, we accepted the species taxonomic status of the Japanese (M. anakuma), European (M. meles) and Southwest Asian (M. canescens) badgers as the fact established in previous studies (Del Cerro et al. 2010; Abramov and Puzachenko 2013; Proulx et al. 2016). We critically examined the cranial variations of the Asian badger to test the hypothesis of a greater morphological separation of the badgers living in the southern Far East, at least within the borders of Russia, from the other continental Asiatic badgers. It is generally accepted that these badgers represent a subspecies (M. l. amurensis) that are characterised by a set of traits that distinguish it from another subspecies, M. l. leucurus, which is distributed in the western part of the species range. Several previous studies have noted the morphological distinctiveness of Far Eastern badgers (Petrov 1953; Heptner et al. 1967; Baryshnikov and Potapova 1990; Abramov and Puzachenko 2013). Analysis of body size, colour, cranial features and sexual dimorphism of badgers from the Russian Far East and North Korea showed significant differences from M. l. leucurus. The badgers of the Russian Primorye have a peculiar colouring (Fig. 9). The whole of their coat is dark, with the brown colour dominating, but a high variability was observed. The snout is also dark brownish, and in some specimens the facial stripes are almost invisible. When they are, the pattern of the facial mask resembles that of typical Asian badgers. North Korean badgers are similarly coloured (Abramov 2003). The dark colour and small size are considered the main characteristics of Far Eastern badgers. The absence of upper and lower first premolars is another characteristic of M. l. amurensis. The first premolars are also sometimes missing in M. l. leucurus, M. canescens, and M. meles (very rarely), but these teeth have never been observed in Far Eastern badgers (Baryshnikov and Potapova 1990; Baryshnikov et al. 2003).

Figure 9. 

The head colour pattern of the Eurasian badgers. A Variation head colour among different species and subspecies adopted from (Abramov 2003); B Far Eastern badger (M. l. amurensis) (published with the permission of the Public Relations Department of the “Land of the Leopard National Park”, Primorsky Krai province, Russia).

In this study, we tested the “neutral” null hypothesis of homogeneity of the Asian badger (M. leucurus) sample. This hypothesis was tested using different clustering techniques to demonstrate the stability/instability of partitioning the sample into clusters that would be essentially invariant with respect to the techniques used. As a result, two subsets were found, in the male and female samples, whose elements were consistently distinguished by all clustering methods. Based on this fact, the hypothesis of sample homogeneity was rejected. One subset included the Far Eastern M. l. amurensis and the Japanese M. anakuma, and the second subset consisted only of the Asian M. l. leucurus.

Here we confirmed statistically significant differences in size of skull between Far Eastern and other continental Asian badgers. The relatively small size of M. l. amurensis brings it closer to the Japanese badger, the smallest of the Meles species. A certain intrigue in the evolutionary relationship between the continental badgers and Japanese badger arises from the assumption that the latter are descended from the continental form, whose descendants are currently represented by the subspecies M. l. amurensis. Unexpectedly, only a small proportion of specimens had skull measurements intermediate between M. l. amurensis and M. l. leucurus. As they came from the area where the ranges of two subspecies of Asian badgers are likely to overlap, it is hypothesised that the “intermediate” specimens may be hybrids between two subspecies. The difference in allometric patterns between M. l. amurensis and M. l. leucurus did not exceed the level of interspecies differences and was smaller than between different samples within the European badger, M. meles. Sexual size dimorphism is weakly expressed in the skulls of Eurasian badgers (Abramov and Puzachenko 2005). In he badgers of continental Asia, the dimorphism varies within the limits characteristic to Meles. The relatively high SSD indices obtained for the Japanese badger need to be verified with a larger sample size than that used in this study.

Finally, we combined morphometric data for all Meles species and assessed the extent of differences between M. l. amurensis and M. l. leucurus, and between M. l. amurensis and M. anakuma. It was shown that M. l. amurensis is indeed morphologically closer to the M. anakuma than to the M. l. leucurus. However, it should be noted that the variability of craniometric traits is primarily continuous within the genus Meles. The morphometrical “niches” of species essentially overlap relative to the main coordinates of model morphospaces (Fig. 7). The absence of clear dichotomies between species had long served as a justification for the monotypic of the genus Meles (Wozencraft 2005).

The divergence of the eastern and western evolutionary lineages has been observed since at least the late Pliocene, but genetic (mtDNA) data do not support substantial divergence between M. leucurus sensu lato and M. meles. Kinoshita et al. (2019) found genetic evidence of hybridisation between both species in the area of modern range overlap. In the sympatry area of M. leucurus with M. canescens in Central Asia (Tien Shan, Turan Lowland), several specimens have been found that can be attributed to “canescens–leucurus” hybrids based on morphological features of the skull (Abramov and Puzachenko 2007). Interestingly, the skulls of the presumed hybrids were more similar to M. leucurus than to M. canescens. Therefore, we do not exclude the possibility that the so-called “deviant” specimens found during this study could be of hybrid origin. However, genetic analyses are needed to test this hypothesis.

The results from the studies of genetic diversity of mtDNA and Y chromosomal genes within M. leucurus did not support genetic differentiation between M. l. leucurus and M. l. amurensis based on the set of markers investigated so far (Tashima et al. 2011a, 2011b; Koh et al. 2014; Lee et al. 2016; Luo et al. 2016; Kinoshita et al. 2017). However, the same studies confirmed the genetic differentiation of the island species M. anakuma. Thus, the results of genetic studies to date contradict the evidence for morphometric geographic variability in continental Asian badgers.

The radiation of the genus Meles occurred during the general climatic changes that took place during the latest Pliocene and the beginning of Pleistocene as a result of environmental shifts across Eurasia (Faggi et al. 2024). We can speculate that the interval of ~3.5–2.6 Ma could be the time of divergence between the western (“thorali”) and eastern (“chiai”) lineages of Meles because Arctonyx (Hog badgers) and Meles diverged in Asia between 4.4 and 3.6 Ma (Koepfli et al. 2008). The find of early “thorali”–like badger (M. cf. thorali) in Europe was dated to the Villafranchian, MNQ16 (Mammal Neogene Quaternary scheme) (Almenara-Casablanca-4, Spain) (Madurell-Malapeira et al. 2009) or slightly later (Luo et al. 2016). The presence of only M. chiai in the Late Pliocene of China (Jiangzuo et al. 2018) is in agreement with this hypothesis. A tentative chronological “boundary” between M. thorali and M. meles was accepted as about 1.5 Ma (Faggi et al. 2024). It is also possible that M. canescens could have evolved in West Asia at about the same time (Luo et al. 2016).

Jiangzuo et al. (2018) suggested that M. magnus could have been the ancestor of M. leucurus and the latter species did not originate until the mid-Early Pleistocene on the territory of modern China. Based on the palaeontological evidence, we can tentatively assume that the Middle Pleistocene ancestors of M. l. amurensis could have migrated to the Japanese islands via the Korean Peninsula during the period of very significant global sea-level decrease, approximately during the MIS12 glaciation (Yang et al. 2023). However, it could be an earlier migration at the end of the Early Pleistocene (~1 Ma, ~MIS28), if the results by Luo et al. (2016) are taken into consideration.

Our results provide indirect evidence for the proximity of modern M. l. amurensis to the hypothetical Pleistocene ancestor of the Japanese badger, which occurred in eastern and southeastern China in the Middle Pleistocene. In particular, they are closer to the Japanese badger in terms of skull size. Besides, the island and Far Eastern badgers are also closer to each other in general skull shape than to the M. l. leucurus.

We proposed that an ancestor of the recent M. l. leucurus arrived to southern Siberia from the Far East in the Middle Pleistocene during the MIS9 interglacial (0.34–0.30 Ma). The origin of M. l. leucurus took place under the conditions of periodic geographical isolation and in harsher mountainous conditions, with more intense climatic fluctuations during the end of Middle Pleistocene (Lisiecki and Raymo 2005) than those observed in the ancestral range of M. l. amurensis. It is possible that the isolation of M. l. leucurus from M. l. amurensis occurred during glacials of the MIS8 (0.30–0.24 Ma) and MIS6 (0.19–0.13 Ma) and in the first phase of the Last Glaciation (MIS4; 0.071–0.057 Ma) or during the Last Glacial Maximum (MIS2; 0.029–0.017 Ma). The isolation of M. l. leucurus from M. l. amurensis could have been interrupted during interglacial periods, for example, during MIS7 (0.24–0.19 Ma BP) and MIS5 (0.130–0.071 Ma), or even during the megainterstadial MIS3 (0.057–0.029 Ma). Due to several breaks in geographical isolation, putative genetic exchanges between M. l. leucurus and M. l. amurensis could reconcile the results of genetic and morphological studies. However, we have no palaeontological evidence for such events.

It is possible to imagine how the geographical isolation between M. l. leucurus and M. l. amurensis could have been broken in the past. The supposed “contact zone” between them in the Amur Basin is quite wide (Fig. 3). A relatively small number of skulls with intermediate values of morphometric characters were found in the investigated sample, which could belong to the hybrids between M. l. leucurus and M. l. amurensis. It is also likely that M. l. leucurus has already penetrated in to the “original” range of M. l. amurensis in the Far East.