Data for 683 grass snakes (Natrix natrix) and 332 dice snakes (N. tessellata) were analyzed for the present study. A total of 181 tissue samples (115 N. natrix, 66 N. tessellata) were newly processed and microsatellite data were generated for 15 additional samples of N. natrix with previously published mtDNA sequences (Kindler et al. 2013, 2017). The remaining data originate from earlier studies (Guicking et al. 2006, 2009; Kindler et al. 2013, 2017; Kyriazi et al. 2013; Asztalos et al. 2021a). Among our material were five samples of big-headed grass snakes from the putative distribution range of N. megalocephala that were, according to their morphology, tentatively identified with this nominal species (Table S1).

Two marker systems (mtDNA, microsatellite loci) were applied, which were successfully used in our previous studies on the phylogeography and hybridization of grass snakes (Kindler et al. 2013, 2017, 2018a, b; Pokrant et al. 2016; Kindler and Fritz 2018; Schultze et al. 2019, 2020; Asztalos et al. 2020, 2021a, b; Ahnelt et al. 2021).

Laboratory procedures followed Kindler et al. (2013) and Pokrant et al. (2016). Samples collected within the frame of TÜBİTAK project 116Z359 in Turkey were processed in the Biodiversity Research Laboratory of Ege University, İzmir. The remaining new material was studied in the molecular laboratory of Senckenberg Dresden.

Two mtDNA fragments were sequenced, the cytochrome b (cyt b) gene (1,117 bp) and the partial ND4 gene plus adjacent DNA coding for tRNAs (tRNA-His, tRNA-Ser, tRNA-Leu; 866 bp). For some samples the obtained sequences were shorter and for a few samples, only one mtDNA fragment could be sequenced (Table S1).

In addition to mtDNA sequencing, the samples were also genotyped at the same 13 polymorphic microsatellite loci as in our previous studies. For N. tessellata, the presence of the microsatellite motifs was verified by sequencing the PCR products in both directions using unlabeled primers. For detailed information on microsatellite loci, primers and PCR conditions for microsatellites and mtDNA, see Tables S2 and S3, Kindler et al. (2013), and Pokrant et al. (2016). Allele lengths of material processed in İzmir were calibrated using 50 samples from the Balkans that were studied in both laboratories. For a few of the genotyped samples, no mtDNA sequences could be obtained (Table S1).

Cyt b sequences could be generated for 110 samples (55 Natrix natrix, 55 N. tessellata), ND4+tRNAs for 112 samples (73 N. natrix, 39 N. tessellata). These sequences were aligned with those from previous studies (Guicking et al. 2006, 2009; Kindler et al. 2013, 2014, 2017, 2018a; Kyriazi et al. 2013; Kindler and Fritz 2018; Schultze et al. 2019, 2020; Asztalos et al. 2021a) using bioedit 7.0.5.2 (Hall 1999). To identify mitochondrial haplotypes, previously characterized haplotypes of the two species (Guicking et al. 2006, 2009; Kindler et al. 2017, 2018a; Schultze et al. 2019, 2020; Ahnelt et al. 2021; Asztalos et al. 2021a, b) were added to the alignments of each mtDNA block, resulting in a 1,117-bp-long alignment of 705 cyt b sequences for N. natrix and 435 cyt b sequences for N. tessellata and an 866-bp-long alignment of 722 ND4+tRNAs sequences for N. natrix and 38 ND4+tRNAs sequences for N. tessellata. The alignment for ND4+tRNAs of N. tessellata included only the sequences generated for the present study, whereas the other alignments also contained all previously published homologous sequence data. For each alignment, parsimony networks were drawn using tcs 1.21 (Clement et al. 2000), with gaps coded as fifth character state. Using the default 95% connection limit, unconnected haplotype clusters were obtained for the different genetic lineages within each species, which is why the connection limit was arbitrarily set to 150 steps. Haplotype terminology for N. natrix follows Kindler et al. (2013) and for the cyt b of N. tessellata, Guicking et al. (2009). For the ND4+tRNAs fragment, which was not studied by Guicking et al. (2009), an analogous system is introduced here in that each haplotype is consecutively numbered and this number follows the code letter for the haplotype cluster (e.g., E1, E2, etc.).

