We focus only on specimens harbouring the most widespread mtDNA lineage 8 (or ‘green’) of N. n. scutata, as defined by Fritz and Schmidtler (2020) and Asztalos et al. (2021). Grass snakes were sampled in different areas across its range from 2005 to 2023 (Fig. 1A; Table S1). Tissue samples consisted of buccal swabs, pieces of ventral scales and tail tips or intercostal muscles (for road-killed animals). All snakes were released at the place of capture immediately after tissue sampling. The samples were stored in 95% ethanol at 4°C. Total genomic DNA was isolated by standard proteinase K and phenol–chloroform protocols (Sambrook et al. 1989) or a salt extraction method (Aljanabi and Martinez 1997).

We determined partial nucleotide sequences of the mitochondrial cyt b gene to assess haplotype diversity and infer phylogeographic patterns. A target fragment (752 bp long) was amplified with the primers L15162Vb and H15914Vb (Ursenbacher et al. 2006). PCRs were carried out in 25 µl reaction mixtures each containing 15–60 ng of DNA, 0.2 mM each dNTP, 2.5 µl of 10× amplification buffer (10 mM Tris-HCl pH 8.5, 50 mM KCl and 2.5 mM MgCl₂), 1 U of Taq DNA polymerase and 5 pmol of each primer. Amplification was performed on a T100 Thermal Cycler (Bio-Rad) with an initial denaturation step of 3 min at 94°C followed by 35 cycles of 60 s at 94°C, annealing at 54°C for 30 s and 60 s extension at 72°C; with a final extension of 5 min at 72°C. We cleaned PCR products either by treatment with a combination of exonuclease I and shrimp alkaline phosphatase or using spin columns (GeneJET PCR Purification Kit, Thermo Scientific). Sequencing was carried out on an ABI 3130xl automated capillary sequencer (Applied Biosystems) with ABI Prism BigDye Terminator Cycle Sequencing Ready Reaction Kit 3.1 by means of the same forward and reverse primers as for PCR.

The cyt b sequences available in GenBank (NCBI) and belonging to the N. n. scutata clade along with published information on haplotypes (Kindler et al. 2017; Asztalos et al. 2021) were included in the analysis (Table S1). The sequences were aligned manually and checked for unexpected stop codons using SEAVIEW 4.5.4 (Gouy et al. 2010). A haplotype network of N. n. scutata was generated via the TCS algorithm implemented in POPART (Leigh and Bryant 2015). Subclades of interest were designated within the resulting network according to criteria of a star-like pattern (a high-frequency ancestral haplotype plus descendants) and to the geographical distribution (presence in areas that were supposedly colonised after the LGM). Average uncorrected genetic distance between sequences and its standard error based on 10,000 replications were calculated in MEGA11 (Tamura et al. 2021).

Geographic structure of cyt b variation in N. n. scutata was examined in BARRIER 2.2 (Manni et al. 2004). This software uses Monmonier’s maximum difference algorithm (Monmonier 1973) to identify genetic barriers, i.e., the zones where differences between pairs of populations are the largest based on the provided distance matrix and geographical coordinates. Although this approach requires the definition of small-sample groups a priori, it allows for the identification of geographical structuring within a group of DNA samples a posteriori. We calculated the first 10 barriers (the significance of barriers was expected to decrease with their rank) within a set of 27 ‘populations’ predefined manually based on their geographical grouping and a minimum allowed sample size of 2. The matrix of distances between the populations was computed by means of toolbox SPADS 1.0 (Dellicour and Mardulyn 2014) as IID1 (interindividual distance 1): an interindividual distance based on allelic frequencies (Miller 2005). To visualise the spatial pattern of cyt b diversity, interpolation surfaces of nucleotide diversity (Pi) were generated with R functions slidingWindowPi and GDivPAL from the SPADS toolbox. As the first step, slidingWindowPi assigned a value to a given raster cell where the Pi estimated for the group of sequences sampled within a circle centred on this cell (was set to 25 km). Then, output was fed into GDivPAL to build interpolating surfaces with an inverse distance interpolation parameter (a) set to 5. The resulting raster was visualised and clipped by means of the approximate limits of the N. n. scutata range in QGIS 3.30 (QGIS Development Team 2023).

