We obtained morphological data from 541 individuals of N. tessellata across its entire species range (Figs 1B, S1 A, B; Tables S1, S2) with a focus directed towards populations from the ‘Europe’ and ‘Greece’ clades. The data were collected from individuals stored in museum collections [ Zoological Research Museum Alexander Koenig in Bonn, Germany (ZFMK), Natural History Museum in Vienna, Austria NHMW(), the Hungarian Natural History Museum in Budapest, Hungary (HNHM)], as well as from our field effort. The morphological and osteological data analysed in this study are available in the Supplementary files or can be obtained from the corresponding author of the study (DJ) at the herpetological collection of Comenius University in Bratislava (CUHC). The obtained morphological dataset was divided according to the ranges of currently recognised evolutionary clades of N. tessellata (Table S1), as defined by Guicking et al. (2009) and Jablonski et al. (2024) (Figs 1A, B and S1A, B).

We obtained morphological data from four individuals of the ‘Iran’ clade, 179 from the ‘Greece’ clade, 27 from the ‘Jordan’ clade, six from the ‘Crete’ clade, 269 from the ‘Europe’ clade, 20 from the ‘Anatolia’ clade (the ‘Turkey’ lineage sensu Guicking et al. 2009) and 36 from the ‘Central Asia’ clade (sensu Jablonski et al. 2024), which includes 16 individuals from the ‘Caspian’ alias the ‘Caucasus’ lineage, six from the ‘Kazakhstan’ lineage and 14 from the ‘Uzbekistan’ lineage as previously defined by Guicking et al. (2009). Basic morphometric and meristic data were collected according to Laňka (1978) and Mebert (1993). We used a calliper to measure head dimensions to the nearest ± 1 mm and a tape to measure body dimensions. The following 13 morphometric or meristic data were taken: head length (HL; from the tip of the rostrum to the posteriormost edge of pariental scales) , head width (HW; the widest part of the head measured behind the eyes) , mouth length (ML; from the tip of the dentary bone to the posteriormost edge of the mouth opening) , snout-vent length (SVL; from the tip of the rostrum to the edge at the cloacal opening) , tail length (TL; from the edge at the cloacal opening to the end of the tail tip) , total length (TotL; the sum of SVL and TL) , number of preocular scales (PREOC; scales bordering the eyes anteriorly and positioned between the supralabials and the supraocular-prefrontal scales , postocular scales (POSTOC; scales bordering the eyes posteriorly and positioned between the supralabials and the suparocular-parietal scales , supralabial scales (SUPL; lateral labial scales of the upper jaw below the eyes) , sublabial scales (SUBL; lateral labial scales of the lower jaw) , ventral scales (VENT; ventral scales beginning where they touch the lowermost dorsal scale row on both sides up to the anal plate, excluding the preventrals) , subcaudal scales (SUBC; subcaudal scales from the first pair posterior the anal plate to the end of the tail, excluding the terminal tail tip) and dorsal scale rows (DORS; the number of dorsal scale rows of the trunk counted mid-body). Bilateral characteristics such as PREOC, POSTOC, SUPL and SUBL were counted and evaluated separately on the right and left sides of the head. To ensure unity and accuracy of the morphometric data, almost all specimens were examined by the first author of this study. We determined the sex of the examined specimens, based on the shape (width) of the tail base across the first 10 subcaudal scales (sensu Mebert 1993). Juveniles (SVL ˂ 240 mm) were not sexed due to the possibility of incorrect sex identification. For the purpose of the study, the following numeric range was defined (Günther 1996; Moravec 2015 and literature therein): juveniles SVL ≤ 240 mm, subadult males SVL ≥ 240.1 ≤ 480 mm, adult males SVL > 480 mm, subadult females SVL ≥ 240.1 ≤ 550 mm, adult females SVL > 550 mm. Individuals with SVL > 480 mm were considered adult males capable of courting females and those with SVL > 550 mm were considered adult females capable of developing eggs (Luiselli and Rugiero 2005). Morphologically investigated individuals were assigned to mitochondrial clades according to previous molecular studies (Guicking et al. 2009; Asztalos et al. 2021; Jablonski et al. 2024). Populations from within the expected contact zone in the south-western Balkans, particularly in Albania, were categorised as the ‘Greece’ clade, based on the genetic data from an upcoming study (Jablonski et al., unpubl. data). The decision resulted from initial analyses, which suggested a morphological similarity between the studied specimens and the ‘Greece’ clade in central-south Albania. It is important to note that, from three evolutionary clades present in Greece, one is representing exclusively an island population (Crete; Figs 1B and S1B).

From each studied individual, we also recorded overall body and head colouration (Figs 2, 3) and apparent anomalies in the scalation. For the colour pattern variations within European populations of N. tessellata, our study used mostly online citizen-science data, particularly iNaturalist.org (www.inaturalist.org), Observation.org (www.observation.org) and Balcanica.info (www.balcanica.info) to analyse the amount of pigmentation of labials in 767 N. tessellata individuals (720 of ‘Europe’ clade and 47 of the ‘Greece’ clades; Fig. 4A, Table S3). The goal was to determine statistical differences in this characteristic between the two clades. For this purpose, we selected photographs with clearly visible labials and precise location information. Then we categorised labial colouration into four types: (i) type 1, where both sublabials and supralabials are coloured similarly to the rest of the head; (ii) type 2, where white pigment is present only on sublabials and/or less than 50% of the surface of each supralabial scale, excluding the last one; (iii) type 3, where white pigment is present on sublabials and approximately 50% or more of the surface of each supralabial scale, excluding the last one; and (iv) type 4, where white pigment covers the entire surface of both sublabials and supralabials, except for the last supralabial scale, with a thin strip of black pigment on the posterior edge of the supralabials (see Fig. 4B).

Descriptive statistics were applied to all 541 examined specimens (Table S2). Individuals with missing values were excluded from the following statistical analyses, as well as hatchlings and juveniles, to remove potential effects of allometry and ontogeny. Ultimately, a total of 382 individuals from all seven clades (inclusive of three lineages from Central Asia separately, as mentioned earlier) were analysed statistically.

Initially, we investigated both metric and meristic characters amongst all clades/lineages using Principal Component Analysis (PCA). Allometric body size corrections implemented in the function allom from the package GroupStruct (Chan and Grismer 2021) were applied to the metric data to eliminate size-related information. To maximise intergroup variation, Discriminant Analysis of Principal Components (DAPC; Jombart and Collins 2015) was performed with a priori selected clusters as genetically defined clades. For the optimal number of PCs retained for DAPC, we followed Thia (2022), where the optimal number was determined as k-1, with ‘k’ representing the number of clusters. This criterion captures maximal amongst-group variation. We adhered to the recommended standard reporting for DAPC, as suggested by Miller et al. (2020).

Subsequently, we used the framework in accordance with Chan and Grismer (2021) for the comparison of the ‘Europe’ and ‘Greece’ clades. Sexual dimorphism was investigated by Multivariate Analysis of Variance (MANOVA). Based on these results, the data were analysed separately for males and females. Boxplots were constructed to visualise the range, median, significance, the data distribution and extent of differences between metric and meristic characteristics, while also identifying any outliers. The normality of the data was tested by the Shapiro-Wilk test, followed by the Mann-Whitney U test used for identification of statistically significant mean differences between both clades in each metric and meristic character. To visualise the morphospatial relationship between the studied clades, Principal Component Analyses (PCA) were conducted separately for males and females (see Figs 1C and S2A), while two subadult males from Alpnach Lake (Switzerland) were omitted from this analysis due to a high rate of abnormalities in the scalation (see Mebert 2011c). DAPC was performed in addition to PCA. To investigate linear relationships between head width (HW) and snout-vent length (SVL), we created sex-specific scatter plots with linear regression lines associated with both clades. Furthermore, a chi-square test was applied to investigate differences in the labial colouration between the ‘Europe’ and ‘Greece’ clades. The nature of the association and contributions to the chi-square score were subsequently visualised by plotting standardised Pearson’s residuals between expected and observed values. All statistical analyses were carried out using R (v.4.3.2; R Core Team 2023) with the following packages: poppr (Kamvar et al. 2014), corrplot (Wei and Simko 2016), ggplot2 (Arnold et al. 2018), tidyverse (Wickham et al. 2019), factoextra (Kassambara and Mundt 2020), vegan (Oksanen et al. 2020), psych (Revelle 2020), GroupStruct (Chan and Grismer 2021), dplyr (Wickham et al. 2021), ggsignif (Ahlmann-Eltze 2022) and MorphoTools2 (Šlenker et al. 2022) packages.

The osteological examination focused solely on the braincase area. The skulls used for geometric morphometrics and micro-CT analyses were collected from dead specimens found in the wild or obtained from the Zoological Research Museum Alexander Koenig, Bonn, Germany (ZFMK) and the Natural History Museum, Vienna, Austria (NHMW; see Table S5).

The skulls that were not analysed through micro-CT (18 specimens) were first cleared by dermestid beetles (genus Dermestes), then placed in cold water at a temperature of 20°C to drain the blood from the bones, followed by immersion in hot water at an initial temperature of 60°C to remove the grease. Subsequently, the skulls were soaked in a 10% hydrogen peroxide (H₂O₂) solution for final cleaning and bleaching with each step lasting 24 hours. The braincases of the disarticulated skulls were then scanned at various depths of field using a digital camera Canon EOS 5D MARK IV attached to a Zeiss AXIO Zoom.V16 microscope. Individual photographs were later composited using the programme Zerene Stacker v. 1.04 (Zerene Systems, Richland, WA, USA).

In our micro-computed tomography methods, we employed the high-resolution micro-CT scanner v|tome|x L 240, located at the Laboratory of Computed Tomography (100 kV maximum voltage; 250 µA maximum current; 25 µm maximum voxel size; 333 ms timing) within the Earth Science Institute of the Slovak Academy of Sciences in Banská Bystrica, Slovakia. Obtained data were reconstructed in Phoenix datos|x 2.0 CT Scanning Software, saved as .vgl files and subsequently micro-CT two-dimensional (2D) images were analysed using Avizo 8.1 software. This step involved generating individual three-dimensional (3D) colour-coded images of the braincase with all its elements of nine evolutionary groups. Moreover, rigid braincase units, such as three bones – frontal, parietal and supraoccipital on the dorsal side and two bones – basioccipital and parabasisphenoid on the ventral side of each braincase, were separated for detailed comparisons between the ‘Europe’ and ‘Greece’ clades. The description of the braincase structures and elements follows the terminology and abbreviations of von Szunyoghy (1932), Szyndlar (1984), Racca et al. (2020), and Seghetti et al. (2020); aljg – anterolateral jugal groove; apc – anterior parietal crest; asc – anterior supraoccipital crest; ascc – channel for the anterior semi-circular canal; aVc – anterior orifice of the Vidian canal (anterior carotid foramen); boc – basioccipital crest; bot – basioccipital tubercle; bpp – basipterygoid process; bsc – basisphenoid crest; cc – cavum capsularis; c.cerf – cerebral carotid foramen; cof – common foramen; de.bcr – basicranial depression; djpp – dorsal jugal process of the parietal; efp – external frontal process; fn – frontal notch; fpp – frontoparietal process; ifp – internal frontal process; lbop – lateral basioccipital process; lep – lateral extension of the parietal; mbop – medial basioccipital process; mc – medial crest; mfp – medial frontal process; oc – occipital crest; ot – occipitocondylar tubercle (basioccipital condyle); pa – ascendens process; pc – parietal crest; pf – parietal foramen; pfs – pituitary fossa; psc – posterior supraoccipital crest; pscc – channel for the posterior semi-circular canal; pspp – parasphenoid process; ptc – pterygoid crest; pVc – posterior orifice of the Vidian canal (posterior carotid foramen); p.vpr – parietal ventral process; seltu – sella turcica; smp – septomaxillary process; snf – sympathetic nerve foramen; soa – supraoccipital area; sof – supraoccipital foramen; som – supraorbital margin; svf – supraoccipital ventral flanges; tg – trabecular groove; tp – trabecular process; tr – trabecular ridge; vjpp – ventral jugal process of the parietal.

Our analysis encompassed a total of 47 N. tessellata specimens (18 specimens with disarticulated skulls, representing only the ‘Europe’ and ‘Greece’ clades and 29 of our own or museum specimens representing other clades as well). Initially, we compared two individuals (DJ4916 and DJ5813) for which we had both 3D models and disarticulated skulls. The positive outcome of this comparison prompted us to extend our evaluation and include all 47 specimens collectively within the morphospace. All specimens were adults of both sexes (Table S5). The photographs obtained from disarticulated skulls of the two clades and the final 3D model images of all seven clades and three lineages within the ‘Central Asia’ clade were transferred to *.tps using TpsUtil32 software (Rohlf 2008) and subsequently loaded to TpsDig2 software (Rohlf 2006). In this software, 21 two-dimensional landmarks (LM) on the dorsal and 17 two-dimensional landmarks on the ventral part of the braincase (Fig. S3) were digitised. After digitising landmarks on the braincases, a Generalised Procrustes Analysis (GPA) was applied to obtain Procrustes coordinates, which were used as input variables in Principal Component Analyses (PCA). The PCA was performed with the programme PAleontological STatistics (PAST3; Hammer et al. 2007). Moreover, using the Thin-plate Spline (TPS) model, deformation grids were created to visualize the dorsal and ventral shape variation of the braincase. Cooler (blue) colors indicate a higher amount of stretching, meaning less distance between the landmarks, while hotter (red) colors indicate higher compression. This model was applied to specimens associated with the first three principal components.

Unlike the molecular phylogeny of the species, we found only little morphometric variation across their entire distribution (see Supplementary data text and Fig. S1). The results of PCA for adult and subadult males from the ‘Europe’ and ‘Greece’ clades explain 84.5% of the total dataset’s variation. PC1 accounts for 53.5% of this total variation and exhibits the highest loadings for snout-vent length (SVL; Fig. 1C). PC2 contributes an additional 31.0% to the overall dataset’s variation and is most strongly associated with head width (HW). Similarly, PCA for adult and subadult females in both clades explains 87.7% of the dataset variation. The most significant differences are found in PC1 (55.9%), which is primarily linked to snout-vent length and PC2 (31.8%), which is mainly associated with head width (Fig. 1C). We observed statistical significance (p < 0.005; p < 0.0005) in metric characteristics when comparing males and females from the ‘Europe’ and ‘Greece’ clades (Fig. 1D). The other differences were visualised in the HW/SVL ratio, where males and females of the ‘Europe’ clade display more robust (wider and larger) heads in comparison to individuals of the ‘Greece’ clade (Fig. 1E). Regarding meristic characteristics (Fig. S2A), individuals of the ‘Greece’ clade show a significantly higher number of right and left sublabials in both sexes (p < 0.005; p < 0.0005), whereby their females exhibit also a higher number of left preoculars (p < 0.005) compared with individuals of the ‘Europe’ clade. No statistically significant variances (p > 0.05) were observed in comparisons to other studied meristic characteristics between the two clades in both sexes (Fig. S2B).

Out of 448 individuals examined (Fig. 1B), comprising 269 from the ‘Europe’ and 179 from the ‘Greece’ clade, 263 (97.77%) and 176 (98.32%), respectively, possessed the so-called ‘standard’ body colouration (Fig. 2). This whole-body standard colouration ranges from olive to beige, grey and various shades of light and dark brown. Typically, black spots or blotches are arranged in four or five dorsal rows, contrasting with the lighter body colouration. Some individuals exhibited whitish or yellowish spots interspersed between the black ones along their dorsal and lateral body, although weakly spotted individuals were observed as well. In the ‘Europe’ clade, only 1.12% displayed spotless grey-olive colouration (concolorous) and an equal proportion were entirely black (melanistic). The ‘Greece’ clade exhibited even fewer unique colours, with only 1.12% displaying the grey-olive concolorous morph and only 0.56% featuring complete melanism amongst the studied populations. Differences in the colouration of the ventral and dorsal parts of the head (Fig. 3Q–T, W–Z) between the two clades showed that the ‘Greece’ clade individuals have a greater amount of black dots/speckles on the dorsal part of the head, mostly observed in juveniles and subadults and the colouration of the body more often exhibits additional pale white spots than dice snakes in the ‘Europe’ clade (Figs 2, 3U, V, Z’, Z’’).

We also noted a significant (p < 2.2e-16) and consistent presence of white labial scales across all age stages in the ‘Greece’ clade (Fig. 4A, B; see definition of all four-colour types in the Material and Methods). A Pearson’s Chi-squared test indicated a significant association (Chi-squared = 205.12, df = 3, p-value < 2.2e-16) between Type 4 labial colouration and the ‘Greece’ clade, using data gathered from online citizen-science projects (Fig. 4C). Amongst the total of 47 individuals examined within the ‘Greece’ clade, an impressive 40 (85.11%) displayed white labial colouration. In contrast, out of 720 individuals examined in the ‘Europe’ clade, only 71 (9.86%) exhibited Type 4 labial colouration, with more than 70% of these individuals being juveniles. Type 1 labial colouration was observed in 210 individuals (29.17%) of the ‘Europe’ clade but only in a single individual (2.13%) of the ‘Greece’ clade. Type 2 labial colouration was recorded in 213 individuals (29.58%) of the ‘Europe’ clade, while none of the individuals within the ‘Greece’ clade displayed this type. The most prevalent labial colouration in the ‘Europe’ clade was Type 3, observed in 226 individuals (31.39%), whereas it was found in six individuals (12.77%) within the ‘Greece’ clade. However, it remains that all four colour types were distributed throughout the entire range of the ‘Europe’ clade (Fig. 4). Finally, based on the results from colour pattern, we also anticipate that populations of the ‘Greece’ clade extend to Albania (see Fig. 4A), a distribution that has not yet been genetically confirmed.

When comparing the ‘Europe’ and ‘Greece’ clades, the first two principal components account for 42.6% of the total shape variation (Fig. 5A) in the dorsal part of the cranium. PC1 (25.5%) describes shape differences (Fig. 5B), such as slight shortening of the fronto-prefrontal suture, anteromedial displacement of parietal-jugal contact, an increase in the distance between the lateralmost points of the parieto-supraoccipital suture and craniocaudal shortening of the supraoccipital in the positive values and anterior elongation of the interfrontal suture, an increase in the distance between the anteriormost and posteriormost point of parietal-jugal contact and slight shortening of the supraoccipital in both craniocaudal and lateromedial directions in the negative values on PCA. PC2, accounting for 17.1% of the variation (Fig. 5A), is characterised by slight craniocaudal elongation of the frontal, anteromedial approaching of the anterior points of the fronto-prefrontal suture and a remarkable shortening of the distance between parietal-supraoccipital-prootic junctions in the positive values. Specimens located at the negative extreme of PC2 tend to have a slight craniocaudal elongation of the frontal bone posteriorly and a decrease the distance between parietal-supraoccipital-prootic junctions (Fig. 5B). PC3 accounts for only 9.5% of the overall variation, while samples with the highest positive values on the axis have the fronto-parietal suture shifted posteriorly and the posteriormost points of the contact between the dorsal jugal process of parietal and jugal are remarkably displaced anterolaterally, while samples with negative values have a slightly elongated supraoccipital in the posterior direction and the anteriormost points of parietal-jugal joint are displaced posterolaterally (Fig. S5C).

On the ventral part of the cranium, the comparison of the ‘Europe’ and ‘Greece’ clades showed that the PC1 and PC2 accounted for 43.6% (Fig. 5C). The PC1 (23.3%) shows changes in shifting the posterior-most point of the basioccipital and the lateralmost points of the parabasisphenoid-basioccipital suture posterolaterally and the fronto-parietal suture is markedly medially shifted in positive values. Conversely, the distance between the lateralmost points of the parabasisphenoid is greater and the fronto-parietal suture is laterally shifted in the negative values (Fig. 5D). The PC2, associated with 20.3% of the total shape variation in the dataset (Fig. 5C), is characterised by shifting the anteriormost point of the parabasisphenoid posteriorly and increasing the distance between the contact point of the basioccipital, exoccipital and prooticum in the positive values. Negative PC2 values relate to a posterolaterally shifted fronto-parietal suture from a ventral view and slight elongation of the parabasisphenoid in the anterior direction with its lateralmost points markedly shifted in the same direction (Fig. 5D). The PC3, accounting for 12.8%, is associated with shifting the anteriormost point of parabasisphenoid process posteriorly, displacement of the lateralmost points of the parabasisphenoid anterolaterally and the lateralmost points of the parabasisphenoid-basioccipital suture posterolaterally in the positive values. In the negative values, as in the PC2, specimens have a fronto-parietal suture from a ventral view shifted posterolaterally and the lateralmost points of the parabasisphenoid shifted anteromedially (Fig. S5D).

The study of 3D models of N. tessellata braincase indicates differences in morphology between the ‘Europe’ and ‘Greece’ clades present both on the dorsal (frontal, parietal, supraoccipital) and ventral (parabasisphenoid, basioccipital) parts of the cranium (Figs 6, S6 and S7) and supported by the results of PCA (Fig. 5A, C). In general, the frontoparietal portion of the braincase is slightly shorter and wider in dice snakes of the ‘Europe’ clade than of the ‘Greece’ clade. However, the most striking differences between the two discussed clades are visible in the parabasisphenoid (Figs 6G, H, S6D and S7D) which substantially contributes to the ventral base of the braincase. In the ventral view, the parabasisphenoid of representatives of the ‘Europe’ clade is posteriorly wide with strongly developed pterygoid crests directed clearly posterolaterally and protruding into distinct basipterygoid processes. The posterior openings of the Vidian canal, occurring on either side of the basisphenoid portion of the bone, are usually placed not far from the anterior openings of the Vidian canal, so the canal for the internal carotid artery passage is rather short. In dice snakes of the ‘Greece’ clade, the narrow parabasisphenoid is typical due to the very short basipterygoid processes. Posterior openings of the Vidian canal occur far from the anterior openings and, therefore, the passage for the internal carotid artery is relatively long. From the dorsal view, distinct foramina for the sympathetic nerve are situated at the posterolateral edge of the pituitary fossa close to the dorsum sellae in both clades. However, the anterior orifices of the abducens nerve (VI) occur within slit-like grooves developed laterally from sympathetic nerve foramina in representatives of the ‘Europe’ clade, whereas in individuals of the ‘Greece’ clade, those grooves are lacking and anterior orifices of abducens nerve (VI) foramina are shifted clearly posteroventrally. A summary of differences in braincase morphology between the ‘Europe’ and ‘Greece’ clades is given in Table 1.

Our results indicate that both males and females belonging to the ‘Greece’ clade have relatively longer and narrower heads unlike the shorter and wider heads uncovered in individuals of the ‘Europe’ clade (Fig. 1C, D, E). We also recorded a relatively narrower head in other phylogenetically old clades, for example, the Iran and Jordan clades. Although Deepak et al. (2023) noted that head morphology amongst aquatic natricid snakes is not consistent with the phylogeny, we suggest that a narrower and elongated head could be plesiomorphic for N. tessellata, evolved during the early evolution of this snake.

It is well documented that similar head shapes have convergently evolved in various snake groups inhabiting similar environmental conditions (Deepak et al. 2023). The head morphology of snakes is also influenced by habitat specificity and other ecological factors, such as prey preferences (Aubret et al. 2004; Brecko et al. 2011; Segall et al. 2020), whereby head morphology may change relatively quickly in reptiles (Herrel et al. 2008a). In aquatic snakes, narrower heads serve to reduce drag underwater, potentially enhancing their ability to hunt prey more effectively (Hibbitts and Fitzgerald 2005; Brecko et al. 2011; Trapp and Mebert 2011a, 2011b). The ‘Greece’ clade of N. tessellata generally differs from the ‘Europe’ clade by its comparatively more pronounced craniocaudal elongation of the braincase and its shape, which is relatively wider and sturdier in the ‘Europe’ clade. While intraspecific variation in head morphology and anatomy has been documented in N. tessellata (as reported by Mebert 1996; Andjelković et al. 2016), we propose that the distinct head morphology and anatomy observed in the ‘Greece’ clade may indeed be beneficial in its habitat of fast-flowing lotic waters (Jablonski, personal observations, but see Herrel et al. 2008b) and/or a diet focused on smaller fish prey that could minimise the costs of prey ingestion (Hampton 2011). On the other hand, the ‘Europe’ clade prefers amphibians and larger fish species in their diet (Luiselli et al. 2007; Vlček and Jablonski 2016) that are available, for example, in large Balkan lakes, such as Prespa, Ohrid and Skadar, where this clade is present. In order to confirm this hypothesis, a study of the palatomaxillar unit morphology associated with the geography is required, as the relative head and skull parameters influence the swallowing performance (Dumont et al. 2023; Segall et al. 2023).

The most significant differences in skull anatomy are observed in the parabasisphenoid (Figs 6G, H, S6D and S7D). Individuals from the ‘Greece’ clade have short basipterygoid processes and the posterior openings of the Vidian canal are positioned far from the anterior openings. However, in both clades, the basisphenoid crest, located between the basipterygoid processes, is much more pronounced compared to the one in the related extant Natrix astreptophora (Seoane, 1884) and in the extinct species Natrix longivertebrata Szyndlar, 1984 (Szyndlar 1984; Rage and Szyndlar 1986; Vasilyan et al. 2022). In contrast, the posterior margins of the basipterygoid processes cover the posterior orifice of the Vidian canal within the extant and fossil species (Pokrant et al. 2016) mentioned above. It is standard within the ‘Europe’ clade, but was exhibited only in two out of ten individuals of the ‘Greece’ clade and, thus, resembles more the condition found in other extant Natrix taxa (Vasilyan et al. 2022). From the dorsal part of the cranium, the most striking difference is the shape of the opening for the forebrain in the parietal bone where cerebral hemispheres are significantly separated from the ventrally located optic tract. It is characteristic of the ‘Greece’ clade, whereas this separation is much less prominent in the ‘Europe’ clade individuals (Figs S6B and S7B). Overall, additional examination and comparisons of other skull bones are needed, for example, a study of trophic bones such as the compound bone, which seems to be relatively longer and slender in the ‘Greece’ clade compared with the ‘Europe’ clade.

The consistently white labial scales within all ontogenetic stages of the ‘Greece’ clade versus the tendency to darker labial scales in the ‘Europe’ clade may relate to different foraging habitat preferences between these clades (Figs 3I–P and 4A–C). While we cannot state due to a lack of data whether similar distinctions in the labial scale’s colouration exist in other clades [but see some examples in Tuniyev et al. (2011) for the ‘Caspian’ clade and Amr et al. (2011), Werner and Shapira (2011) for the ‘Jordan’ clade], it appears that the ‘Europe’ and ‘Greece’ clades can be consistently differentiated by this feature. Speculatively, this colouration might serve as a form of camouflage in waters they inhabit, where the white colouration may alter the outline of the snake’s head, making it difficult for prey/predators to discern (Savage and Slowinski 1992). This adaptation could be particularly advantageous in fast-flowing, rocky waters that naturally create foam as a certain case of countershading. On the other hand, a head lacking prominent white labial scales, found exclusively in adult individuals of the ‘Europe’ clade, may be more advantageous in muddy rivers and lakes. The concept of pattern diversity resulting in camouflage (darker colouration on the top or upper side, lighter on the underside of the body) or mimicry is widespread in many animals and aquatic snakes would be no exception (Allen et al. 2013). Conversely, the possibility of bright labials serving an anti-predator function is also a potential explanation (Savage and Slowinski 1992; Wüster et al. 2004), an aspect that warrants further exploration. In general, further research combining genotypes with morphometric phenotypes can be applied to contact zones between different recognised clades of N. tessellata for a study of their character expressions of morphotypes inside, outside and respective cline across the contact zone (see Benkovský et al. 2021; Majtyka et al. 2022; Fritz et al. 2023).

The Miocene evolution in the region of today’s south-western Balkans (west of the Hellenides) likely had a significant isolating effect on the local biota. This has been substantiated by numerous studies on the historical biogeography of invertebrates (Jesse et al. 2010; Kamilari et al. 2014; Sfenthourakis and Hornung 2018), amphibians (Plötner et al. 2012; Pabijan et al. 2017; Jablonski et al. 2021) or reptiles (Gvoždík et al. 2010, 2023; Jablonski et al. 2016; Mizsei et al. 2017; Psonis et al. 2017; Kiourtsoglou et al. 2021). Currently, the Hellenides contribute to more than 60% of the overall continental Balkan herpeto-endemism (Jablonski et al. 2021) and the ‘Greece’ clade (confirmed by genomic fingerprinting) having the Miocene origin identified in the Balkans (Guicking et al. 2009), support this trend. Furthermore, both clades from continental Europe exhibit disproportionately large ranges. The ‘Greece’ clade is limited to the small area of south-western Balkans, while the ‘Europe’ clade is found in most other European regions as far as the Dnieper River, Ukraine, where it is replaced by the ‘Caspian’ lineage for the remainder of Europe as far east as the Ural Mountains (Guicking et al. 2009; Asztalos et al. 2021 and Figs 1B and 4A). These ranges are reminiscent of currently recognised endemics that diverged during the Miocene and exhibit similar distribution patterns with their closely-related species, particularly Lissotriton graecus vs. L. vulgaris (Pabijan et al. 2017), Pelophylax epeiroticus vs. P. ridibundus (Sofianidou and Schneider 1989), Podarcis ionicus vs. P. tauricus (Psonis et al. 2017), Anguis graeca vs. A. colchica (Jablonski et al. 2016, 2021; Gvoždík et al. 2023) or Vipera graeca vs. V. ursinii (Mizsei et al. 2017). At the same time, these endemics are characterised by relatively uniform morphology (e.g., Nilson and Andrén 1988; Papežík et al. 2021; Thanou et al. 2021) which, prior to DNA investigations, posed challenges for zoologists in differentiating possible cryptic taxa.

A similar situation is observed for the ‘Greece’ clade of N. tessellata, as demonstrated by morphological and osteological data presented herein. Although the ‘Greece’ clade diverged approximately 8 million years ago (Jablonski et al. 2024), it does not significantly differ in external morphology from populations representing other clades or lineages (see above and Fig. S1). While some differences may be discerned when the ‘Greece’ and ‘Europe’ clades are analysed separately (Fig. 1C–E), they are not sufficiently diagnostic to be used for taxonomic conclusions. On the other hand, variation in labial colouration could serve as a distinguishing character between these two clades (Figs 3, 4) and may have taxonomic implications, when combined with other data. Differences in cranial anatomy between the two studied clades, as presented here, are also apparent (Figs 5, 6, S6 and S7), though we advise caution due to potential morphological (phenotype) plasticity and encourage confirmation with extensive comparative material.

In conclusion, although the ‘Greece’ clade represents one of the oldest clades of N. tessellata within the species range and has been defined by mitochondrial and a limited amount of nuclear data, there are still several questions that should be investigated and addressed before any taxonomic decision can be made (see Dufresnes et al. 2023). Particularly, it is essential to delineate the contact zone between clades and assess allopatry through assessing the level of reproductive isolation or restricted gene flow using an extensive genomic dataset to detect potential historical and contemporary hybridisations known in the genus Natrix and other well-studied watersnakes, for example, Nerodia (see Mebert 2008, 2010; Kindler et al. 2017; Asztalos et al. 2020, 2021; Schöneberg et al. 2023).