Comparing morphology and cranial osteology in two divergent clades of dice snakes from continental Europe (Squamata: Natricidae: Natrix tessellata )

Comparing morphology and cranial osteology in two divergent clades of dice snakes from continental Europe


Introduction
The genus Natrix, described by Josephus Nicolaus Laurenti in 1768, currently includes five snake species widely distributed across the Western Palearctic Region (Speybroeck et al. 2020).The dice snake Natrix tessellata (Laurenti, 1768), with its type locality on the Istrian Peninsula in present-day Croatia, has the most extensive distribution, ranging from central Europe and northern Egypt to north-western China and Pakistan (Gruschwitz et al. 1999;Mebert 2011a and references therein;Mebert et al. 2013).Its phylogenetic position and phylogeogra-phy have been studied for more than two decades, revealing a high level of phylogenetic diversity across its range (Guicking et al. 2002(Guicking et al. , 2009;;Guicking and Joger 2011;Kyriazi et al. 2013;Rastegar-Pouyani et al. 2017;Asztalos et al. 2021;Jablonski et al. 2024).
Initially, nine evolutionary lineages corresponding to distinct mitochondrial clades, each identified according to their geographic distribution (Europe, Crete, Greece, Turkey, Jordan, Iran, Kazakhstan, Uzbekistan and Caucasus), were distinguished in the mitochondrial phylogeny of N. tessellata and supported by genomic fingerprinting (Guicking et al. 2009).While these clades are well-defined, a recent study re-defined the phylogeny to seven clades that diverged significantly during the past climatic fluctuations and environmental changes (Jablonski et al. 2024) with further diversification into Pleistocene lineages.Most of these clades have a Late Miocene origin, particularly those from the Middle East, which might be the subject to further taxonomic evaluation (Dufresnes et al. 2023).
In continental Europe, two clades, namely 'Greece' and 'Europe', occur in parapatry.Their divergence is estimated to have originated during the Miocene epoch, approximately 8 and 6 million years ago, respectively (Jablonski et al. 2024; Fig. 1A).Such ancient divergence implies that various complex factors have influenced this and other reptilian species' evolutionary history and are often linked to the Messinian salinity crisis, specific geomorphological developments and the current topography of the southern Balkans (e.g., Guicking et al. 2009;Jablonski et al. 2016;Thanou et al. 2023).One of these two clades subsequently expanded its range to encompass broader areas of Europe and western Anatolia (the 'Europe' clade), while the second clade ('Greece') likely evolved in the southern Balkans, west of the Hellenides (the mountain range roughly spanning from northern-central Albania to southern Greece; Fig. 1B).The precise location of the contact zone between these two clades remains unknown, but it is potentially located between southern Albania and central Greece (Jablonski et al., unpubl. data).In addition, Greece is the only European country where three evolutionary clades occur, with the third one being endemic to the island of Crete (the 'Crete' clade).Additionally, Greece is renowned for its relatively rich records of natricid snake fossils dating back to the Late Miocene (Georgalis et al. 2016(Georgalis et al. , 2017)).Potential older natricids from the area are even known from the Early Miocene (Vasileiadou et al. 2017).Moreover, the genus Natrix in Greece also includes extinct forms dating from the Neogene (Georgalis et al. 2019) and encompasses several occurrences in the late Neogene and Quaternary (Szyndlar 1991;Georgalis and Delfino 2022).
Nonetheless, the question of whether these deeply divergent clades display consistent differences in external morphology or osteology, which may also reflect their natural history, has never been investigated.The currently available morphological data from the previous studies, analysed without information on DNA variation, indicates minimal diversity in the morphology of European populations of N. tessellata.(Laňka 1973;Kramer et al. 1982;Rehák 1992;Lenz and Gruschwitz 1993), although Mebert (2011b) showed a close phenotypic relationship of mainland Greek populations with Dalmatian populations (mainly Adriatic coast of Croatia), based on the body proportions, positions of scale-row reductions and numbers of ventral, subcaudal and cephalic scales.Hence, even though varying ecological conditions across different spatial and temporal scales can exert evolutionary pressure on populations and potentially result in morphological diversification (as suggested by Ricklefs et al. 1981;Vitt et al. 1997;Tulli et al. 2011), we presently do not possess enough morphological data that can be examined in conjunction with the molecular phylogeny of N. tessellata.Additionally, only some general information exists about the cranial and postcranial differences and anatomy of the species at the inter-population level (Szyndlar 1984;Jandzík and Bartík 2004;Ratnikov and Mebert 2011;Andjelković et al. 2016Andjelković et al. , 2017;;Pokrant et al. 2016).Therefore, the primary objective of the current study is to evaluate morphological and osteological differences between the two principal N. tessellata groups of continental Europe, namely 'Europe' and 'Greece' clades.

Morphological data collection
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. tes sellata (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).

Morphological data analyses
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.

Osteology Skull preparation, micro-CT projection, landmark digitalisation and statistical analyses
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 2 O 2 ) 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).
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 Tps-Dig2 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.

