Research Article
Print
Research Article
Recent recolonisation of West Siberia and northern cryptic refugia in the grass snake Natrix natrix scutata (Pallas, 1771) (Squamata: Natricidae)
expand article infoEvgeniy Simonov§|, Artem Lisachov#, Spartak Litvinchuk¤«, Anastasia Klenina», Polina Chernigova˄, Alexander Ruchin˅, Andrey Bakiev», Kazhmurat Akhmedenov
‡ M. Utemisov West Kazakhstan University, Uralsk, Kazakhstan
§ A. N. Severtsov Institute of Ecology and Evolution, Russian Academy of Sciences, Moscow, Russia
| Western Caspian University, Baku, Azerbaijan
¶ Kasetsart University, Bangkok, Thailand
# Institute of Cytology and Genetics, Siberian Branch of RAS, Novosibirsk, Russia
¤ Dagestan State University, Makhachkala, Russia
« Institute of Cytology, RAS, St. Petersburg, Russia
» Institute of Ecology of Volga River Basin RAS, Togliatti, Russia
˄ Unaffiliated, Tyumen, Russia
˅ Joint Directorate of the Mordovia State Nature Reserve and National Park Smolny, Saransk, Russia
Open Access

Abstract

Most of the geographic range of the grass snake Natrix natrix is occupied by the subspecies N. n. scutata, which occurs from Eastern Europe to East Siberia. This study addressed the phylogeography of this subspecies via analysis of variation of partial cytochrome b (cyt b) sequences from 135 individuals sampled at 91 localities across its range, in addition to species distribution modelling. A haplotype network was constructed and identified a major star-like haplogroup that harbours most of the analysed specimens and is considered the main source of recolonisation of vast territories of northern Eurasia after the Last Glacial Maximum (LGM). Most of this subspecies’ current range is occupied by haplotypes from a single haplogroup, with probable refugia in the south of the European part of Russia or the North Caucasus. The most frequent (ancestral) haplotype in this group is the only one discovered in West Siberia among 25 specimens from 13 sampling sites, indicating relatively recent colonisation of Siberia. According to species distribution modelling, no relatively suitable areas were present in Central Asia or Siberia during the LGM. Nonetheless, the only two specimens examined from the easternmost area of the species’ geographic range, near Lake Baikal, had unique haplotypes that differed from the most common haplotype by one and two substitutions. The most probable explanation for this pattern is that N. n. scutata colonised the area during a previous interglacial period and survived here during the last glacial. The Mugodzhar Range in western Kazakhstan also showed cyt b differentiation as compared to surrounding areas but warrants further sampling to test competing hypotheses.

Keywords

Last Glacial Maximum, Palaearctic, phylogeography, postglacial expansion, species distribution modelling

Introduction

The cyclic climatic changes during the Quaternary led to iterative contraction and expansion of geographical ranges of species (Hewitt 1999; Stewart et al. 2010; Feliner 2011). These fluctuations dramatically shaped the population structure and distribution of many species in the temperate zone. Species’ modern distribution that we can observe in the Palaearctic and Nearctic has been caused by expansion of geographic ranges of thermophilic species and retraction of ranges of cold-adapted species after the Last Glacial Maximum (LGM, ~21 thousand years ago [ka]) (Shafer et al. 2010; Stewart et al. 2010). The phylogeographic patterns associated with Pleistocene glaciation cycles are especially well studied in Europe (reviewed by Hewitt, 1999; Joger et al. 2007; Stewart et al. 2010; Varga 2010). Along with well-documented southern refugia for temperate flora and fauna (Hewitt 1999; Joger et al. 2007; Feliner 2011), there is growing evidence for the existence of so-called cryptic northern refugia (Stewart and Lister 2001; Stewart et al. 2010; Tzedakis et al. 2013; Kindler et al. 2018). This type of glacial refugia was situated at higher latitudes as compared to the expected areas of suitable habitats to the south. One of the most prominent examples of such refugia are the Carpathians. Phylogenetic studies have revealed that they have served as an extra-Mediterranean glacial refugium for various species, including poikilotherms: the moor frog Rana arvalis Nilsson, 1842 (Babik et al. 2004), the yellow-bellied toad Bombina variegata (Linnaeus, 1758) (Fijarczyk et al. 2011), newts Triturus cristatus (Laurenti, 1768) (Wielstra et al. 2015) and Lissotriton vulgaris (Linnaeus, 1758) (Babik et al. 2005) and the snake Vipera berus (Linnaeus, 1758) (Ursenbacher et al. 2006). The integration of species distribution models (SDMs) based on bioclimatic variables with phylogeography provides further insights into species’ responses to warming–cooling cycles and assists with the identification of climatic refugia (e.g., Wielstra et al. 2015; Jablonski et al. 2024). Guillon et al. (2024) demonstrated that the use of physiologically meaningful climatic variables may enhance the accuracy of such inferences. Furthermore, those researchers proposed models of northern refugia for the cold-tolerant reptile species V. berus and Zootoca vivipara (Lichtenstein, 1823).

