Research Article |
Corresponding author: Evgeniy Simonov ( ev.simonov@gmail.com ) Academic editor: Uwe Fritz
© 2024 Evgeniy Simonov, Artem Lisachov, Spartak Litvinchuk, Anastasia Klenina, Polina Chernigova, Alexander Ruchin, Andrey Bakiev, Kazhmurat Akhmedenov.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Simonov E, Lisachov A, Litvinchuk S, Klenina A, Chernigova P, Ruchin A, Bakiev A, Akhmedenov K (2024) Recent recolonisation of West Siberia and northern cryptic refugia in the grass snake Natrix natrix scutata (Pallas, 1771) (Squamata: Natricidae). Vertebrate Zoology 74: 565-576. https://doi.org/10.3897/vz.74.e123485
|
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.
Last Glacial Maximum, Palaearctic, phylogeography, postglacial expansion, species distribution modelling
The cyclic climatic changes during the Quaternary led to iterative contraction and expansion of geographical ranges of species (
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 (
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 (
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.
We focus only on specimens harbouring the most widespread mtDNA lineage 8 (or ‘green’) of N. n. scutata, as defined by
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 (
The cyt b sequences available in GenBank (NCBI) and belonging to the N. n. scutata clade along with published information on haplotypes (
Geographic structure of cyt b variation in N. n. scutata was examined in BARRIER 2.2 (
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 (
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 (
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 (
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 (
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.
Most haplotypes formed a star-like cluster in the haplotype network (Fig.
Two other minor star-like groups of haplotypes could be identified in the haplotype network (Fig.
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
Several measures of genetic diversity were calculated for the population groups inferred by BARRIER (Table
Genetic diversity estimates of cyt b for star-like haplogroups and geographical groups of Natrix natrix scutata identified with BARRIER (see Fig.
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
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.
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.
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.
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.
To date, the easternmost genotyped specimens of N. n. scutata have been found in upper reaches of the Emba River in western Kazakhstan (
Previous studies on mtDNA variability in this species have been based mainly on cyt b and ND4 sequences (
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.
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 (
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,
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.
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
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.
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.
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.
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).
Tables S1–S3
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.