For the total dataset, for the haplogroups discovered in the haplotype network and for the geographical groups of samples suggested by BARRIER, we calculated a number of diversity and population subdivision parameters. The number of haplotypes, haplotype diversity and nucleotide diversity (Hd and Pi, respectively) were computed in DNASP 6.10.04 software (Librado and Rozas 2009). The number of private haplotypes (h_P) was assessed with ARLEQUIN 3.5.2.2 (Excoffier and Lischer 2010). Demographic expansion was evaluated by Tajima’s (1989) D test and Fu’s (1996) F_S test of neutrality as implemented in DNASP and via examination of the frequency distribution of pairwise differences between sequences (mismatch distribution) (Harpending 1994) in ARLEQUIN. The observed distribution was compared with the expected distribution under a model of sudden demographic expansion using the sum of squared deviations (SSD) test and Harpending’s raggedness index tested with 10,000 permutations.

To evaluate the hypothesis of suitable environmental conditions (i.e., potential refugia) for N. n. scutata in Siberia and adjacent areas during the LGM, SDMs for current and past climatic conditions were generated using MAXENT 3.4.4 (Phillips et al. 2006). We applied recent methodological recommendations to compute a robust SDM, e.g., occurrence filtering (with MAXENT), the use of multiple combinations of model parameters (features and regularisation multipliers), and multiple statistical criteria for model selection. A total of 4493 localities with known presence of N. n. scutata were utilised, comprising our own records, museum collections and previously published (Litvinchuk et al. 2013; Zima and Fedorenko 2019) and highly precise open-source data, i.e., GBIF https://doi.org/10.15468/dl.gnm233 and https://reptilia.club (Table S2).