Although many Palaearctic species of European origin have wide distributions across Eurasia (from West/Central Europe to Siberia and shores of the Pacific Ocean), an analysis of their phylogeography is often limited to the specimens collected in Europe. The concept of cryptic northern refugia implies that some of these species might have survived in such refugia located in Central Asia, Southern Siberia or the Far East. High continentality of the climate in this area will cull many European species adapted to a less seasonably variable more humid climate, but some more cold-adapted species may have had a chance to survive. For instance, it has been shown that some relatively cold-tolerant reptiles – Vipera renardi (Christoph, 1861) and V. ursinii (­Bonaparte, 1835) – can survive close to the periglacial zone and even increased effective population size during the LGM (Zinenko et al. 2015). Given the overall lack of phylogeographical research in the eastern part of the Palaearctic, there is a clear need for such studies in this region. The species that are well investigated in this regard in Europe but rarely sampled in the Asian part of their ranges are a good choice to start closing this gap in our knowledge of Palaearctic phylogeography.

The grass snake Natrix natrix (Linnaeus, 1758) has a wide range in the Palaearctic: It extends from Central Europe to the south of Eastern Siberia and adjacent regions of northern Mongolia (Fritz and Schmidtler 2020). The phylogeny and phylogeography of the grass snake in Europe have been extensively characterised, by means of both mitochondrial DNA (mtDNA) sequences (cytochrome b = cyt b and ND4 genes) and nuclear DNA markers (microsatellites) (Guicking et al. 2006; Kindler et al. 2013, 2014, 2017, 2018; Fritz and Schmidtler 2020; Asztalos et al. 2021; Jablonski et al. 2023). These articles revealed incongruence between morphological delineation of subspecies of N. natrix and its genetic structuring (Kindler et al. 2013, 2014, 2017). Most of this species’ range is occupied by mtDNA ‘lineage 8’ (sensu Kindler et al. 2013) corresponding to the subspecies N. n. scutata (Pallas, 1771) (Fritz and Schmidtler 2020), which inhabits parts of eastern Europe, northern Turkey, Transcaucasia and the Urals and is supposed to be present in Siberia, eastern Kazakhstan, northern Kyrgyzstan and northern Mongolia. Moreover, Asztalos et al. (2021) demonstrated that there is no differentiation between the N. n. scutata primary lineage and mtDNA lineages 1 and 2 at 13 microsatellite loci. As a result, the geographic range of this subspecies has been extended to include eastern Transcaucasia and Iran. Some studies also revealed the Caucasian origin of this subspecies (Kindler et al. 2013), implying post-glacial northward recolonisation from this area. On the other hand, only a few specimens from the East European Plain were analysed there, with almost no DNA samples from Asia.

To provide the first insight into the post-glacial expansion history of N. n. scutata across most of its range, we aimed (i) to analyse mitochondrial genetic diversity based on a cyt b fragment in a phylogeographic context, (ii) to model species distributions for current and past conditions on the basis of bioclimatic variables and (iii) to evaluate the possibility of cryptic northern refugia for this subspecies using combined evidence from phylogeography and an SDM.

Materials and Methods