Quantitative analysis of external morphology
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 snoutvent 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).

Differences in colouration and pattern of 'Europe' and 'Greece' clades
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.Dif-ferences 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.

Osteology Skull morphology
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 neg-ative 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 lateral-   most 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).

Braincase bones comparisons between the 'Europe' and 'Greece' clades
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.(Gygax 1968(Gygax , 1971) and the skin sensory organs (Walztöhny and Ziswiler 1979).However, the relatively great morphological uniformity across the entire distribution, particularly the ubiquitous olive-grey snake with black spots as the most common morph, is the reason why the species has not been taxonomically partitioned into many subspecies, unlike many other European reptiles (Mertens and Wermuth 1960).For over 250 years since N. tessellata was described, a number of local morphs have been named (e.g., Natrix tessellata var.flavescens Massalongo, 1853, Tropidonotus tessellatus var.sparsus Dürigen, 1897, or Natrix tessellata var.cyréni Sochurek, 1956; see list in Uetz et al. 2023).However, only one subspecies remained for an extended period (Tropidonotus tessellatus heinrothi Hecht, 1930), but was later synonymised (Gruschwitz et al. 1999).This contrasts substantially with the closely-related congeneric grass snake Natrix natrix (Linnaeus, 1758) complex (Fritz and Schmidtler 2020), both exhibiting similarly extensive and largely sympatric ranges from central Europe eastwards.
Earlier assessments of N. tessellata morphology, including morphometric data, pholidosis and colouration, were provided by Dürigen (1897), Schreiber (1912) and Hecht (1930).However, their sample sizes and data quality were insufficient for general conclusions.Subsequent investigations into the morphology of the species then became more regionally focused (e.g., Laňka 1973Laňka , 1978;;Kminiak and Kalúz 1983;Rehák 1992;Lenz and Gruschwitz 1993;Mebert 1993Mebert , 1996;;Göçmen and Böhme 2002;Herczeg et al. 2005;Dincaslan et al. 2011;Werner and Shapira 2011;Savasari et al. 2019), rendering it challenging to draw large-scale comparisons.The first primary morphological division was proposed by Kramer et al. (1982) who categorised N. tessellata into a Mediterranean and an East European-Asian group.Later, Mebert (2011b) published relatively extensive data on morphological characteristics from the entire range, mainly focusing on the scalation.According to his results, the number of ventrals, subcaudals, oculars and labials increased substantially, but gradually from west to east and from south to north, but the overall external morphology remained uniform.On the other hand, on the microgeographic scale, using 50 morphological characters was sufficiently powerful to detect statistically significant and consistent geographic variation amongst natural and introduced populations in northern Italy and Switzerland (Mebert 2011c).Moreover, a cluster analysis, based on 27 morphological characters (scalation and body proportions), correctly resulted in the partition of N. tessellata per sex into a clade for Italy and Switzerland and another clade for individuals from the western Balkan and eastern Mediterranean, whereas the latter clade was further split into Croatia+Greece versus Anatolia+Levant (cluster diagrams in Mebert 2011b, 2011c).In addition, Brecko et al. (2011) found that frog-eating N. tessellata populations have developed a broad-headed morph, whereas fish-feeding populations produced a narrow-headed morph.Finally, morphology appears to contrast significantly with the observed molecular phylogeny of the species (Guicking et al. 2009;Kyriazi et al. 2013;Jablonski et al. 2024) and indicates that local or regional selection is influential in morphological expression.
Similar to previous studies, our study showed low intraspecific morphological variation amongst studied populations/clades of N. tessellata, especially in the 'Europe' and 'Greece' clades (Fig. S1).While the quantity of morphological data from the eastern part of the species' range is not equivalent to that obtained from 'Europe', our analytical comparisons represent the first broad geographic scale view of the species' morphology (see Supplementary data text).Coupled with the context of molecular phylogenetic hypotheses (Guicking et al. 2009;Kyriazi et al. 2013;Jablonski et al. 2024), this opens pathways for further research.Overall, the body proportions of N. tessellata do not exhibit a clear geographic pattern in any cardinal direction, albeit individuals from Greece show relatively longer tails (TL/SVL) that are about 1.54-2.65%longer compared to dice snakes from western and central Europe, consistent with the results of Mebert (2011d).So far, the longest European individual measuring 130 cm in total length, was found on Serpilor Island in the Black Sea (Ukraine, formerly Romania; Calinescu 1931; Mebert 2011b).However, locally good conditions appeared to have produced particularly many large-grown individuals in Lake Geneva and the Caspian Sea (Mebert 2011d;Tuniyev et al. 2011).Indeed, adult N. tessellata from the introduced populations in the fish-rich lakes Geneva, Alpnach and Brienz of central and western Switzerland normally achieved greater lengths and weights than from their parental population (source of translocation) in the Maggia Valley in southern Switzerland (Mebert 2011b).In conclusion, our findings as well as previous data concur that geographic variation in metric characteristics amongst N. tessellata populations depends rather on local or regional selection than on affiliation to any of the seven clades (see Fig. S1) and, thus, should not be used as determining characteristics for clade/lineage classification.
When it comes to the variation in meristic data (see Supplementary data), the most striking differences were observed in the number of ventral scales.Individuals from the western part of the range had 8-9 fewer ventral scales compared to those from the easternmost populations in 'Europe' in this study.However, Mebert (2011b) recorded even higher mean ventral values for N. tessellata from north of the Black Sea between the Crimea Peninsula and the Volga Delta, representing a mix of the 'Europe' clade and 'Caspian' lineage.For subcaudal counts, the highest mean values (~ 75) in males were obtained from the western Balkan and (~ 67) in females from Greece (Mebert 2011b).In labial and ocular scales, there is a trend to increase from west to east as suggested by Mebert (2011b).However, our results indicate an opposite trend.For instance, the most common combination of preocular scales is 2/2 (61.4%) and for postoculars 4/4 (50.0%) within western and central European populations (Table S4).However, the bilateral proportions (ratio) of preoculars of 2 or less vs. 3 or more scales is 70.3/29.7 in the 'Europe' clade and 38.6/61.4 in the 'Greece' clade.In postoculars, the bilateral proportions of 3 or less vs. 4 or more scales is 29.6/70.4 in the 'Europe' clade and 34.4/65.6 in the 'Greece' clade.This generally confirms the trend observed by Mebert (2011b); however, the comparatively small decrease of postocular numbers in the Greek samples herein likely relates to our inclusion of a much greater proportion of Eastern European dice snakes with somewhat higher values compared to the much higher proportions of more western dice snakes from Italy and Switzerland in the samples analysed by Mebert (2011b).No differences were observed in ocular scales in the south-north directions as previously suggested by Laňka (1973Laňka ( , 1978)), Szczerbak and Szczerban' (1980) and Rehák (1992).
Across the entire distribution range, the most frequent combinations of supralabial scales were 8/8 and 9/9 for sublabial scales, whereas 10 subalabials is also reasonably common, a situation also confirmed herein for eastern and central European N. tessellata (see also Laňka 1978;Mebert 2011b).Yet, it remains that the frequencies of eight or less supralabials and sublabials is lower in the 'Europe' clade than in the 'Greece' clade, paralleling results in Mebert (2011b).Generally, geographic variation in subcaudals, oculars and labials is less pronounced than the variations observed in ventral scales (Laňka 1973(Laňka , 1978;;Rehák 1989Rehák , 1992;;Gruschwitz 1993;Mebert 1993).
Consistent with our results, overall body colouration typically ranges from grey, olive to dark brown with 4-5 rows of blotches throughout the whole species range (Mebert 2011b).We also recorded colour aberrations and morphs (melanism, concolorous; see Table S1); however, it seems that this colour polymorphism is distributed without any geographic pattern and can be locally absent or marginally present (Mebert 2011b).The syntopic triple colour polymorphism (standard spotted, concolorous or unicolour, melanistic) appears particularly common in very fish-rich water bodies, for example, Prespa Lake in North Macedonia, Beyşehir Lake in Türkiye, Caspian Sea or fish hatcheries in Syria (Dinçaslan et al. 2011 and other articles in Mebert 2011a).The 'Europe' and 'Greece' clades show fewer differences in body colouration (Fig. 2), but we observed small white spots on the dorsal and lateral parts of the body (e.g., Figs 2K, P, Q and 3Z') more frequently in dice snakes of the 'Greece' clade than in the 'Europe' clade.
Overall, these findings prompt the question: Why does N. tessellata exhibit such pronounced morphological uni-formity?One possible hypothesis could be rooted in the aquatic environment and overall natural history of the species, likely originating around the former Parathethys area (see Jablonski et al. 2024).The aquatic environment is very stable and often leads to convergent evolution of the body, probably due to adaptations related to moving in the water (Deepak et al. 2022(Deepak et al. , 2023;;Hallas et al. 2022).On the other hand, the dice snake may have developed the successful basic morphological model early on in their evolutionary history that geographically radiated and persisted to the present day.Given that N. tessellata is primarily a (semi)aquatic snake with a specialisation in fish consumption throughout its entire range, it is reasonable to assume that different populations and clades have survived under relatively similar conditions of aquatic habitats.This stable environment might result in reduced selective pressures on the colouration phenotype, thus maintaining the standard spotted morph across the entire geographic range.Although N. tessellata appears to exhibit a very successful and ubiquitous colour pattern model, genetic variation increased because of allopatric divergence and selection for more specific local traits not relating to colouration aside from some local to regional differences (e.g., frequencies of melanistic and concolorous morphs).Therefore, additional research pertaining to the taxonomy of this snake should also place emphasis on thorough investigations to identify potential other morphological distinctions amongst the defined taxa.
Hence, it would be intriguing to investigate and regionally compare the hemipenes, the body structure that displays high morphological variation amongst snakes (Andonov et al. 2017) and may have played a pivotal role in reproductive isolation and speciation in general.While brief descriptive interpretations of the internal male reproductive organ morphology of N. tessellata have been provided from Switzerland (Cadle 2011), Iran (Savasari et al. 2015), Bulgaria (Andonov et al. 2017) and by Gruschwitz et al. (1999) from an unknown locality, there is a need for comprehensive and comparative research to fill these gaps and facilitate future investigations.

Can differences amongst the 'Europe' and 'Greece' clades be explained by environmental and ecological factors?
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 hab-itat 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 andMebert 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).

Does the 'Greece' clade represent possible cryptic snake species from the Balkans?
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(Gvoždík et al. , 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' lin-eage 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(Jablonski et al. , 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 2008Mebert , 2010;;Kindler et al. 2017;Asztalos et al. 2020Asztalos et al. , 2021;;Schöneberg et al. 2023).

Figure 1 .
Figure 1.A The dated mitochondrial phylogeny of Natrix tessellata, as adopted from Jablonski et al. (2024), which served as the basis for the division of the morphological and osteological dataset; B The geographic origin of the morphologically investigated specimens of Natrix tessellata from the 'Greece' and 'Europe' clades, divided, based on the mitochondrial distribution as per Guicking et al. (2009), Asztalos et al. (2021) and Jablonski et al. (2024); C Principal Component Analysis (PCA) of the 'Europe' (green) and the 'Greece' (purple) clades, based on metric characteristics compared in males and females separately.Most of the variation in the dataset are explained by the first two axes given in %; D Boxplot comparisons of metric characteristics between both clades.Black circles represent means.Asterisks denote the degree of significance ** p < 0.005, *** p < 0.0005; E Linear regression indicating HW rates against SVL rates between males and females of both clades.

Figure 5 .
Figure 5.The position of Natrix tessellata specimens from the 'Europe' and 'Greece' clades within the morphospace from the dorsal (A) and ventral (C) parts of the cranium with the three-dimensional models and dissected braincases correlated with the PC1 and PC2 axes, along with the coloured deformation grids (B, D).Cooler colours indicate a higher amount of stretching, hotter colours indicate a higher amount of compression.Borderline specimens are in bold and the numbers assigned to the specimens correspond to those listed in TableS5.The scale bar indicates 1 mm.

Table 1 .
The most significant differences between the studied bones on the dorsal and ventral parts of the braincase observed within the 'Europe' and 'Greece' clades.