For each of the two species and for individual lineages within each species, uncorrected p distances (averages) were computed based on the haplotypes of each mtDNA fragment using mega 10.2.3 (Kumar et al. 2018) and the pairwise deletion option.

Our microsatellite data were subjected to unsupervised Bayesian cluster analyses using structure 2.3.4 (Pritchard et al. 2000; Falush et al. 2003). micro-checker 2.2.3 (van Oosterhout et al. 2004) revealed no evidence for null alleles, necessitating no corrections for null alleles. structure was run using the admixture model and correlated allele frequencies. All calculations were repeated 10 times for each K ranging from 1 to 10, using a Monte Carlo Markov chain of 1,000,000 generations and a burn-in of 250,000. To infer the most likely number of clusters (K), structure harvester (Earl and vonHoldt 2012) was used, a software that implements the ΔK method of Evanno et al. (2005). Population structuring and individual admixture were visualized using the R package pophelper 2.2.9 (Francis 2017). Since structure is known to detect only the uppermost level of population differentiation (Evanno et al. 2005), different hierarchical runs were executed: A first calculation included all data from 639 Natrix natrix and 66 N. tessellata. Then, data for all individuals with evidence for interspecific hybridization were removed and a pruned dataset for each species alone was processed. Within each species, further analyses were run for lower hierarchical levels to examine for finer substructuring. For each of these analyses, data for snakes with evidence for admixture between the relevant interspecific clusters were excluded (Table S1).

To infer whether individual genotypes in structure analyses represent hybrids or pure individuals, thresholds for hybrid identity were determined using hybridlab 1.0 (Nielsen et al. 2006). For doing so, from the first structure run 25 pure individuals each of N. natrix and of N. tessellata were selected as parental genotypes and the cluster membership coefficients (Q values) from the four clusters corresponding to N. natrix were summed up and opposed to the values of N. tessellata (for the rationale, see under Results and Discussion). Then, 25 genotypes of each hybrid class (F₁, F₂, two backcrosses) were modeled and used for inferring the thresholds for pure genotypes of each parental species (Table S4). The same approach was used for finding thresholds for pure clusters within N. natrix (Table S5). Inferred hybrid thresholds are summarized in the Supplementary Tables S4 and S5. For the lowermost hierarchical analysis for N. natrix and for N. tessellata, the ΔK method suggested the presence of three and four distinct clusters, respectively. Since hybridlab cannot infer hybrid thresholds for more than two clusters, an arbitrary threshold was applied here and genotypes with proportions of cluster membership below 80% were treated as admixed (Randi 2008).

structure analyses are based on population-genetic presumptions (Hardy-Weinberg equilibrium, linkage equilibrium; Pritchard et al. 2000) and therefore prone to bias from uneven sample sizes (Puechmaille 2016). To cross-examine structure results for potential artifacts, additional Principal Component Analyses (PCAs) were run using the R package adegenet 2.1.1 (Jombart 2008). PCAs are independent from any population-genetic presumptions and exclusively based on genetic information, rendering them less sensitive to different sample sizes. PCAs were computed for the same datasets as used for structure analyses, and the individual results were depicted either according to their morphological or mitochondrial and nuclear genomic identity and compared to our structure results.

For selected populations, allelic richness was calculated using fstat 2.9.4 (Goudet 1995). The number of alleles (n_A), the average number of alleles per locus (n_Ā), average observed and expected heterozygosity (H_O and H_E) and the inbreeding coefficient (F_IS) were inferred using arlequin 3.5.2.2 (Excoffier and Lischer 2010).

In parsimony network analyses (cyt b, ND4+tRNAs), each mitochondrial lineage of Natrix natrix and N. tessellata revealed by phylogenetic analyses (Guicking et al. 2009; Guicking and Joger 2011; Kindler et al. 2013) corresponded to a distinct haplotype cluster.

Previous investigations focused either on the phylogeography of N. natrix (e.g., Kindler et al. 2013, 2017) or N. tessellata (Guicking et al. 2009; Guicking and Joger 2011), but the two codistributed species were never examined together. Consequently, these studies were unable to find signals for interspecific hybridization. However, there is a growing body of evidence that interspecific hybridization in animals is more frequent than traditionally thought and may contribute to the enrichment of the gene pools of the involved taxa, in particular via adaptive introgression (Taylor and Larson 2019; Pfennig 2021). Thus, it comes not entirely unexpected that our study revealed evidence for hybridization when using samples of both species, irrespective of whether the microsatellite data were analyzed with structure or PCAs (Figs 4 and 5). Why we found interspecific hybridization largely restricted to Turkey remains unclear and calls for further research.