Sampling, DNA extraction and sequencing

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).

Figure 1. 

A Approximate distribution of Natrix natrix (green) and N. n. scutata (lineage 8) and sampling localities of the N. n. scutata used in this study. Colours for haplogroups and haplotypes correspond to the geographical groups in panel B. B TCS haplotype network of N. n. scutata using cyt b sequences. Specimens marked with an asterisk in panel A were excluded from this analysis because of missing data. C Spatial distribution of haplotype groups (positions of overlapping points were randomly adjusted to improve representation).

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 MgCl2), 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.

Phylogeographic analyses

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 (hP) 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) FS 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.

Species Distribution Modelling (SDM)

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 10th 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).

Results

MtDNA variability and phylogeography

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).

Figure 2. 

Heatmap of interpolated nucleotide diversity (Pi) and first 10 barriers identified by Monmonier’s algorithm in BARRIER within the approximate range of lineage 8 of Natrix natrix scutata.

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).

Table 1.

Genetic diversity estimates of cyt b for star-like haplogroups and geographical groups of Natrix natrix scutata identified with BARRIER (see Fig. 2).

Group n H hP Hd±SD Pi±SD K Tajima’s D Fu’s FS
Star-like haplogroups vs. all other haplotypes
S1 94 18 n/a 0.41±0.07 0.0009±0.0002 0.59 –2.478** –23.5***
S2 13 5 n/a 0.53±0.16 0.0010±0.0003 0.62 –1.775 –3.1
S3 6 3 n/a 0.60±0.22 0.0010±0.0004 0.67 –1.132 –0.9
Other 19 11 n/a 0.87±0.07 0.0052±0.0009 3.45 –1.101 –3.5
Groups with one barrier
Northern 103 23 22 0.51±0.06 0.0012±0.0002 0.77 –2.497*** –30.8***
Southern 31 16 15 0.90±0.04 0.0040±0.0006 2.59 –1.566 –9.2***
Groups with 10 barriers
EEP 65 15 12 0.52±0.08 0.0012±0.0002 0.76 N/A N/A
MUG 2 2 2 1.00±0.50 0.0046±0.0023 3.00 N/A N/A
ENC 7 4 2 0.81±0.13 0.0017±.00004 1.14 N/A N/A
LC 5 4 3 0.90±0.16 0.0052±0.0022 3.40 N/A N/A
NBS 2 2 2 1.00±0.50 0.0076±0.0038 5.00 N/A N/A
TB 2 2 2 1.00±0.50 0.0046±0.0023 3.00 N/A N/A
TC 23 10 8 0.87±0.05 0.0028±0.0005 1.82 N/A N/A
WS 25 1 0 N/A N/A N/A N/A N/A
WT 3 3 3 1.00±0.27 0.0030±0.0012 2.00 N/A N/A
Overall 134 38 N/A 0.71±0.04 0.0024±0.0003 1.54 –2.394** –34.2***
n – sample size; H – number of haplotypes; hP – number of private haplotypes; Hd – haplotype diversity; Pi – nucleotide diversity (per site); K – the mean number of nucleotide substitutions; Tajima’s D and Fu’s FS test results are also presented; *** P < 0.001, ** P < 0.01. Acronyms for specimens’ geographical groups designated by Monmonier’s algorithm as having 10 barriers: Northern – all samples north of the MCR; Southern – all samples south of the MCR; EEP – East European Plain and Caspian Depression (west Kazakhstan); MUG – Mugodzhar Range, upper reaches of the Emba river (west Kazakhstan); WS – West Siberia; TB – Transbaikalia; ENC – eastern North Caucasus; NBS – North Black Sea region; TC – Transcaucasia; LC – Lesser Caucasus; WT – western Turkey.

The results of Tajima’s D test and Fu’s FS 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.

Figure 3. 

Mismatch distribution analysis for the cyt b sequences of two major geographic groups (to the north and south of the MCR) of Natrix natrix scutata. Lines with dots are expected frequencies of pairwise differences under a model of sudden demographic expansion, while the grey bars are observed values.

Species Distribution Modelling (SDM)

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.

