It takes two to tango – Phylogeography, taxonomy and hybridization in grass snakes and dice snakes (Serpentes: Natricidae: Natrix natrix , N. tessellata )

U It takes two to tango – Phylogeography, taxonomy and hybridization in and dice Abstract Using two mitochondrial DNA fragments and 13 microsatellite loci, we examined the phylogeographic structure and taxonomy of two codistributed snake species ( Natrix natrix, N. tessellata ) in their eastern distribution area, with a focus on We found evidence for frequent interspecific hybridization, previously thought to be extremely rare, and for backcrosses. This underscores that closely related sympatric species should be studied together because otherwise the signal of hybridization will be missed. Further-more, the phylogeographic patterns of the two species show many parallels, suggestive of a shared biogeographic history. In general, the phylogeographies follow the paradigm of southern richness to northern purity, but the dice snake has some additional lineages in the south and east in regions where grass snakes do not occur. For both species, the Balkan Peninsula and the Caucasus region served as glacial refugia, with several mitochondrial lineages occurring in close proximity. Our results show that the mitochondrial divergences in both species match nuclear genomic differentiation. Yet, in the former glacial refugia of grass snakes there are fewer nuclear clusters than mitochondrial lineages, suggesting that Holocene range expansions transformed the glacial hotspots in melting pots where only the mitochondrial lineages persisted, bearing witness of former diversity. On the other hand, the deep mitochondrial divergences in N. tessellata across its entire range indicate that more than one species could be involved, even though lacking microsatellite data outside of Turkey prevent firm conclusions. On the contrary, our microsatellite and mitochondrial data corroborate that N. megalocephala is invalid and not differentiated from sympatric populations of N. natrix . For Cypriot grass snakes, our analyses yielded conflicting results. A critical assessment of the available evidence suggests that N. natrix is a genetically impoverished recent invader on Cyprus and taxonomically not distinct from a subspecies also occurring in western Anatolia and the southern Balkans. Based on combined mitochondrial and nuclear genomic evidence we propose that for grass snakes the following subspecies should be recognized in our study region: (1) Natrix natrix vulgaris Laurenti, 1768, southeastern Central Europe and northern Balkans; (2) Natrix natrix 1882), southern Balkans, western Anatolia, and Cyprus; and (3) Natrix natrix scutata (Pallas, eastern Anatolia, Caucasus region, Iran, northeastern distribution range (from eastern Poland and Finland to Kazakhstan and Lake Baikal Thus, Natrix

Natrix maura is a Western Mediterranean and Western European species, ranging from northwestern Africa and the Iberian Peninsula to central and southern France and northwestern Italy. Introduced populations also occur on the Balearic Islands and Sardinia (Schätti 1999;Guicking et al. 2008;Geniez 2015). Further east, the ecologically and morphologically similar N. tessellata replaces N. maura. Natrix tessellata occurs across most of mainland Italy and ranges from the Czech Republic and southern and eastern Austria eastward through the south of Eastern Europe to northwestern China (Xinjiang); a few isolated Central European populations (Germany, Switzerland) are witness of a formerly wider Holocene distribution in the northwest. Along the Levantine coast, N. tessellata is distributed southward to Israel and northeastern Egypt (Bannikov et al. 1977;Zhao and Adler 1993;Gruschwitz et al. 1999;Geniez 2015;Speybroeck et al. 2016). While the ranges of N. maura and N. tessellata are largely mutually exclusive and only overlap in northwestern Italy, these two water snakes occur in wide sympatry with the three grass snake species (N. astreptophora, N. helvetica, N. natrix sensu stricto).
Until a few years ago, the three grass snake species were lumped together in one polytypic species, N. natrix sensu lato (Mertens and Wermuth 1960;Kabisch 1999;Kreiner 2007). However, recent studies using mtDNA sequences and microsatellite loci revealed only limited gene flow across narrow geographic contact zones and led to their recognition as three distinct species (Pokrant et al. 2016;Kindler et al. 2017;Schultze et al. 2019Schultze et al. , 2020Asztalos et al. 2020Asztalos et al. , 2021aSpeybroeck et al. 2020). Natrix astreptophora occurs in the northern Maghreb region, the Iberian Peninsula and southwestern France (Mertens and Wermuth 1960;Kabisch 1999). In southwestern France lies a narrow contact zone of this species and N. helvetica, where limited bidirectional gene flow occurs (Pokrant et al. 2016;Asztalos et al. 2020). This contact zone is a bimodal hybrid zone (see Jiggins and Mallet 2000), with pure parental genotypes and hybrid genotypes occurring together or in close proximity. Natrix helvetica is distributed in France, Great Britain, the Benelux countries, western Germany, Switzerland, Italy and across the Alps to southern Bavaria and Tyrol (Geniez 2015;Schultze et al. 2020;Asztalos et al. 2021a). The third grass snake species is Natrix natrix sensu stricto, which occurs from the Rhine region in the west to Lake Baikal in the east. Its vast distribution range also includes Fennoscandia, northeastern Italy, the Balkan Peninsula and Anatolia, the Transcaucasus, and northern Iran (Geniez 2015; Kindler et al. 2017;Fritz and Schmidtler 2020;Schultze et al. 2020). Three narrow bimodal hybrid zones between N. helvetica and N. natrix have been described. One is located in the Rhine region, where gene flow is mainly unidirectional from N. h. helvetica into N. natrix Schultze et al. 2019). In the two other hybrid zones another subspecies of the barred grass snake is involved, N. h. sicula (Cuvier, 1829). It hybridizes with N. natrix in northeastern Italy ) and, beyond the Alps, in southern Bavaria and Tyrol (Asztalos et al. 2021a).
Each of the three grass snake species is polytypic. Recently, Fritz and Schmidtler (2020) reconciled traditional morphology-based subspecies with genetically differentiated units. These authors presented a classification scheme recognizing geographically vicariant, evolutionarily significant genetic clusters as distinct subspecies. Only for the Balkan Peninsula, Anatolia and adjacent regions including the Caucasus and Transcaucasus, the intraspecific variation of N. natrix could not be assessed because no nuclear-genomic evidence was then available. While intraspecific variation in grass snakes is now, with the mentioned notable exception, well understood, for N. maura and N. tessellata the situation is strikingly different. Both species contain genetically deeply divergent parapatric lineages (Guicking et al. 2008(Guicking et al. , 2009Guicking and Joger 2011) suggestive of taxonomic divergence. Yet, in contrast to grass snakes, traditionally no subspecies have been recognized for N. maura (Mertens and Wermuth 1960;Schätti 1999;Geniez 2015;Speybroeck et al. 2016). For N. tessellata, only the micro-endemic subspecies Natrix tessellata heinrothi (Hecht, 1930) from the small island of Serpilor in the Black Sea has been recognized by some authors. However, its validity has been disputed (see Gruschwitz et al. 1999;Geniez 2015).
In the present study we aim to resolve the taxonomic and genetic differentiation of grass snakes (N. natrix sensu stricto) in the Balkan Peninsula, Anatolia, the Caucasus region and Iran. In addition, we scrutinize whether the phylogeographic and taxonomic differentiation of the dice snake (N. tessellata) parallels that of N. natrix in our study region, and we examine, for the first time, whether and, if so, to which extent, natural hybridization occurs between these two sympatric species. For doing so, we apply an established toolbox consisting of combined evidence from mitochondrial DNA (mtDNA) sequences and 13 nuclear genomic microsatellite loci (Kindler et al. 2013(Kindler et al. , 2018aKindler and Fritz 2018;Pokrant et al. 2016;Asztalos et al. 2020Asztalos et al. , 2021aSchultze et al. 2019Schultze et al. , 2020 to examine 683 grass snakes and 332 dice snakes from our immediate study region and beyond. We also reassess the status of N. megalocephala Orlov and Tuniyev, 1987 based on nuclear genomic microsatellite loci. This western Caucasian species was described based on slight morphological differences to putatively syntopic N. natrix (see Orlov andTuniyev 1987, 1992). However, in the face of inconsistent morphological characters and the lack of any distinctiveness in mtDNA, the validity of N. megalocephala was soon challenged and the name was synonymized with N. natrix. Natrix megalocephala is currently identified as a junior synonym of the subspecies N. n. scutata (Pallas, 1771) (Kindler et al. 2013; see the review in Fritz and Schmidtler 2020).

Sampling and laboratory procedures
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. The remaining data originate from earlier studies (Guicking et al. 2006(Guicking et al. , 2009Kindler et al. 2013Kindler et al. , 2017Kyriazi 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. megaloce phala 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(Kindler et al. , 2018aPokrant et al. 2016;Kindler and Fritz 2018;Schultze et al. 2019Asztalos et al. 2020Asztalos et al. , 2021aAhnelt 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), andPokrant 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).

Mitochondrial DNA sequences and network analyses
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(Guicking et al. , 2009Kindler et al. 2013Kindler et al. , 2014Kindler et al. , 2017Kindler et al. , 2018aKyriazi et al. 2013;Kindler and Fritz 2018;Schultze et al. 2019Schultze et al. , 2020Asztalos 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(Guicking et al. , 2009Kindler et al. 2017Kindler et al. , 2018aSchultze et al. 2019Schultze et al. , 2020Ahnelt 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.

Microsatellite cluster analyses, infer ring hybrid status and PCAs
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 1 , F 2 , 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.

Inference of populationgenetic para meters
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).

Haplotype networks for grass snakes (Natrix natrix)
The network for cyt b of Natrix natrix contained several loops (Fig. 1). The haplotype cluster corresponding to lineage 7 was located in the center of the network and connected to the haplotype clusters of lineages 1, 2, 3, 5, 6, and 8. Lineage 7 comprised 31 haplotypes that differed by a maximum of 16 mutations. The haplotype cluster of lineage 7 was connected by a minimum of 27 mutation steps to the haplotypes of lineage 2 and by a minimum of 29 mutations to those of lineage 1. The haplotypes of lineages 1 and 2 differed from one another by a minimum of 12 mutations. Among the five haplotypes of lineage 2 a maximum of 15 mutations occurred; among the six haplotypes of lineage 1, a maximum of eight mutations. Haplotypes of lineage 5 were connected to haplotypes of lineages 1, 2, and 7 by a minimum of 58, 56, and 39 mutations, respectively. Lineage 5 comprised 25 haplotypes that differed by a maximum of seven mutations. Haplotypes of lineage 4 were only connected to lineage 5 by a minimum of 11 mutations. The 71 haplotypes of lineage 4 differed by a maximum of 18 mutations. Haplotypes of lineage 6 were connected to haplotypes of lineages 5 and 7 and differed by a minimum of 50 mutations from haplotypes of lineage 5 and by a minimum of 37 mutations from those of lineage 7. Lineage 6 comprised two haplotypes that differed by one mutation. Haplotypes of lineage 7 were also connected with haplotypes of lineages 3 and 8. Haplotypes of lineage 3 differed by a minimum of 44 mutations from haplotypes of lineage 7 and by a minimum of 57 mutations from those of lineage 8. Among the 49 haplotypes of lineage 3 occurred up to seven mutations. Haplotypes of lineage 8 were connected by a minimum of 21 mutations to haplotypes of lineage 7. Lineage 8 consisted of 36 haplotypes that differed by a maximum of 19 mutations.
The network for ND4+tRNAs sequences of N. natrix (Fig. S1) resembled the cyt b network, but with fewer mutations, and the clusters of lineages 4 and 5 were not clearly distinct. Again, the haplotype cluster of lineage 7 was located in the center of the network, with direct connections to haplotypes of lineages 1, 3, 5, 6, and 8. The 30 individual haplotypes of lineage 7 differed by up to 14 mutations. The haplotypes of lineage 7 were connected to haplotypes of lineage 1 by a minimum of 35 mutations. Lineage 1 comprised three haplotypes that differed by a maximum of four mutations. Haplotypes of lineage 2 were connected to haplotypes of lineage 1 by a minimum of four mutations. The four haplotypes of lineage 2 differed by a maximum of three mutations. In addition, haplotypes of lineage 7 were connected over at least 25 mutations to haplotypes of lineage 5. The latter lineage comprised 20 individual haplotypes that differed by a maximum of eight mutations. Lineage 5 was connected with haplotypes of lineage 4 by a minimum of only 2 mutations. Lineage 4 comprised 39 haplotypes that differed by a maximum of ten mutations. The only haplotype of lineage 6 was connected to lineage 7 by a minimum of ten mutations. Furthermore, haplotypes of lineage 7 were also connected to those of lineages 3 (minimum 29 mutations) and 8 (minimum 16 mutations). Haplotypes of lineages 3 and 8 were separated by a minimum of 37 mutations. The 48 haplotypes of lineage 3 and the 21 haplotypes of lineage 8 differed both by a maximum of nine mutations.
The following haplotypes for ND4+tRNAs were newly identified in the present study: gn17-gn21, gy13-gy31, and r39 (European Nucleotide Archive accession numbers OU862605-OU862629). All currently known haplotypes for N. natrix and their accession numbers are listed in the Supplementary Information (Table S6).

Haplotype networks for dice snakes (Natrix tessellata)
For cyt b of Natrix tessellata (Fig. 2), the 32 haplotypes of lineage T were located in the interior of the network of all haplotype clusters. Direct connections existed to the haplotypes of lineages A, E, G, I, and K. Among the haplotypes of lineage T occurred a maximum of 54 mutations. The haplotypes of lineage T were connected by a minimum of 46 mutations to those of lineage G and by a minimum of 43 mutations to the haplotypes of lineage I. The haplotypes of lineages G and I were separated by a minimum of 87 mutations. The 15 haplotypes of lineage G differed by a maximum of 47 mutations, while the three haplotypes of lineage I differed by a maximum of six mutations. Lineage E comprised 66 haplotypes, which differed by a maximum of 19 mutations. The haplotypes of lineage E were connected by a minimum of 35 mutations to lineage T. The haplotypes of lineage C were connected to those of lineage E by a minimum of 26 mutations. The four haplotypes of lineage C differed by a maximum of five mutations. Lineage T was also connected to the haplotypes of lineages A and K and differed by a minimum of 19 mutations from lineage A and by a minimum of 21 mutations from lineage K. The haplotypes of lineages A and K differed by a minimum of 24 mutations from one another. Lineage K comprised 12 haplotypes that differed by a maximum of 12 mutations. Lineage A was comprised of 12 haplotypes that differed by a maximum of 11 mutations. Lineage A differed from lineages J and U by a minimum of 25 and 19 mutations, respectively. The three haplotypes of lineage J differed by a maximum of three mutations, the six haplotypes of lineage U by a maximum of 12 mutations.
The following cyt b haplotypes were newly identified in the present study: E52-E61 and T16-T24 (European Nucleotide Archive accession numbers OU862672-OU862690).
The second network (ND4+tRNAs) for N. tessellata was restricted to two lineages, E and T (Fig. S2). The haplotypes of the two lineages were connected over two pathways by a minimum of 44 mutations each. Lineage E corresponded to eight haplotypes, which differed by a maximum of six mutations. The 10 haplotypes of lineage T differed by up to 28 mutations.
All haplotypes for ND4+tRNAs were newly identified in the present study: E1-E8 and T1-T10 (European Nucleotide Archive accession numbers OU862630-OU862647). These haplotypes and all currently known cyt b haplotypes for N. tessellata and their accession numbers are listed in the Supplementary Information (Table S6).

Mean uncorrected p distances
Using distinct haplotypes, the average uncorrected p distance between dice snakes and grass snakes was for the cyt b gene 11.50%, and for the DNA coding for the partial ND4 gene plus tRNAs 12.10% (Table S7). For the eight mitochondrial lineages of Natrix natrix, divergences of 0.96-4.21% were observed for the cyt b gene, and for ND4+tRNAs, 0.60-5.17%. The nine mitochondrial lineages of N. tessellata were more divergent. In the cyt b gene, they differed on average by 1.34-6.15%, and for ND4+tRNAs haplotypes, a divergence of 6.93% was observed for lineages E and T (Table S8). The highest divergence values within N. tessellata resemble those between the three grass snake species N. astreptophora, N. helvetica, and N. natrix (Kindler et al. , 2018a.

Geographic distribution of mitochondrial lineages
When the geographic distribution of the individual mitochondrial lineages of the two Natrix species is compared (Fig. 3), some similarities and some differences become apparent. In both species, there is a widely distributed Balkan lineage (lineage 4 in N. natrix; lineage E in N. tessellata). Beyond the Balkan Peninsula and beyond the scope of the present study, these lineages also encroach on regions further north and west (Guicking et al. 2009; Guicking and Joger 2011;Kindler et al. 2013Kindler et al. , 2017. On the Peloponnese, and in N. natrix also northwest of the Peloponnese, this Balkan lineage is replaced by another lineage (lineage 5 in N. natrix; lineage G in N. tessellata). However, in N. natrix, there are two other lineages (3 and 7) in the southeastern Balkan Peninsula, western Anatolia and Cyprus, whereas in N. tessellata the widely distributed Balkan lineage E also occupies the southeastern Balkan Peninsula and western Anatolia. In both species there is a divide in Anatolia between a western lineage (lineage 7 in N. natrix; lineage E in N. tessellata) and an . Sampling sites and geographic distribution of mtDNA lineages for 653 grass snakes (Natrix natrix; top) and 325 dice snakes (N. tessellata; bottom) used in this study (eastern parts of the distribution ranges). For N. tessellata, data outside Turkey are from Guicking et al. (2006Guicking et al. ( , 2009 and Kyriazi et al. (2013). Insets: Natrix natrix (top; Ropotamo, Bulgaria) and N. tessellata (bottom; Ömeroba, Edirne, Turkey). Photos: Daniel Jablonski. eastern lineage (lineage 8 in N. natrix; lineage T in N. tessellata). Yet, in N. tessellata the eastern lineage is found on Cyprus (albeit based on only a single specimen), and in N. natrix the western lineage occurs on the island. Furthermore, in N. natrix there is a unique endemic lineage in the Gulf of İskenderun (lineage 6), whereas the range of N. tessellata reaches further south to northeastern Egypt, where another distinct lineage (J) occurs. In both species, several mitochondrial lineages meet in the Caucasus region 1 , with a unique endemic lineage (2) in N. natrix in the eastern Caucasus region. In both species, a widely distributed lineage occurs north of the Black Sea and the Caspian Sea (lineage 8 in N. natrix; lineage A in N. tessellata), however, lineage 8 in N. natrix encroaches on eastern Anatolia and occupies there the area where lineage T in N. tessellata occurs. East of the Caspian Sea, two additional lineages occur in N. tessellata (lineages K and U), in regions where N. natrix is absent.

Interspecific hybridization and the status of Cypriot snakes
For each of the two species, the 13 studied microsatellite loci were highly polymorphic, with 9-32 alleles per locus for Natrix natrix and 3-19 alleles per locus for N. tessellata. Total allele numbers for N. natrix were 229 and for N. tessellata, 158 (for allele size range and individual allele numbers per locus, see Table S3).
After initial analysis of our total dataset (n = 705) using structure 2.3.4, hidden structuring was examined using hierarchically nested subsets corresponding to each of the two species and genetic clusters within either species, with admixed individuals in each step removed. The respective results were then compared to the mitochondrial identity of each sample. In addition to the structure analyses, PCAs were calculated for the same datasets to scrutinize the structure results using an approach free of population-genetic assumptions.
For the structure run using all 705 genotypes, the ΔK method identified K = 2 as the optimal number of clusters. However, a second pronounced peak occurred for K = 5 (Fig. S3A). Under K = 2, one cluster matched Natrix natrix and the other N. tessellata (Fig. 4); 555 genotypes were assigned to N. natrix, 77 to N. tessellata, and the remaining 73 were identified as interspecific hybrids. The genotypes of 'N. megalocephala' were either not dif-1 However, an isolated record of lineage 7 in the western Caucasus region of Turkey (Borçka, Artvin; ZDEU DNA1284) most likely refers to a translocated grass snake or a confused sample. The microsatellite identity of ZDEU DNA1284 matches its mtDNA haplotype, whereas four other snakes from Borçka (ZDEU DNA1283, DNA1287, DNA1288, ZFMK 71145; Table S1) correspond to the expectation for this region and represent the local microsatellite cluster. For one of these snakes the mitochondrial haplotype is known and it belongs to the expected lineage 8. ferentiated from N. natrix (n = 2) or identified as hybrids (n = 2).
Most of the interspecific hybrids were from Turkey, but hybrid signatures were also found in the Balkans, Transcaucasia and Russia. Sixty-six hybrids were morphologically identified as grass snakes (64 as N. natrix, 2 as 'N. megalocephala') and the remaining seven hybrids were morphologically identified as N. tessellata. The majority of hybrids (49) held mitochondrial haplotypes of N. natrix, three had haplotypes of N. tessellata, and for the remaining snakes the mitochondrial identity is unknown (Table S1). Under K = 2, structure assigned most samples from Cyprus to the cluster for pure N. tessellata, even though the snakes were morphologically unambiguous N. natrix. A few other Cypriot samples were genotypically assigned to N. natrix or identified as interspecific hybrids. Yet, under K = 5, the genotypes of the controversial Cypriot snakes were not lumped together with N. tessellata but constituted an extra cluster ( Fig. 4 bottom: yellow cluster). Otherwise, N. tessellata remained as a distinct cluster, but the former cluster of N. natrix was further subdivided, so that there was now one cluster for N. tessellata, another one for Cypriot snakes, and three clusters for the remaining N. natrix. Under K = 5, there was also evidence for interspecific hybridization ( Fig. 4; Table S1), even though fewer putative hybrids occurred and all four genotypes of 'N. megalocephala' were now placed into the cluster of pure N. natrix. In a PCA using the same dataset as for the structure analyses, N. tessellata and N. natrix were clearly distinct, although intermediate genotypes occurred. The genotypes of all 'N. megalocephala' clustered within N. natrix. However, most Cypriot snakes were placed completely outside of N. tessellata and N. natrix; a few Cypriot samples clustered within N. natrix (Fig. 5). Population-genetic parameters for all Cypriot samples (n = 27) revealed very low genetic variability with evidence for inbreeding compared to other populations harboring the same mitochondrial lineage (i.e., lineage 7; Table S9).
In the face of these conflicting results and the morphological identity of the Cypriot snakes, we used the K = 5 scenario as the basis for further calculations. For calculating Q values in hybridlab to determine interspecific hybrids, the individual Q values of the three clusters for N. natrix and of the Cypriot cluster were merged and opposed to the Q values of the cluster for N. tessellata. This resulted in thresholds of 93% for pure N. natrix and 94% for pure N. tessellata.

Genotypic structuring within Natrix natrix
Based on these thresholds for pure Natrix natrix and N. tessellata, a second structure run was conducted (Fig.  6A) for which all samples of N. tessellata or with genetic impact of N. tessellata were excluded. For this dataset (n = 626), the ΔK method inferred K = 2 as the optimal number of clusters (Fig. S3B). One cluster (green in Fig.  6A) corresponded to populations in the west of our study region. Its geographic distribution matched well that of mtDNA lineage 4 of N. natrix (red in the bar for the mi-  Figure 1. However, the two mitochondrial lineages of N. tessellata are shown in black to contrast them from N. natrix. White indicates missing mtDNA data. Samples are arranged from west to east. In the structure diagrams, an individual is represented by a vertical segment that reflects its ancestry. The black arrows highlight samples from Cyprus; the dark red arrows, samples identified as N. megalocephala; and the pink arrow indicates the only grass snake with mitochondrial lineage 6 of N. natrix. tochondrial identity in Fig. 6A). The second cluster (pink in Fig. 6A) represented populations in the south and east of our study region and corresponded to several mtDNA lineages distributed there (mtDNA lineages 1, 2, 3, 5, 6, 7, 8; remaining colors in the bar for the mitochondrial identity in Fig. 6A). This cluster also included the four samples from grass snakes morphologically identified as N. megalocephala for which microsatellite data are available (ZDEU DNA1287, DNA1288, ZFMK 58029, 60732; Table S1). Admixture between the two clusters occurred mainly in the southern Balkan Peninsula (173 of 205 hybrids) and matched the geographic contact zone of the two clusters. The corresponding PCA showed only weak differentiation, with the western cluster matching mtDNA lineage 4 only slightly differentiated and with massive admixture (Fig. 7A). Yet, many Cypriot samples were again highly distinct.
In the third structure analysis, the cluster of N. natrix from the south and east, matching several distinct mtDNA lineages, was examined. Samples representing pure N. natrix of the western cluster and admixed snakes were removed. For the 159 processed samples, the ΔK method revealed K = 3 as the optimal solution (Fig.   S3C). A mainly western cluster (purple in Fig. 6B) largely corresponded to the mitochondrial lineages 3, 5 and 7, whereas an eastern cluster (beige in Fig. 6B) largely matched mitochondrial lineages 1, 2 and 8. A third cluster corresponded to most of the samples from Cyprus (21 of 26 snakes). Ten individuals were identified as admixed.
Notably, not all admixed snakes were from the geographic contact zone of the two clusters (except Cyprus cluster). In addition, one sample from Transcaucasia (Azerbaijan) was assigned to the purple cluster (MTD T 3680) and another one from there was admixed (MTD D 48550). The PCA generally also confirmed for this dataset weak differentiation and admixture, but again 21 of the 26 Cypriot samples were clearly distinct (Fig. 7B).
The results of the second and third structure analyses largely mirrored the outcome of the first structure calculation for K = 5 (Fig. 4) in that the four clusters for grass snakes from the first calculation match those also identified in the subset data. The same dataset was used as for structure analyses. Samples are colored according to their mitochondrial lineage or structure cluster membership; admixed individuals in structure analyses are shown in grey. Thresholds for admixed samples were derived from hybridlab simulations and the same as used for structure. Oval outlines correspond to 95% confidence intervals. PC1 explains 3.96% of variance, PC2 2.13%, and PC3 1.70%. Some snakes morphologically matching N. tessellata are highlighted. The pink arrows denote the only grass snake corresponding to mitochondrial lineage 6 of N. natrix. Figure 6. Genotypic structuring of Natrix natrix in the study region. Symbol colors indicate structure cluster membership; admixed ancestries are depicted as pie chart sectors according to Q values. In the bar plots below the maps, samples are represented by two bars, indicating their mitochondrial lineage (top) and their inferred cluster membership (bottom). Mitochondrial lineages are color-coded as in Figure 1; white bars indicate missing mtDNA data. Samples are arranged from west to east. The black arrows highlight samples from Cyprus; dark red arrows, samples identified as N. megalocephala; and the pink arrow indicates the only grass snake with mitochondrial lineage 6 of N. natrix. A most likely translocated grass snake from Borçka (Artvin, Turkey) representing the western cluster is not shown in the map (see footnote 1).

Figure 7.
Principal Component Analyses for microsatellite data of (A) pure Natrix natrix (n = 626) and (B) N. natrix of the southern and eastern cluster (n = 159). The same datasets were used as for structure analyses. Samples are colored according to their mitochondrial lineage or structure cluster membership. Thresholds for admixed samples were the same as used for structure. Oval outlines correspond to 95% confidence intervals. For (A) PC1 explains 2.82% of variance, PC2 2.18%, and PC3 1.66%; for (B) PC1 explains 4.78% of variance, PC2 3.52%, and PC3 2.57%. The pink arrows highlight a snake corresponding to mitochondrial lineage 6 of N. natrix.

Genotypic structuring within Natrix tessellata
In a fourth structure run, only our 52 samples of pure Natrix tessellata from Turkey were processed (Fig. 8). For this dataset, the ∆K method inferred K = 4 as the optimal number of clusters (Fig. S3D). In the following, we refer to the cluster colors in Figure 8. A purple cluster corresponded to samples from Çanakkale, a green cluster to samples from the Lake Işıklı region (Denizli), a yellow cluster to samples from the Karamık Swamp (Afyonkarahisar), and an orange cluster to most samples from Lake Çıldır (Ardahan) and other sites in eastern Anatolia. The distribution of the orange cluster seems to be disjunct and interrupted by the occurrence of the yellow and green clusters. Admixture seems to be widespread. The purple and the green clusters largely match mtDNA lineage E (bluish green in Fig. 8) and the yellow and orange clusters largely match mtDNA lineage T (blue in Fig. 8), so that two structure clusters each correspond to one mito-chondrial lineage. This differentiation pattern is also confirmed by the PCAs (Fig. 9), which supports both weak differentiation and admixture.

Discussion
The most prominent results of our study are (1) the evidence for hybridization between grass snakes (Natrix natrix) and dice snakes (N. tessellata), (2) the pronounced divergence of Cypriot grass snakes, (3) the similarity of phylogeographic differentiation patterns of N. natrix and N. tessellata, and (4) that the eight mitochondrial lineages of N. natrix correspond to less microsatellite clusters. These results have implications both for phylogeography in general and for the taxonomy of the two Natrix species. In addition, our results confirm that N. megalocephala is not a valid species.

Interspecific hybridization
Previous investigations focused either on the phylogeography of N. natrix (e.g., Kindler et al. 2013Kindler et al. , 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.

Status of Cypriot grass snakes
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 Figure 9. Principal Component Analyses for microsatellite data of pure Natrix tessellata (n = 52). The same dataset was used as for structure analyses. Samples are colored according to their mitochondrial lineage or structure cluster membership. Thresholds for admixed samples were the same as used for structure. Oval outlines correspond to 95% confidence intervals. PC1 explains 6.65% of variance, PC2 5.23%, and PC3 5.08%. 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 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.

Comparative phylogeographies of Natrix natrix and N. tessellata
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 wellknown 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, andN. 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 . 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).

Taxonomic implications
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. 2 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 val-id 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).

Conclusions
Our study showcases that phylogeographic investigations should include samples of codistributed congeners because otherwise evidence for interspecific hybridization will be missed. Genomic approaches are expected to be more informative than investigations based on microsatellite loci alone. In our study system, we identified hybridization between Natrix natrix and N. tessellata, two species that were previously thought to hybridize only rarely (Mebert et al. 2011). Our study also exemplifies that the study of codistributed taxa allows the direct comparison of phylogeographic differentiation patterns and more general inferences about biogeography. Furthermore, our findings support that the software structure is prone to errors regarding the number of clusters when no additional evidence, independent from population-genetic assumptions, is considered. With respect to N. tessellata our results suggest that the previously identified mitochondrial lineages reflect taxonomic divergence, as in N. natrix. Further research is needed in N. tessellata to examine nuclear genomic variation beyond Turkey. The deep mitochondrial lineages in this species (Guicking et al. 2009;Guicking and Joger 2011) and the nuclear genomic fingerprints presented in Guicking et al. (2009) suggest that more than only a single species could be involved, although this is not confirmed for the two mitochondrial lineages studied herein. For N. natrix, we propose a refined intraspecific classification and recognize three subspecies from our study region, bringing the total number of subspecies to four (Fig. 10): Natrix natrix vulgaris Laurenti, 1768 from southeastern Central Europe and the northern Balkans, Natrix natrix moreotica (Bedriaga, 1882) from the southern Balkans, western Anatolia and Cyprus, and Natrix natrix scutata (Pallas, 1771) from eastern Anatolia, the Caucasus region, Iran, and a wide area from eastern Poland and Finland to Kazakhstan and the Lake Baikal region in Russia and adjacent Mongolia. All of these three subspecies comprise individuals with and without back stripes, with an increasing tendency to the presence of back stripes and other coloration variations in the south of the range. The status of Natrix natrix syriaca (Hecht, 1930) from the Gulf of İskenderun, Turkey, cannot be resolved due to insufficient material. Beyond our study region, Natrix natrix natrix (Linnaeus, 1758) inhabits the northwestern distribution range of the species (Central Europe and Scandinavia).  Sindaco et al. (2013), GBIF (www. gbif.org), and iNaturalist (https://www.inaturalist.org), combined with our genetically verified records, following the approach described in Asztalos et al. (2020). Insets: Natrix natrix moreotica (bottom left; Peloponnese, Greece) and N. n. scutata (top right; Białowieża, easternmost Poland). Photos: Henrik Bringsøe.