When the putative hybrids from Turkey are compared to the parental species (Table S10), it is obvious that the hybrids share with each parental species many otherwise species-specific alleles. This corroborates their hybrid status, as do specimens that are assigned in the PCAs to the N. natrix cluster although they were morphologically identified as N. tessellata (Fig. 5; Table S1). Also, three admixed snakes that were morphologically determined either as N. natrix (ZDEU DNA1280) or N. tessellata (ZDEU DNA1269, DNA1276) harbored a mitochondrial haplotype of the other species (Table S1), supportive of hybridization. Furthermore, our data indicate that backcrosses occurred (Fig. 4; Table S1), i.e., that the hybrids are obviously not sterile.

Yet, our study also showed the limitations of microsatellites. Using these markers, only recent hybrid and backcross generations can be detected (Sanz et al. 2009). Thus, we cannot trace back older or ancient interspecific gene flow and its extent. Also, a caveat is the possible homoplasy of microsatellites in different species (Pujolar et al. 2014), so that our results need to be scrutinized using other marker systems. Currently, a follow-up investigation using whole genomes is underway. It will also allow the examination of older hybridization events and clarify which parts of the genomes have been exchanged. In any case, our present study exemplifies that phylogeographic investigations disregarding codistributed congeners might miss evidence for interspecific hybridization due to the study design.

Regarding Cypriot grass snakes, our study revealed a known weakness of the ∆K method and the structure software: The misidentification of K = 2 as top level of the hierarchical structure, even when more clusters should be present (Puechmaille 2016; Janes et al. 2017). This flaw is further exaggerated by uneven sampling (Puechmaille 2016), which is frequently unavoidable. In our case, this situation led to a grave misplacement of most grass snakes from Cyprus. Under K = 2, they were lumped together with dice snakes, and only the inspection of K = 5 (Fig. 4) and the parallel application of PCAs (Fig. 5) unveiled that this assignment was an artifact. There is only a single allele of the 3TS locus exclusively shared between Cypriot grass snakes and N. tessellata (Table S11), so that neither hybridization nor homoplasy can be invoked as reason for the misplacement of Cypriot N. natrix. Consequently, the structure result under K = 2 is clearly erroneous.

Cypriot grass snakes are showing an excess of homozygotes, indicating inbreeding (Table S9). This could have caused their distinctiveness in structure analyses under K = 5 (Fig. 4) and in PCAs (Fig. 5). Also, mtDNA corroborates that Cypriot grass snakes are genetically impoverished. For cyt b and ND4+tRNAs only two haplotypes each were recorded on Cyprus (cyt b: gy3 in 18 individuals, gy9 in 6 individuals; ND4+tRNAs: gy6 in 6 individuals, gy11 in 17 individuals), whereas in other sites a greater diversity occurred. For instance on Gökçeada island, being represented by a similar sample size (n = 17), four haplotypes each were recorded both for cyt b and ND4+tRNAs and in the Karamık Swamp (n = 10), four haplotypes for cyt b and five for ND4+tRNAs (Table S1). Grass snakes from Gökçeada and the Karamık Swamp share their mitochondrial lineage with Cypriot N. natrix (Fig. 3: lineage 7).

Grass snakes are extremely rare on Cyprus and seriously endangered (Blosat 2008; Baier and Wiedl 2010; Zotos et al. 2021). In light of this situation, Cypriot N. natrix are best regarded as genetically impoverished, and their distinct microsatellite profile does not imply deep genetic or taxonomic divergence.