Figure 4. 

Predicted regions suitable for Natrix natrix scutata under mid-Pleistocene (~787 ka), last interglacial (~130 ka), LGM (~21 ka), latest Pleistocene (17.0–11.7 ka) and Holocene (11.7–0.0 ka) climatic conditions. The intensity of the blue colour reflects the degree of suitability (darker means more suitable). Territories covered by glaciers are grey. The red curves outline the limits of inhabited sites employed in SDM.

Figure 5. 

Hypothetical post-LGM northward colonisation routes of mitochondrial lineage 8 of Natrix natrix scutata, according to the phylogeographic evidence and SDM. Pink dashed arrows indicate possible pre-LGM events of migration into hypothesised northern refugia, while the black dashed arrow denotes the hypothetical route in the unsampled area.

Discussion

Distribution of Natrix natrix scutata

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.

Low genetic variation in West Siberia

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).

Location of glacial refugia and possible dispersal routes

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.

Acknowledgements

We thank Anatoliy Chulisov, Andrey Kiselev and Boris Levin for providing tissue samples and Anton Svinin for help in the field. We would like to extend our gratitude to Sylvain Ursenbacher and an anonymous reviewer for their valuable feedback and suggestions, which have helped us to improve the manuscript. This research was funded by the Science Committee of the Ministry of Science and Higher Education of the Republic of Kazakhstan (grant No. AP 19675960).