To compute the model under the current climatic conditions, 19 bioclimatic layers representative of the climatic data for the ~1950–2000 period were extracted from the WorldClim database (http://www.worldclim.org). Layers had 2.5-arc-minute spatial resolution, and all analyses were carried out under the WGS84 projection with a species-specific mask covering the area of occurrence of this subspecies and adjacent regions (from 36°N to 65°N and 12°E to 120°E). To eliminate predictor collinearity before generating the model, we calculated Pearsons’s correlation coefficients for all pairs of 19 bioclimatic variables by means of ENMTOOLS (Warren et al. 2010). For correlated pairs (|r| > 0.75), we excluded the variable that appeared the least biologically important for grass snakes. The resulting dataset contained eight bioclimatic variables: Bio 1 (annual mean temperature, °C, ×10), Bio 3 (isothermality; Bio 2/Bio 7 × 100), Bio 4 (temperature seasonality; CV × 100), Bio 8 (mean temperature of wettest quarter, °C, ×10), Bio 13 (precipitation of wettest month, mm, ×10), Bio 14 (precipitation of driest month, mm), Bio 15 (precipitation seasonality; CV) and Bio 19 (precipitation of coldest quarter, mm). A jackknife analysis was performed to estimate relative contributions of variables to the MAXENT model.

We ran MAXENT for 10 replicates with 30% random test percentage testing. Model calibration consisted of the evaluation of models created with distinct regularisation multipliers (0.5 to 6.0 at intervals of 0.5) and feature classes (that resulted from all combinations of linear, quadratic, product, threshold and hinge response types). The number of background points was 20,000. The best parameter settings were selected considering statistical significance (partial receiver-operating characteristic; ROC), predictive power (omission rates E = 5%) and a complexity level (Akaike Information Criterion, corrected; AICc) obtained using R package kuenm (Cobos et al. 2019). Additionally, model performance was evaluated with the help of the area under the curve (AUC; ranging between 0 and 1) and the true skill statistic (TSS; ranging from –1 to +1) of the 10^th percentile training omission threshold (Allouche et al. 2006). The ClogLog output format (ranging between 0 and 1) was chosen for processing resulting maps (Phillips et al. 2017). To project a current ecological niche onto climatic conditions during the Late Quaternary, we applied bioclimatic layers with 2.5-arc-minute spatial resolution from the PaleoClim database (http://www.paleoclim.org), which represent the climatic data during 4.2–0.3 ka (late Holocene), 8.3–4.2 ka (mid-Holocene), 11.7–8.3 ka (early Holocene), 12.9–11.7 ka (late Pleistocene, Younger Dryas Stadial), 14.7–12.9 ka (late Pleistocene, Greenland Interstade 1), 17.0–14.7 ka (late Pleistocene, Heinrich Stadial 1), ~21 ka (LGM), the last interglacial (~130 ka) and the mid-Pleistocene (~787 ka; Bio 3 is lacking in the dataset).

We successfully obtained cyt b sequences for newly collected tissue samples from 69 N. n. scutata from 49 geographical locations. The total dataset consisted of 135 individuals from 91 sampling sites across the range of this subspecies (Fig. 1A; Table S1). There are 95 genotyped localities for grass snakes belonging to lineage 8 of N. n. scutata, including ND4 data from GenBank (Fig. 1A). The final length of the trimmed alignment was 658 bp. A total of 38 haplotypes were found (Hd = 0.705 ± 0.044), while the number of polymorphic sites was 40, the average number of nucleotide differences 1.54 and nucleotide diversity 0.0024 ± 0.0003. The uncorrected average p distance between N. n. scutata haplotypes in the dataset was 0.26% ± 0.08% (range: 0.00–1.52%). Because of shorter cyt b sequences compared to previous studies (Guicking et al. 2006; Kindler et al. 2013, 2014, 2017; Asztalos et al. 2021) on N. natrix, it was not possible to directly assign haplotype names to the previously discovered ones. Therefore, new names were introduced (see Table S1). Nonetheless, during a comparison of the full-length sequences to the shortened ones, only a small decrease in the number of haplotypes was observed (4.8%; see ­Table S1).

Most haplotypes formed a star-like cluster in the haplotype network (Fig. 1B). The most frequent haplotype n1 (corresponding to gn1 in Asztalos et al. 2021) proved to be shared by 53.7% of all DNA samples and was found on a region from Gotland Island (Sweden), Kaliningrad Oblast (Russia) and eastern Poland in the west to the east bank of the Yenisei River in Siberia (Krasnoyarsk Krai, Russia) in the east, corresponding to ~4400 km in the direction of longitude. This is the only haplotype found in West Siberia among 25 specimens from 13 sampling sites (Fig. 1C). In contrast, the only two analysed specimens from Transbaikalia (East Siberia) had unique haplotypes that differed from the most common haplotype by one and two substitutions. The haplotypes mostly found on the East European Plain differ from n1 by one to three nucleotide substitutions and are hereafter referred to as haplogroup S1. In contrast, the Black Sea Region proved to not be dominated by any single haplotype and actually harbours many unique and distant haplotypes. Moreover, the most common haplotype n1 was not detected south of the Main Caucasian Range (MCR).

Two other minor star-like groups of haplotypes could be identified in the haplotype network (Fig. 1B). Haplogroup S2 was predominantly found in west Transcaucasia and along the south-eastern coast of the Black Sea but also in southern Finland, Central Russia and southern Ukraine. Haplogroup S3 was identified from the North Caucasus, South Russia, West Kazakhstan and Lithuania. In both cases, the northernmost haplotypes are the most frequent and central haplotypes (=ancestral) in each cluster (Fig. 1C). Additionally, the structure of the haplotype network suggested that both most frequent haplotypes in S2 and S3 are candidates for being ancestral haplotypes for S1. This result is concordant with their more southern distribution.

Monmonier’s maximum difference algorithm implemented in BARRIER indicated that the MCR is the most important barrier structuring diversity within N. n. scutata, thus dividing it into two groups (hereafter referred to as northern and southern). The results of the analysis involving 10 inferred barriers are given in Figure 2. It shows that most of the genetic structure within lineage 8 is associated with southern parts of its range. Nonetheless, three barriers were detected outside this area: between the Caspian Depression and the Mugodzhar Range, between the Mugodzhar Range and West Siberia and between West and East Siberian populations (Fig. 2).

Several measures of genetic diversity were calculated for the population groups inferred by BARRIER (­Table 1). The southern group had significantly higher cyt b diversity as compared to the northern group. Notably, these two groups have only one common haplotype. Nucleotide diversity was interpolated across the range of this subspecies, thereby revealing the highest diversity levels in West Transcaucasia, southern Ukraine and Crimea (Fig. 2). Areas with elevated diversity were also noticed in western Kazakhstan in the Mugodzhar Range and in Transbaikalia. It should be noted that these areas are poorly sampled, and the level of nucleotide diversity is likely an overestimate. Nonetheless, diversity was clearly higher in these areas than in neighbouring sampled areas (Table 1).

The results of Tajima’s D test and Fu’s F_S test were negative and significant for the northern group and all specimens combined (Table 1). For the southern group, both statistics were negative, but Tajima’s D was not significant. Mismatch distributions were unimodal for both groups (Fig. 3), with no significant SSD and raggedness index values. These observations possibly indicate a recent sudden demographic expansion of this subspecies.

To calibrate SDMs, we assessed 372 replicates, all of which were statistically significant when compared with a null model of random prediction. Of these significant models, all met the omission criterion of 5% and only two models were statistically significant among models meeting the AICc. Performance metrics for parameter settings used for creating final models are given in ­Table S3. Among the parameters included in the model, Bio 14 (precipitation of driest month), Bio 1 (annual mean temperature), Bio 3 (isothermality) and Bio 15 (precipitation seasonality) made the largest percentage contributions (31%, 21%, 14% and 14%, respectively). The average AUC and TSS values were 0.918 (SD = 0.002) and 0.727 (SD = 0.007), indicating a very good predictive power of the final model.

The averages of the selected current model identified areas with a high level of suitability across the temperate zone and southern part of eastern Europe, the western Caucasus and some regions in northern and eastern Kazakhstan and in the southern part of West Siberia (Fig. 4). The projection of the current ecological niches of N. n. scutata under Holocene and late Pleistocene climatic conditions enabled us to consistently analyse shifts of range boundaries in the studied period. During the LGM (~ 21 ka), the presumed range of this subspecies covered territories coinciding with the northern and southern coasts of the Black Sea, regions in southern Ukraine and south of European Russia, western and central regions of the North Caucasus and Western Transcaucasia (Fig. 4). No relatively suitable areas were found in Central Asia or Siberia. During the warmer Heinrich Stadial 1 (17.0–14.7 ka), living conditions for these snakes in the south of the Russian Plain, along the coasts of the Black Sea and the Caucasus improved dramatically (Fig. 4). During Greenland Interstade 1 (14.7–12.9 ka) and the Younger Dryas Stadial (12.9–11.7 ka), climatic conditions became suitable for penetration into the more northern regions of eastern Europe and West Siberia and then into eastern ­Kazakhstan (Fig. 4). Since the early Holocene (11.7–­8.3 ka), the climatic conditions for the snakes have become most suitable throughout their entire current range. Thus, during this period, the distribution range could have reached its maximum extension to the Baikal region and Mongolia. Nevertheless, to the end of the Holocene, the living conditions in the east of the range began to gradually to deteriorate, contributing to a fragmentation of the range. The SDM projection to the mid-Pleistocene (~787 ka) and the last interglacial (~130 ka) suggest a distribution range that resembles the current one (Fig. 4). Hypothetical colonisation routes of the mitochondrial lineage 8 of N. n. scutata, according to the phylogeographic and SDM evidence, are shown in Figure 5.

To date, the easternmost genotyped specimens of N. n. scutata have been found in upper reaches of the Emba River in western Kazakhstan (Guicking et al. 2006). Thus, it has been supposed that only the eastern part of the range is occupied by the same mtDNA lineage. Our data undoubtedly confirm the distribution of mitochondrial lineage 8 of N. n. scutata across Siberia to Transbaikalia. Additional sampling around the terra typica of this subspecies further confirmed the presence of only one phylogenetic lineage in this area. It is noteworthy that according to the SDM, climatic conditions in the north and east of Baikal Lake are poorly suitable for the life of snakes, thereby confirming that thermal springs are necessary for the existence of local populations (Litvinchuk et al. 2013).

Previous studies on mtDNA variability in this species have been based mainly on cyt b and ND4 sequences (Kindler et al. 2013, 2017; Asztalos et al. 2021; Jablonski et al. 2023). They have detected similar levels of variation and concordant patterns of subdivision both between and within different lineages of N. natrix, including lineage 8 of N. n. scutata. By contrast, the haplotype networks reconstructed with ND4 contain fewer mutations (Asztalos et al. 2021). Jablonski et al. (2023) identified an ND4 haplotype – at the northern edge of the Zagros Mountains (Turkey) – that turned out to be the hypothetical ancestral haplotype for the other ND4 haplotypes found in Transcaucasia and Turkey. Our newly obtained cyt b sequences are in agreement with the results from previous papers on both cyt b and ND4 markers and suggest that Transcaucasia and eastern Turkey are the centre of mtDNA diversity in N. n. scutata.

The observed mtDNA diversity of N. n. scutata in West Siberia proved to be even more uniform than in other reptile species with similar geographic distributions. For instance, in the sand lizard Lacerta agilis Linnaeus, 1758 several distinct cyt b haplotypes were found in West Siberia, although it was reported that the subspecies L. a. exigua Eichwald, 1831 colonised the eastern part of the Eastern European Plain, Siberia and Kazakhstan from a Caucasian refugium (Kalyabina-Hauf and Ananieva 2004). Notably, identical haplotypes of L. a. exigua were also registered in localities separated by more than 3000 km. Ursenbacher et al. (2006) found different cyt b and control region haplotypes even within only three specimens of the common adder V. berus from Siberia (the Altai Mountains and Krasnoyarsk Krai), while Cui et al. (2016) identified only a single widespread haplotype of V. berus in the Southern Altai Mountains (China, n = 9), but the latter specimens were sampled in an apparently isolated population at the very edge of the species range. In contrast, all 25 N. n. scutata individuals (from 13 localities in West Siberia) analysed in the present study share a single haplotype. Assuming high similarity in generation times, metabolic rates and other factors that could lead to possible differences in substitution rates between Vipera and Natrix species, we suppose that historical and ecological factors that are linked to population expansion, rather than differences in substitution rates, have played a leading role in the formation of the observed differences in spatial genetic diversity.

There are examples of other wide-ranging Palaearctic animals for which common mitochondrial haplotypes have been found at sites separated by thousands of kilometres. Identical cyt b haplotypes of the red squirrel Sciurus vulgaris Linnaeus, 1758 have been registered in Western Europe and Northeastern China (Grill et al. 2009; Liu et al. 2014). In a caudate amphibian – the Siberian salamander Salamandrella keyserlingii Dybowski, 1870 – a common haplotype was found between the Urals and the Kamchatka Peninsula (Malyarchuk et al. 2011). The moor frog R. arvalis has identical haplotypes throughout its range from Central Europe to East Siberia (Babik et al. 2004). A similar pattern has been observed in purely aquatic vertebrates, for example, identical cyt b haplotypes of the common roach Rutilus rutilus sensu lato have been recovered from western parts of Russia to the Lena River basin in East Siberia (Levin et al. 2017).

The results of neutrality tests, of the mismatch analysis are concordant with the distribution analysis of the haplotype frequencies and the observed network pattern, implying a relatively recent fast expansion of N. n. scutata across the vast territory of the Eastern European Plain and Siberia. The structure and geographical distribution of mtDNA genetic diversity within lineage 8 along with the results of species distribution modelling provide evidence that the North Caucasus, southern Ukraine and the south of the European part of Russia are the main refugia areas for the post-glacial colonisation of the Eastern European Plain and Siberia. A similar pattern of spatial mtDNA variation is well documented for a number of other Palaearctic species colonising northern areas from glacial refugia in Southern Europe and the Caucasus (for example, Ursenbacher et al. 2006; Joger et al. 2007; Grill et al. 2009; Lajbner et al. 2011). For instance, it has been shown that the steppe viper V. renardi has colonised most of its extant range from refugia in the North Caucasus, whereas two distinct lineages were found in southern Ukraine and Crimea (Zinenko et al. 2015).

In our study, haplotypes from the three different star-like clusters (S1, S2 and S3) were detected in the northern parts of the distribution range of N. n. scutata (i.e., areas that have been post-glacially colonised). Combined with the distinct geographical distributions of these groups, this observation might indicate northward dispersal of this subspecies from different areas (“microrefugia”) within the Caucasian refugium after the LGM. Most of the current range of this subspecies is occupied by haplogroup S1. The spatial distribution of the haplotypes and their diversity suggest that this haplogroup dispersed from the southern part of European Russia, most likely somewhere north of the Kuma–Manych depression. The SDM is also consistent with the existence of suitable conditions in this area during the LGM (Fig. 4). The last glacial refugium of haplogroup S3 was probably located in the area between the MCR and the Kuma–Manych depression. The large brackish/saline strait/bay that periodically existed in the depression during Caspian Sea transgressions (including the early Khvalynian transgression at the end of the Pleistocene; Fig. 1C) should have limited the northward dispersal of its ancestral haplotype to some extent. Conversely, the simultaneous regression of the Black Sea basin (with disappearance of the Azov Sea) may facilitate a north-westward dispersal from this area, as suggested by the discovery of this haplogroup in Lithuania.

Western Transcaucasia seems to be the main microrefugium for the S2 group, as do areas along the northern shore of the Black Sea. The retreat of the sea in the Late Pleistocene probably helped with the dispersal of this group in a north-westerly direction. Reconstruction of ancestral areas for the Caspian lineage of congeneric species N. tessellata by Jablonski et al. (2024) also points to several areas of post-LGM dispersal in this geographical region, along with a strong barrier role of the MCR. Alternatively, the pattern observed by us for N. n. scutata might have arisen due to expansion from a narrower area with the Caucasian refugium being located in the north of the Black Sea and/or Transcaucasia regions. Expansion might have happened earlier in the north-western direction, as suggested by the SDM (Fig. 4), thereby promoting somewhat higher haplotype diversity in this area. Denser sampling of these areas and an increase in sequence length per specimen are needed to achieve better discrimination between mitochondrial haplotypes.

No evidence of glacial refugia was revealed for Central Asia or Southern Siberia by our species distribution modelling, thus supporting the hypothesis that most of the current range was not suitable during the LGM. Nevertheless, two regions with distinct genetic structure relative to surrounding areas were identified, suggesting that some populations may have survived outside the main refugia during the LGM while retaining some genetic diversity, whereas concurrently, species disappeared completely from large areas around them. In Transbaikalia, the two analysed individuals possess unique haplotypes of haplogroups 1, whereas for the much more extensively sampled Siberian populations west of Lake Baikal no diversity was revealed. The most likely explanation for this pattern is that N. n. scutata colonised the area in a previous interglacial period (as evidenced by the SDM, Fig. 4). During the last glacial period, the survival of the species in this area may have been facilitated by the possible presence of thermal waters. This survival was followed by a rapid re-colonisation of West Siberian territories from the southern refugium during the Holocene, thereby resulting in uniform mtDNA diversity across thousands of kilometres.

The same scenario is worth considering for the Mugodzhar Range region in western Kazakhstan, where unique haplotypes of haplogroups S1 and S3 were detected (Fig. 1C). Nevertheless, this area may have been more accessible for grass snakes during the mid-Holocene (Fig. 4), and the presence of an S3 haplotype could be explained by the eastward dispersal of this group along the Caspian Sea shore and the Emba River. This appears to be the simplest explanation for the observed pattern, but the absence of S3 haplotypes along the Ural River does not favour this hypothesis. Therefore, until more specimens are analysed, both hypotheses are plausible.

In summary, the identified spatial patterns of ­mtDNA diversity provide evidence supporting the existence of cryptic northern refugia of N. n. scutata in Transbaikalia and less likely in the Mugodzhar Range. There is an obvious need for further sampling in these and surrounding areas to test these hypotheses within a model-based framework. Additional sampling is required in the south-western part of this subspecies’ range to clarify migration routes and the haplogroup distribution. Furthermore, a substantial increase in mtDNA sequence data per specimen is necessary to improve resolution of inferred relationships and enable reliable divergence time estimates. The easternmost populations of N. n. scutata in East Siberia and southeastern Kazakhstan have yet to be examined. Despite these limitations, this study offers new insights into (and intriguing hypotheses about) the phylogeographic history of reptiles in the northern Palaearctic.