The herpetofauna of Cyprus is characterized both by old endemics and recent invaders (Böhme and Wiedl 1994; Poulakakis et al. 2013). The most prominent case of an old endemic species is perhaps the little-known Cyprus racer, Hierophis cypriensis (Utiger and Schätti 2004; Kyriazi et al. 2013). Also, the genetic divergences of a toad (Bufotes cypriensis, Dufresnes et al. 2019), four lizards (Ablepharus budaki, Ophisops elegans, Poulakakis et al. 2013; Mediodactylus orientalis, Kotsakiozi et al. 2018; Phoenicolacerta troodica, Tamar et al. 2015) and two snakes (Telescopus fallax, Šmíd et al. 2019; Xerotyphlops vermicularis, Kornilios 2017) qualify them as old Cypriot endemics, whereas the shallow genetic differences of some other amphibians (Hyla savignyi, Pelophylax sp., Poulakakis et al. 2013) and reptiles (Mauremys rivulata, Vamberger et al. 2017; Acanthodacytlus schreiberi, Poulakakis et al. 2013; Chalcides ocellatus, Bisbal-Chinesta et al. 2020; Hemidactylus turcicus, Rato et al. 2011; Rhynchocalamus melanocephalus, Tamar et al. 2020) suggest that these taxa are recent invaders, even though this view has been challenged for the water frogs (Pelophylax sp., Plötner et al. 2012; but see Speybroeck et al. 2020). The immigration of the recent invaders seems to be frequently human-mediated (Böhme and Wiedl 1994; Poulakakis et al. 2013).

The local grass snakes have been treated as another candidate for an old Cypriot endemic (Böhme and Wiedl 1994). However, Cypriot grass snakes share their mitochondrial lineage with N. natrix from the southeastern Balkan Peninsula and western and southern Anatolia (Fig. 3: lineage 7). When the individual haplotypes are inspected, it becomes obvious that three of the four haplotypes found on Cyprus also occur in Bulgaria, Greece, and Turkey (cyt b: gy9; ND4+tRNAs: gy6, gy11; Table S1). If the grass snakes would represent an old Cypriot endemic, we would expect a deep mitochondrial separation from the mainland populations. This is not the case and the present evidence instead supports recent immigration and subsequent genetic impoverishment.

In general, the mitochondrial phylogeographies of N. natrix and N. tessellata show some similarities (Fig. 3), both with respect to the distribution of some mitochondrial lineages and the location of contact zones, suggestive of a shared biogeographic history. For both species, the Balkans, Anatolia or Western Asia in general, and the Caucasus have undoubtedly served as glacial refugia, where genetic lineages survived that diverged prior to the Pleistocene. For N. natrix, a fossil-calibrated molecular clock estimated that six lineages separated 3.4–7.2 million years ago, long before the Pleistocene glaciations, and only two of these lineages diverged then further during the early Pleistocene (1.5 and 1.7 million years ago; Kindler et al. 2018a). Yet, the current distribution pattern undoubtedly has been massively shaped by repeated glacial range restrictions and interglacial range expansions, and the wide distribution of one mitochondrial lineage in the northeast of the distribution ranges of the two species (Fig. 3: lineage 8 in N. natrix; lineage A in N. tessellata) clearly results from Holocene range expansions. Thus, the general phylogeographic patterns in both species match the well-known paradigm of southern genetic richness to northern purity (Hewitt 2000), reflecting the survival of distinct lineages in former southern glacial refugia and the Holocene colonization of more northerly regions by only one or a few of these lineages.

In contrast to N. natrix, N. tessellata shows additional lineages in regions where grass snakes do not occur, i.e., on the island of Crete, in the southern Levant and Egypt or in southern Central Asia. It remains a challenge to study the disjunct populations of N. natrix from Siberia and the Lake Baikal region to find out how they match the general pattern.

As in many other taxa (Joger et al. 2007; Schmitt 2020), multiple lineages, which most likely can be traced back to distinct glacial microrefugia, occur in the Balkan Peninsula and in the Caucasus region, both in N. natrix and N. tessellata. In these two regions, and further north in the Balkan Peninsula and in Anatolia, also lie secondary contact zones that can be understood as having formed during Holocene range expansions.

In dice snakes, we are lacking nuclear genomic evidence for most lineages. Our sample is restricted to Turkey. Yet, it shows that the two mitochondrial lineages occurring there (Fig. 3: E, T) match two microsatellite clusters each, i.e., that mitochondrial divergence parallels nuclear genomic differentiation. Further, there is clear evidence of admixture (Figs 8 and 9), suggestive of conspecificity. However, the deep mitochondrial divergences among all mitochondrial lineages of N. tessellata (Tables S7, S8) resemble those observed between N. astreptophora, N. helvetica, and N. natrix (Kindler et al. 2017, 2018a), i.e., three species with much restricted gene flow. This implies that dice snakes could also represent more than only one species. For a better understanding of their taxonomy, microsatellite data or other nuclear genomic information is needed beyond Turkey.