References

  • Aljanabi S, Martinez I (1997) Universal and rapid salt-extraction of high quality genomic DNA for PCR-based techniques. Nucleic Acids Research 25: 4692–4693. https://doi.org/10.1093/nar/25.22.4692
  • Allouche O, Tsoar A, Kadmon R (2006) Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS). Journal of Applied Ecology 43: 1223–1232. https://doi.org/10.1111/j.1365-2664.2006.01214.x
  • Asztalos M, Ayaz D, Bayrakcı Y, Afsar M, Tok CV, Kindler C, Jablonski D, Fritz U (2021) It takes two to tango – Phylogeography, taxonomy and hybridization in grass snakes and dice snakes (Serpentes: Natricidae: Natrix natrix, N. tessellata). Vertebrate Zoology 71: 813–834. https://doi.org/10.3897/vz.71.e76453
  • Babik W, Branicki W, Crnobrnja-Isailovic J, Cogălniceanu D, Sas I, Olgun K, Poyarkov NA, Garcia-París M, Arntzen JW (2005) Phylogeography of the two European newt species – discordance between mtDNA and morphology. Molecular Ecology 14: 2475–2491. https://doi.org/10.1111/j.1365-294X.2005.02605.x
  • Cobos ME, Peterson AT, Barve N, Osorio-Olvera L (2019) kuenm: An R package for detailed development of ecological niche models using Maxent. PeerJ 7: e6281. https://doi.org/10.7717/peerj.6281
  • Cui S, Luo X, Chen D, Sun J, Chu H, Li C, Jiang Z (2016) The adder (Vipera berus) in Southern Altay Mountains: Population characteristics, distribution, morphology and phylogenetic position. PeerJ 4: e2342. https://doi.org/10.7717/peerj.2342
  • Fijarczyk A, Nadachowska K, Hofman S, Litvinchuk SN, Babik W, Stuglik M, Gollmann G, Choleva L, Cogălniceanu D, Vukov T, Džukić G, Szymura JM (2011) Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories. Molecular Ecology 20: 3381–3398. https://doi.org/10.1111/j.1365-294­X.2011.05175.x
  • Fritz U, Schmidtler JF (2020) The Fifth Labour of Heracles: Cleaning the Linnean stable of names for grass snakes (Natrix astreptophora, N. helvetica, N. natrix sensu stricto). Vertebrate Zoology 70: 621–665. https://doi.org/10.26049/VZ70-4-2020-07
  • Gouy M, Guindon S, Gascuel O (2010) SeaView version 4: A multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Molecular Biology and Evolution 27: 221–224. https://doi.org/10.1093/molbev/msp259
  • Guillon M, Martínez-Freiría F, Lucchini N, Ursenbacher S, Surget-Groba Y, Kageyama M, Lagarde F, Cubizolle H, Lourdais O (2024) Inferring current and Last Glacial Maximum distributions are improved by physiology-relevant climatic variables in cold-adapted ectotherms. Journal of Biogeography: 1–16. https://doi.org/10.1111/jbi.14828
  • Harpending HC (1994) Signature of ancient population-growth in a low-resolution mitochondrial-DNA mismatch distribution. Human Biology 66: 591–600.
  • Jablonski D, Asztalos M, Yılmaz C, Avcı A (2023) The range-wide mitochondrial lineage of Natrix natrix scutata (Pallas, 1771) presented in the northern Zagros Mountains. Evolutionary Systematics 7: 67–71. https://doi.org/10.3897/evolsyst.7.89662
  • Jablonski D, Mebert K, Masroor R, Simonov E, Kukushkin O, Abduraupov T, Hofmann S (2024) The Silk roads: Phylogeography of Central Asian dice snakes (Serpentes: Natricidae) shaped by rivers in deserts and mountain valleys. Current Zoology 70: 150–162. https://doi.org/10.1093/cz/zoad008
  • Joger U, Fritz U, Guicking D, Kalyabina-Hauf S, Nagy ZT, Wink M (2007) Phylogeography of western Palaearctic reptiles – Spatial and temporal speciation patterns. Zoologischer Anzeiger 246: 293–313. https://doi.org/10.1016/j.jcz.2007.09.002
  • Kalyabina-Hauf SA, Ananjeva NB [Калябина-Хауф CA, Ананьева НБ] (2004) Филогеография и внутривидовая структура широ­коареального вида ящериц Lacerta agilis L., 1758 (Lacertidae, Sauria, Reptilia) (опыт использования митохондриального гена цитохрома b). Труды Зоологического Института РАН 302: 1–108 [with English summary].
  • Kindler C, Böhme W, Corti C, Gvoždík V, Jablonski D, Jandzik D, Metallinou M, Široký P, Fritz U (2013) Mitochondrial phylogeography, contact zones and taxonomy of grass snakes (Natrix natrix, N. megalocephala). Zoologica Scripta 42: 458–472. https://doi.org/10.1111/zsc.12018
  • Kindler C, Bringsøe H, Fritz U (2014) Phylogeography of grass snakes (Natrix natrix) all around the Baltic Sea: Implications for the Holocene colonization of Fennoscandia. Amphibia-Reptilia 35: 413–424. https://doi.org/10.1163/15685381-00002962
  • Kindler C, Chèvre M, Ursenbacher S, Böhme W, Hille A, Jablonski D, Vamberger M, Fritz U (2017) Hybridization patterns in two contact zones of grass snakes reveal a new Central European snake species. Scientific Reports 7: 7378. https://doi.org/10.1038/s41598-017-07­8­47-9
  • Levin BA, Simonov EP, Ermakov OA, Levina MA, Interesova EA, Kovalchuk OM, Malinina YA, Mamilov NS, Mustafayev NJ, Pilin DV, Pozdeev IV, Prostakov NI, Roubenyan HR, Titov SV, Vekhov DA (2017) Phylogeny and phylogeography of the roaches, genus Rutilus (Cyprinidae), at the eastern part of its range as inferred from mtDNA analysis. Hydrobiologia 788: 33–46. https://doi.org/­10.1007/s10750-016-2984-3
  • Liu Z, Li B, Ma J, Zheng D, Xu Y (2014) Phylogeography and genetic diversity of the red squirrel (Sciurus vulgaris) in China: Implications for the species’ postglacial expansion history. Mammalian Biology 79: 247–253. https://doi.org/10.1016/j.mambio.2014.02.004
  • Malyarchuk BA, Berman DI, Derenko MV, Perkova MA, Bulakhova NA, Grzybowski T, Leirikh AN (2011) Polymorphism of the mitochondrial cytochrome b gene, phylogeography, and molecular adaptation of the Siberian salamander (Salamandrella keyserlingii, Amphibia, Caudata). Biology Bulletin Reviews 1: 366–380. https://doi.org/10.1134/S2079086411040049
  • Manni F, Guerard E, Heyer E (2004) Geographic patterns of (genetic, morphologic, linguistic) variation: How barriers can be detected by using Monmonier’s algorithm. Human Biology 76: 173190.
  • Miller MP (2005) Alleles In Space (AIS): computer software for the joint analysis of interindividual spatial and genetic information. Journal of Heredity 96: 722–724. https://doi.org/10.1093/jhered/esi119
  • Monmonier MS (1973) Maximum-difference barriers: An alternative numerical regiodization method. Geographical Analysis 3: 245–261.
  • Phillips SJ, Anderson RP, Dudík M, Schapire RE, Blair ME (2017) Opening the black box: An open-source release of Maxent. Ecography 40: 887893. https://doi.org/10.1111/ecog.03049
  • QGIS Development Team (2023) QGIS Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.osgeo.org
  • Sambrook J, Fritsch E, Maniatis T (1989) Molecular Cloning: A Laboratory Manual. Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY, xxviii + 1546 pp.
  • Shafer ABA, Cullingham CI, Côté SD, Coltman DW (2010) Of glaciers and refugia: A decade of study sheds new light on the phylogeography of northwestern North America. Molecular Ecology 19: 4589–4621. https://doi.org/10.1111/j.1365-294X.2010.04828.x
  • Stewart JR, Lister AM, Barnes I, Dalén L (2010) Refugia revisited: Individualistic responses of species in space and time. Proceedings of the Royal Society B 277: 661–671. https://doi.org/10.1098/rspb.2009.1272
  • Ursenbacher S, Carlsson M, Helfer V, Tegelström H, Fumagalli L (2006) Phylogeography and Pleistocene refugia of the adder (Vipera berus) as inferred from mitochondrial DNA sequence data. Molecular Ecology 15: 3425–3437. https://doi.org/10.1111/j.1365-294X.2006.03031.x
  • Varga Z (2010) Extra-Mediterranean refugia, post-glacial vegetation history and area dynamics in eastern Central Europe. In: Habel JC, Assmann T (Eds) Relict Species: Phylogeography and Conservation Biology. Springer, Heidelberg, Dordrecht, London, New York, 57–87. https://doi.org/10.1007/978-3-540-92160-8
  • Wielstra B, Babik W, Arntzen JW (2015) The crested newt Triturus cristatus recolonized temperate Eurasia from an extra-Mediterranean glacial refugium. Biological Journal of the Linnean Society 114: 574–587. https://doi.org/10.1111/bij.12446
  • Зима ЮА, Федоренко ВА [Zima YuA, Fedorenko VA] (2019) О новых находках амфибий и рептилий в Акмолинской области. Sele­vinia 27: 51–60.
  • Zinenko O, Stümpel N, Mazanaeva L, Bakiev A, Shiryaev K, Pavlov A, Kotenko T, Kukushkin O, Chikin Y, Duisebayeva T, Nilson G, Orlov NL, Tuniyev S, Ananjeva NB, Murphy RW, Joger U (2015) Mitochondrial phylogeny shows multiple independent ecological transitions and northern dispersion despite of Pleistocene glaciations in meadow and steppe vipers (Vipera ursinii and Vipera renardi). Molecular Phylogenetics and Evolution 84: 85–100. https://doi.org/10.1016/j.ympev.2014.12.005

Supplementary material

Supplementary material 1 

Tables S1–S3

Simonov E, Lisachov A, Litvinchuk S, Klenina A, Chernigova P, Ruchin A, Bakiev A, Akhmedenov K (2024)

Data type: .xlsx

Explanation notes: Table S1. Natrix natrix scutata samples used in the study. — Table S2. Presence localities of Natrix natrix scutata used for species distribution modelling with MAXENT. — Table S3. Performance metrics for parameter settings used for creating the present time species distribution model and relative contribution (%) of variables.

This dataset is made available under the Open Database License (http://opendatacommons.org/­licenses/odbl/1.0). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (148.33 kb)
login to comment