For grass snakes, range-wide microsatellite data are available (Figs 6 and 7). As in N. tessellata, there is evidence for admixture in the contact zones of distinct microsatellite clusters. Contrary to N. tessellata, in N. natrix the number of mitochondrial lineages, however, exceeds the number of microsatellite clusters (Fig. 3). Eight mitochondrial lineages correspond to only four microsatellite clusters, and if the Cyprus cluster is disregarded, only to three (Figs 3, 6 and 7). The uppermost hierarchical level of nuclear genomic divergence is represented by one microsatellite cluster from the west of our study region that matches mtDNA lineage 4. Another southern and eastern microsatellite cluster corresponds to several distinct mtDNA lineages (Figs 6A and 7A). This latter southern and eastern microsatellite cluster is further structured in three clusters when analyzed alone: Two western clusters (Cyprus versus southern Balkans plus western Turkey) largely match mtDNA lineages 3, 5, and 7, and a third eastern cluster corresponds to three other mtDNA lineages (1, 2, 8). This third cluster also includes all four samples morphologically identified as N. megalocephala (Fig. 6B). Further sampling from the southern Balkans and the eastern distribution range is needed to disentangle this situation in more detail. We cannot rule out that more samples will unveil more structure than our present data set. However, the PCA for the grass snakes from the south and east (Fig. 7B) shows quite compact scatter plots that rather argue for genetic homogeneity of each of the three clusters (with the inbred Cypriot population being most distinct).

How can this be explained? A similar situation is known for the closely related barred grass snake (N. helvetica). Grass snakes are highly mobile animals, with males having larger home ranges than females (Madsen 1984; Wisler et al. 2008; Reading and Jofré 2009; Meister et al. 2012). In Italy, five in part deeply divergent mtDNA lineages of N. helvetica correspond to only a single nuclear genomic cluster, due to largely male-mediated gene flow and female philopatry (Schultze et al. 2020). We hypothesize that the situation in our study region parallels that in Italy, i.e., that the distinct mtDNA lineages corresponding to one and the same microsatellite cluster are vestiges of former genetic divergence that is increasingly eradicated due to recent gene flow, which commenced when the snakes dispersed out of distinct glacial microrefugia in the Holocene. So to speak, the Holocene range expansions transformed the glacial ‘hotspots’ into ‘melting pots,’ in which only ‘mitochondrial ghost lineages’ bear witness of the formerly distinct evolutionary lineages. Such a scenario explains both the situation on the southern Balkan Peninsula and in the Caucasus region. There, distinct mitochondrial lineages occur in close geographic proximity but correspond to only one microsatellite cluster each (Figs 3, 6 and 7: Balkans – mtDNA lineages 3, 5, 7; Caucasus – mtDNA lineages 1, 2, 8).

For the Balkan Peninsula and southern Central Europe, this also suggests that a ‘rolling expansion and hybridization wave’ established when grass snakes with mtDNA lineage 4 (N. n. vulgaris) spread northwestward. While N. n. vulgaris successively invaded more northerly regions and ‘genetically swamped’ the resident grass snakes characterized by mtDNA lineage 3 (N. n. natrix; see Asztalos et al. 2021) there, the populations of N. n. vulgaris in their former southern refuge became increasingly admixed when other lineages, having independent glacial refuges in the southern Balkan Peninsula, expanded their ranges with Holocene warming. This lead to the current, at first glance counterintuitive, pattern that no pure populations of N. n. vulgaris exist anymore in their former glacial refuge, in contrast to their Holocene expansion area, where relatively pure populations occur in southern Central Europe and the adjacent northern Balkans.

Among the mtDNA lineages that had glacial refugia in the southern Balkan Peninsula is also mtDNA lineage 3, otherwise characteristic for the nominotypical subspecies N. n. natrix (see Fritz and Schmidtler 2020). However, the disjunct Balkan haplotypes slightly differ from the northern ones found within the range of N. n. natrix, suggesting that lineage 3 represents on the southern Balkan Peninsula an old genetic relict of a formerly wider distribution during an earlier interglacial. In contrast, the distinct Central European and Scandinavian haplotypes of lineage 3 provide evidence for a northern glacial refuge in Central Europe (Kindler et al. 2018b).

Our results corroborate that the species Natrix megalocephala is invalid (Kindler et al. 2013; Fritz and Schmidtler 2020) because it is not genetically differentiated from populations of N. natrix in the same geographic region. With respect to N. tessellata, the congruence of nuclear genomic and mitochondrial divergence in our study region suggests that the earlier described deeply divergent mitochondrial lineages (Guicking et al. 2009; Guicking and Joger 2011) reflect taxonomic differentiation. Evidence for gene flow between the two lineages in our study region supports their conspecificity, if large-scale gene flow and reproductive isolation are used as criteria for species delineation (cf. Fritz and Schmidtler 2020; Speybroeck et al. 2020). Also for grass snakes the admixture revealed in our analyses supports conspecificity.

For dice snakes additional genetic and morphological investigations are needed for a comprehensive taxonomic revision. For N. natrix the present evidence and the nomenclatural foundation laid by Fritz and Schmidtler (2020) allow assigning names to the conspecific clusters and refining the currently accepted subspecies classification. In doing so, we follow Kindler and Fritz (2018) and understand subspecies as evolutionarily significant units that differ in their nuclear genomic and mitochondrial identity and that are, in contrast to species, fully capable of extensive gene flow that can even lead to complete genetic and taxonomic amalgamation. If mitochondrial ghost lineages (i.e., higher mitochondrial diversity), mitochondrial introgression or mitochondrial capture (i.e., lower mitochondrial diversity and mitonuclear discordance) occur, nuclear genomic differentiation is decisive for subspecies delimitation.²

Within this framework (Fig. 10), the microsatellite cluster from the west of our study region (characterized by mtDNA lineage 4) is to be identified with the subspecies Natrix natrix vulgaris Laurenti, 1768 (Fritz and Schmidtler 2020). The nomenclatural situation for the south and east of our study region is less straightforward and complicated by extensive past and current hybridization, leading to the survival of many mitochondrial ghost lineages that conflict with nuclear genomic differentiation. In this region, traditionally many subspecies were recognized (Kabisch 1999; Kreiner 2007; Geniez 2015), among them also Natrix natrix cypriaca (Hecht, 1930). However, acknowledging that Cypriot grass snakes are distinct on the nuclear genomic level only due to genetic impoverishment, we propose that Cypriot populations should be taxonomically lumped together with grass snakes from the southern Balkans and western Anatolia characterized by mtDNA lineage 7. Using microsatellites, these populations cluster together with grass snakes from the southern Balkans also harboring mtDNA lineages 3 and 5. The oldest available name for these populations is Natrix natrix moreotica (Bedriaga, 1882) with a type locality on the northern Peloponnese (Greece; see Fritz and Schmidtler 2020). This renders Natrix natrix cypriaca (Hecht, 1930) from Cyprus a junior synonym of N. n. moreotica.

Grass snakes of the eastern microsatellite cluster, characterized by mtDNA lineages 1, 2, and 8, are currently assigned to two subspecies (Fritz and Schmidtler 2020). Natrix natrix scutata (Pallas, 1771), with a type locality of Atyrau (Kazakhstan), is identified with the northern populations, while the name Natrix natrix persa (Pallas, 1814), type locality Gilan and Mazandaran (Iran), is assigned to the southern populations. Our present results suggest that these two taxa should be synonymized and that the older name, N. n. scutata, should be used for all populations.

Furthermore, our results are inconclusive with respect to another subspecies that was recognized as valid by Fritz and Schmidtler (2020), namely Natrix natrix syriaca (Hecht, 1930), with a type locality of Zincirli in southeastern Turkey, near the Gulf of İskenderun. Fritz and Schmidtler (2020) tentatively identified this subspecies with grass snakes harboring mtDNA lineage 6. We could study only microsatellite data of a single individual of this lineage and this situation prevents any taxonomic conclusions, even though this sample did not seem to be distinct from other material (Figs 4–7).