Contact zone of slow worms Anguis fragilis Linnaeus, 1758 and Anguis colchica (Nordmann, 1840) in Poland

Sampling and species identification

The materials used in the study (n = 251) originated from two sources: field-collected individuals from April 2015 to August 2017, representing 35 populations of Poland, and specimens from museum collections (79 populations), representing mostly neighbouring areas in Europe (Table 1; Fig. 1). The museum specimens were taken for morphological examination only. Species affinity of collected material was set in a three-way approach. Individuals from the field were identified genetically using two molecular markers: (1) mitochondrial DNA gene NADH dehydrogenase subunit 2, (ND2) (n = 89) (Gvoždík et al., 2013) and (2) nuclear DNA gene RAG1 (n = 71) (Szabó & Vörös, 2014). At first, the individuals were classified using the BLAST tool on the ND2 sequence and then verified according to the results of phylogenetic analyses. The remaining part of the field material was identified as A. fragilis vs. A. colchica based on distributional criterion (Jablonski et al., 2021). Individuals from grey zone of incomplete molecular identification (by less than two markers) were assigned to GZ group and as such subjected to morphological analyses (Table S5A, Fig. 2). Moreover, we checked the species identification using standard morphological characters for European slow worms, including scale rows in the central part of the body and ear opening presence (Table S1A). The sexes of field-collected individuals were determined by dissections, and museum specimens were determined based on coloration pattern and other sexually dimorphic traits (Juszczyk, 1987; Sos & Herczeg, 2009).

Individuals were collected in the field (Poland) thanks to the permission of the General Director of Environment Protection in Poland (No. DZP WG. 6431.02.4.2015.JRP). The number of 78 individuals were euthanized according to the protocol of Conroy et al. (2009); no in vivo experiments were performed. In the first step of euthanasia, 1% MS222 solution was injected into the specimen’s coelom. In the second step, after the body lost the righting reflex completely, the 50% MS222 solution was injected. The injection volume was adjusted to the specimen weight. Next, secondary physical euthanasia by removing internal organs was performed (obtained organs were used for dietary and parasite infection survey). Details on euthanasia protocol, including solutions preparation, timing, and outcomes are in Conroy et al. (2009). The used specimens were deposited in the Museum of Natural History of the University of Wrocław.

Molecular laboratory procedures and phylogenetic analyses

Total genomic DNA was extracted using GeneMATRIX TISSUE extraction kits following the manufacturer’s protocols (EURX). Two protein-coding gene fragments were amplified: an 732 bp fragment of the mitochondrial NADH dehydrogenase 2 gene (ND2), and a 1010 bp portion exonic sequence of the nuclear recombination-activating gene 1 (RAG1). The protocol described by Szabó & Vörös (2014) was used.

PCR products were cleaned using a PCR/DNA Clean Up Purification Kit according to the manufacturer’s protocols (EURX), then secondary PCR with forward and reverse primers and sequencing was done by Macrogen Inc. (Amsterdam Netherland; http://www.macrogen.com). Obtained trace files of mitochondrial and nuclear sequences were automatically assembled in CodonCodeAligner (CodonCode Corporation, http://www.codoncode.com) with limited manual correction according to the trace file.

Assembled sequences were aligned in MEGA X (Kumar et al., 2018) using the default settings for gap opening and extension penalties. The same program was used to estimate genetic diversity and uncorrected p-distance. All samples were translated into amino acids, which revealed no stop codons.

The dataset for ND2 (n = 89) was supplemented with selected and previously published sequences from Genbank of A. fragilis (21), A. cephallonica (9), A. colchica (19), A. graeca (8), A. veronensis (9) and 1 Pseudopus apodus as the outgroup species (Table S3A) (Gvoždík et al., 2010; Gvoždík et al., 2013; Szabó & Vörös, 2014; Thanou, Giakas & Karnilios, 2014; Jablonski et al., 2017; Jandzik et al., 2018).

The dataset for RAG1 (n = 72) was supplemented with Ophisaurus attenuatus (AY662602) as the outgroup species (Townsend et al., 2004) (Table S4A).

The approximate models of sequence evolution were estimated using Partitionfinder2 software (Lanfear et al., 2017) to find the best partition model-based results of AICc. Phylogenies were constructed using Bayesians inference (BI) performed in MrBayes (Huelsenbeck & Ronquist, 2001) and maximum likelihood (ML) performed in PhyML 3.3 (Guindon et al., 2010) with Shimodaira-Hasegawa approximate likelihood-ratio test support of branches measured (SH-aLRT) (Anisimova et al., 2011). Due to limited DNA evolution models in MrBayes, the following models were used according to the Partitionfinder2 results: ND2 HKY + G for first position of the codon, GTR for second position and GTR + I for third. The same protocol was used for ML analyses. For RAG1, for all tree codon positions model GTR + I+G was used in Bayesians analyses. ML analyses were performed with K80 + I models, following the “Smart Model Selection” module (Lefort, Longueville & Gascuel, 2017).

Bayesian analyses for both markers were performed with two independent runs of Metropolis-coupled Markov chain Monte Carlo analyses. Each of the four Markov chains in temperature 0.2 ran for 700,000 generations and were sampled every 100 generations, except 25% of the first trees, which were excluded as burn-in (Hall, 2008)

Aligned sequences were collapsed into haplotypes in DnaSP 6 (Rozas et al., 2017). The same program was used to phase RAG1 samples to gametal haplotypes (PHASE module performed with default settings), and estimated numbers of haplotypes (h), haplotype diversity (hd), number of segregating sites (S), nucleotide diversity (π), parsimony informative sites (P) and proportion between synonymous and non-synonymous mutation for mitochondrial and nuclear markers of samples from Poland. Haplotype networks were constructed using the TSC method implemented in PopART (Leigh & Bryant, 2015).

New nucleotide sequences have been deposited in GenBank (Table 1).

Morphology and statistical analyses

Snout-vent length (SVL) and head proportions, the complex phenotypic characters that are known to differentiate the two studied lizards, were of special focus in this study. Along with the classical scheme of distance-morphometry, 10 characters were measured with an electronic calliper on the right side of the body, repeated three times, and a mean was used in further analyses (Tables S1A, S2A) (Kaczmarek, 2015; Kaczmarek, Skawiński & Skórzewski, 2016). As the classical technique of straight-line measurements was employed, the top-down omission of minimal and non-linear perturbations was assumed (Humphries et al., 1981). Well-shaped animals were exclusively taken for morphological comparisons and thus different sample sizes were used for each analysis (the numbers of specimens are written in parentheses each time). Adult slow worms of both sexes (with SVL above 120 mm, (Sos & Herczeg, 2009) were analysed separately to avoid the impact of clearly marked sexual dimorphism and ontogenetic development in these lizards.

Analyses of head shape were performed using transformed measurements for allometry according to the formula by Elliott, Haskard & Koslow (1995). This method was used to verify the putative differences in head morphology of the two species and their hybrids. Pearson’s correlation test of transformed dimensions and SVL confirmed the lack of correlations between them (Table S2B; preceded by Kolmogorov–Smirnov test).

Then, PCA was calculated on the transformed data, and a MANOVA with a Tukey post-hoc test was performed for principal component scores to compare the head shape of the two species and hybrids. MANOVA was preceded by Kolmogorov–Smirnov test and Levene’s test of variance homogeneity to fulfill test’s assumptions. Moreover, discriminant function analyses (DFA) were run to verify the correction of specimen classification to each group, i.e., A. fragilis, A. colchica sensu lato (specimens from different A. colchica subspecies were considered as an operational unit), or GZ. In all PCAs, the components were extracted based on a correlation matrix (Falniowski, 2003). The snout-vent length (SVL) between A. fragilis and A. colchica was compared with a Student’s t-test (Levene’s test of variance homogeneity was conducted before).

To describe hybrid individual morphology, identified in this study by molecular markers, in relation to the parental species, ten standard taxonomic features were evaluated including scalation, type of prefrontal scales contact, ear opening and colouration pattern (Tables S1A and S2A). The frequencies of observed variants were analysed with chi-squared or Student’s-t tests. All calculations were performed in IBM SPSS Statistics 20 (IBM Corp., Armonk, NY, USA).

Genetic differences of slow worm species and hybrids/contact zone in Poland

Two distinct diversification modes of ND2 haplotypes were seen in slow worm species in Poland. In A. fragilis as many as 12 haplotypes were identified, and only one haplotype was widespread (identical to f1; Gvoždík et al., 2010). Some rare haplotypes were sparsely represented within the occurrence range of f1. Samples originating from eight populations carried 10 new haplotypes (Table 3). All analysed A. fragilis haplotypes belong to one Illyrian-Central European haplogroup (Jablonski et al., 2016). Nine synonymous and four nonsynonymous mutations were detected.

Within A. colchica, 10 ND2 haplotypes were detected. The most frequent was c2 (Gvoždík et al., 2010). Samples from four populations carried five new haplotypes. In contrast to A. fragilis, the sequences were recognized as members of three haplogroups (I, III and IV; Jablonski et al. (2016) of the Carpathian lineage (Figs. 4–5; see also Figs. S1A and S1B). The uncorrected p of paired distances between groups equalled from I to IV 0.6%, from I to III 0.5% and from IV to III 0.8%. Greater haplotype diversity was noticed for A. colchica than A. fragilis (0.77 vs. 0.42, respectively). It is worth mentioning that specimens belonging to two A. colchica lineages (Carpathian I and III) co-occur in population 17 (south-eastern Poland). Eight synonymous and seven nonsynonymous mutations were detected.

Sympatric occurrence of ND2 haplotypes of both species was noticed in two populations (nos. 29 and 40) in the central part of southern Poland (contact zone in Upper Silesia).

For RAG1, low genetic differentiation was noticed. Eleven haplotypes were detected, eight belonging to samples of A. colchica and four to A. fragilis based on ND2 haplotype identification of the same specimens (Fig. S1A and S1B). Haplotype diversity for the dataset containing samples of both species were estimated to 0.61 (Hd = 0.058), and six polymorphic sites were detected (5 parsimony information sites), nucleotide diversity (Pi) were estimated to 0.00142, Tajima’s D statistic: 0.404. Proportion of synonymous to nonsynonymous mutation was 1:1. Thus, in both phylogenetic analyses, samples of each taxon did not form well-separated clades/clusters, however, two main groups of sequences that correspond well with ND2 classification were clearly shown (Fig. 6).

In both species, most of the specimens were homozygous (heterozygotes were recognized within samples of A. colchica from population 17 in southwestern Poland). For the two species, specific main haplotypes of wide distributions were detected: Hap_1 for A. fragilis and Hap_8 for A. colchica (Fig. 6).

Hap_1 samples form a not well-supported group (0.84/0.74) are identical to the A. fragilis RAG1 haplotype “AfR01” from Hungary (Szabó & Vörös, 2014) in the analysed sequence. Within AfR01 populations a limited presence displays Hap_4 in population 11 (new haplotype) and Hap_11 in population 18 (AfR03; Szabó & Vörös, 2014). Hap_8 in the analysed part is identical to AcR02 (Szabó & Vörös, 2014). AcR02 is widespread within all analysed A. colchica populations except the north-eastern part of Poland (populations 19 and 40). AcR02 was also found within samples in southern Poland (populations 18 and 17).

Most complex phylogenetic relationships of RAG1 haplotypes were noticed for samples from population 27; all specimens represented ND2 haplotypes, typical for A. colchica. A single specimen (C2017_4) represented an A. colchica haplotype c2 (Gvoždík et al., 2021) described from Finland and Lithuania (Hap_5). The other three specimens of not fully clear phylogeny represented new haplotypes (C2017_1; C2017_3: Hap_3, and C2017_2: Hap_4; Fig. 6)

The presence of two RAG1 haplotypes common to both slow worms was detected: AfR01 (Hap_1) and Hap_7 (new haplotype). As noticed above, AfR01 was found in most specimens classified in ND2 analyses as A. fragilis, and in two specimens from population 29 (1 ♀) and 40 (1 ♂), which belong to the A. colchica mitochondrial clade. These two specimens (TG2017_1, KK2017_1) originated from the “mitochondrial contact zone in Upper Silesia”. Thus, they fit well to be considered hybrids of both species. Their hybrid origin can also be supported by values of p distance between the two specimens to A. fragilis (0.018%), which is 14 times smaller than to A. colchica (0.25%) (Fig. 6, Table S1B). Similar p distances between the group three specimens with unclear phylogeny from population 27 and A. fragilis is about 1.8 times smaller than that between this population and A. colchica (0.17% vs. 0.32%). Thus, these three specimens are also recognized as hybrids (C2017_1 - C2017_3), or at least display gene flow between the two species (Mazovia contact zone).

The second RAG1 haplotype (Hap_7) shared by both species was found in four other slow worms, a single homozygous specimen from population 7 (ND2: A. fragilis) and three heterozygotes from population 17 (ND2: A. colchica). In two of them, their second sequence was nested within Hap_2 with A. colchica samples from population 17 (north-eastern Poland) near the main group of A. colchica sequences (Fig. 6). The position of the third sample (Lu2015_1_A; Hap_9) in the ML tree is unclear because this branch is not supported (0/0).

In contrast to specimens possessing Hap_1, the hybrid status of slow worms with Hap_7 is not promising. The long geographical distance between populations 7 and 17 (over 300 km air-distance), and from the potential contact zone makes their crossing unlikely.

Snouth-venth length (SVL)

In total, 116 specimens of A. fragilis (♂: 64; ♀: 52) and 81 specimens of A. colchica (♂: 41; ♀: 40) were used for comparison. The body sizes of A. colchica males and females (♂: M: 189 mm, SD: 27.67; ♀: M: 185 mm, SD: 37.42) are significantly larger than those of A. fragilis (♂: M: 168 mm, SD: 22.89; ♀: M: 162 mm, SD: 20.87) (♂: t[103] = −4.237, p < 0.001; ♀: t[57.383] = −3.590, p < 0.001; results of t-test with correction for variance heterogeneity).

Head shape and interspecific classification

In separate PCA analyses performed in sex groups (A. fragilis: ♂: 64; ♀: 52; A. colchica: ♂: 41; ♀: 40; GZ: ♂: 27; ♀: 13), three main principal components with a summarized eigenvalue over 72% of the total variance were analysed (matrix of components in Table S3B). Variability patterns showed by PC-graphs (Figs. 7A and 7B) were quite similar for the two sexes, as no clear clusters were formed. Nevertheless, two overlapping groups that corresponded to species affiliation occurred, mostly along the PC2 axis, and more separate in females (Fig. 7B). GZ specimens are located within the variability of A. fragilis. Such a pattern was reflected in the MANOVA results (♂: F(6; 254) = 7.005; p < 0.001; Wilk’s lambda = 0.736; ♀: F(6; 200) = 8.268; p <0.001; Wilk’s lambda = 0.642) followed by Tukey’s post hoc tests (Table S4B). The test results depicted the similarity of all males on PC1 and PC3, and a significant difference between the compared species on PC2 (p < 0.001). Moreover, a difference between A. colchica and GZ specimens was detected this way (p < 0.001). Importantly, the MANOVA assumption of variance homogeneity was not met once—for PC1 males’ comparison (F(2;129) = 7.235; p < 0.001), so this particular result should be taken with caution. MANOVA and post hoc tests performed for females depicted significant differences on the three PCs (PC1 p < 0.001; PC2 p < 0.001; PC3 p < 0.006). The Tukey’s test confirmed significant differences on all three PCs between A. fragilis and A. colchica females, and the distinctiveness of the GZ group on PC1 (vs. A. colchica, p = 0.034), on PC 2 (vs. A. colchica, p = 0.004; vs. A. fragilis, p = 0.001) and on PC3 (vs. A. fragilis p = 0.005) (Table S4B). It is worth noting that specimens of different A. colchica subspecies (A. c. incerta, A. c. orientalis, A. c. colchica and A. colchica Pontic) are distributed within variability range of A. colchica incerta in the PC1 vs. PC2 graphs (Figs. 7A, 7B).

More pronounced interspecies divergence in head morphometry was found using DFA, mainly on the first canonical function, which accounted for 95.3% and 82.7% of total variability (males and females, respectively) (eigenvalues and matrix of functions in Tables S5B and S7B) (♂: Wilk’s lambda = 0.635; df = 20; p < 0.001; ♀: Wilk’s lambda = 0.554; df = 20; p < 0.001) (Figs. 7C and 7D). The DFA performed for males showed that two morphological characters (HL1, and FL) influence the model most; for females there were five characters (HL1, HL2, HL3, Or-N and FL) (Table S8B). The highest correlation with canonical function 1 (CF1) was noticed for HW, HL1 and FL in males, and HL1, HL2, HL3, HW, OR_N and FL in females (Table S7B).

Cross-validated classification showed that A. fragilis slow worms (♂: 81.3%; ♀: 76.9%) were slightly more correctly classified than A. colchica (♂: 65.9%; ♀: 65%) to their respective taxa. Specimens from contact zones, as a rule, were mostly classified incorrectly (♂: 0%; ♀: 7.7%) and samples from populations 27, 29 and 40, which were identified as molecular hybrids, were assigned to A. fragilis in the DFA.

Morphological relationships of hybrids to parental species

The most prominent differences between A. fragilis and A. colchica among the analysed taxonomic characters (Table S1A) were observed in (1) number of scales on the ventral part of the body (V), higher in A. colchica, both males and females (M_ACmales = 133 (SD ± 5.34), M_ACfemales = 134 (SD ± 5.993); M_AFmales = 128 (SD ± 5.427), M_AFfemales = 128 (SD ± 5.376)); (2) number of scale rows around the central part of the body (SRC), the same mode of differences as above (M_ACmales = 28 (SD ± 0.6), M_ACfemales = 28 (SD ±1.01); M_AFmales = 25 (SD ± 0.715), M_AFfemales = 25 (SD ± 0.836)), (3) types of prefrontal shield contact (P), in both sexes of A. colchica the most frequent type was C (♂: 75.6%; ♀: 80%), and the least was type A (♂: 4.9%; ♀: 2.5%); in A. fragilis, types A (♂: 61.9%; ♀: 71.7%) and C (♂:11%, ♀:C = 11.3%); (4) presence of ear opening (EO), absent in A. fragilis males and females, and much more frequent in A. colchica (♂: 75.6%, ♀: 93%) (results of all performed statistical tests are summarized in Table S9B). Moreover, most of these specimens were characterized by two-side ear openings. The colour pattern is another feature that differs between the two species (Table S2A). The ventral part of A. colchica was clearly darker than in A. fragilis on the three-level scale of darkness (CV₁ = 35%, CV₂ = 47.5%, CV₃ = 15%; CV₁ = 15.1%, CV₂ = 34%, CV₃ = 50.9%, respectively). What’s more, in females, the presence of brown spots behind the head (HP₃) was noticed more often in A. fragilis (92%) than A. colchica (73.7%) (Table S9B).

Considering the five hybrid specimens identified with molecular tools (populations 27, 29, and 40), the number of scales around the central body (SRC) varied between 26 and 28. All types of prefrontal shield contact were observed, and the ear opening was present in all hybrids, at least on one side of the head (two specimens). (Table 4). The colour pattern of the ventral part of the body varied somehow: two males with the light stage (CV3) (population 27), one intermediate stage (population 40), one of blue coloration (CV₄) (27) and a female of black coloration (CV₁). Brown spots behind the head were not noticed on males (HP₁), but they were present on females (HP₂).

The hybrid zone of slow worm species in Poland

Hybrid zone is defined as a zone where genetically distinct groups of individuals meet and mate, resulting in at least some offspring of mixed ancestry (Harrison, 1993). In the Polish part of the grey zone between slow worms which extends from Bulgaria in the south up to Finland in the north (Jablonski et al., 2021), several hybrids were detected. Hybrid individuals were detected in three populations out of 35 studied: nos. 40 (Tarnowskie-Góry; 1 ♂), 29 (Kędzierzyn-Koźle; 1 ♀) and 27 (Celestynów; 3 ♂). The first two belong to the maternal genotype of A. colchica (according to ND2 haplotypes) and parental A. fragilis (main RAG1 haplotype AfR01). Both are homozygotes in nuclear DNA. Likewise, the three slow worms of population 27 belong to the A. colchica mitochondrial clade and are homozygous in nuclear marker RAG1. They carry new haplotypes that are genetically closer to A. fragilis than A. colchica, especially individuals C2017_1 and C2017_3 (Table S1B, Fig. 6). It is important to note that our findings report for the first-time hybrids of maternal ND2 haplotypes of A. colchica (Šifrová, 2017; Benkovský et al., 2021). Hybrid specimens between A. fragilis and A. colchica are quite rare; they were seldom reported in Hungary (Szabó & Vörös, 2014), Czechia and Slovakia (Gvoždík et al., 2015; Šifrová, 2017; Benkovský et al., 2021). Nevertheless, hybridization events can be quite common in these legless lizards because of the conservative karyotype (same number and structure of chromosomal sets) in this genus (Altmanová et al., 2024). Considering the microevolution of the studied lizards, the finding of homozygous RAG1 alleles in hybrid specimens, especially those from Upper Silesia, suggests they are further level hybrids. A similar situation was reported for other lizards, i.e., Iguana iguana x I. delicatissima hybrids. They inherited nuclear alleles exclusively from I. iguana whereas their mitochondrial haplotypes were specific to I. delicatissima, so they were considered secondary hybrids rather than F1 (Vuillaume et al., 2015). Interestingly, the results presented by Benkovský et al. (2021) might deliver further examples of such individuals in the genus Anguis from a hybrid zone in Bratislava. All these findings imply the retained ability for successful reproduction of A. colchica x A. fragilis hybrids and lead to the conclusion that the mechanisms of reproductive isolation between those two species are insufficient or do not work efficiently. Thus, the current formal taxonomic status of the two slow worms based on the evolutionary species concept (Benkovský et al., 2021) seems suitable for this case. Nonetheless, the limited gene flow is observed in hybrid zones in a wide array of taxa, and it is also accepted according to the contemporary understanding of various species definitions including the biological species concept (BSC) (Wang et al., 2019). Transfer of some genes can be expected within A. fragilis—A. colchica species pair, because they belong to the so-colled A. fragilis species complex, containing four closely related European slow worm taxa: A. fragilis, A. veronensis, A. colchica, and A. graeca which began to diversify around 7 Mya and currently live in parapatry (Gvoždík et al., 2023). Alternatively, it may suggest the ongoing speciation of the two forms. Moreover, it is necessary to consider other explanations for the obtained results. First, a singular mutation could change a specific sequence and mislead the result because of low diversity within RAG1. It might explain the divergence observed in specimens with Hap_7. Second, ancestral polymorphism in RAG1—the existence of a common haplotype of a progenitor inherited by the two taxa because of incomplete segregation of the new evolutionary lineages (Avise, 2004). This mechanism explains the presence of A. veronensis PRLR haplotypes from Czechia (Šifrová, 2017). It is worth noting that if ancestral polymorphism or single mutations are reliable explanations for the observed RAG1 diversity, its use for Anguis species identification is much more limited than currently believed (Szabó & Vörös, 2014). On the contrary, the analysis of singular nucleotide polymorphism in populations from the Baltic region proved its use for identifying Anguis species (Gvoždík et al., 2021). Our results showed that the description of genetic differentiation between species should include a larger number of individuals from each population. More numerous study samples allowed us to recognize a higher number of ND2 haplotypes in A. fragilis and in A. colchica from Poland (including rare sequences) than were identified in previous studies (A. fragilis 11 vs. 2, A. colchica 10 vs. 4; Jablonski et al., 2017). Similar results were obtained by Harca (2021): 103 ND2 haplotypes were detected (A. colchica—59, A. fragilis—44) in Czech and Slovak slow worms’ populations, based on over 1,300 samples. It is worth noting that 39 haplotypes were identified from single individuals (Harca, 2021). Three new locations of slow worm hybrids confirm the hypothesized course of the Europe-wide A. fragilis vs. A. colchica hybrid zone inside grey zone (Jablonski et al., 2021). Populations 29 and 40 (Upper Silesia) probably constitute an extension of a hybrid zone identified in Slovakia and Czechia (Gvoždík et al., 2015; Šifrová, 2017; Benkovský et al., 2021; Harca, 2021). According to the obtained distribution data of both species’ ND2 and RAG1 haplotypes, the Upper Silesia hybrid zone seems to be as narrow as 30–50 km. A similar size for the Anguis hybrid zone was found in France (Dufresnes et al., 2023), and even narrower (estimated to 11 km) in the Czech-Slovak contact zone (Harca, 2021). Putative hybrid specimens found in population 27 confirm Mazovia as the northernmost region of gene flow between the two slow worm species. Hybrid zone dynamics can be driven by more complex mechanisms e.g., level and direction of gene flow for which the explanation requires further investigation, especially focusing on natural hybrid individuals.

Morphological relationships of species and their hybrids

Although differences in size and shape of slow worms have been discussed for years (e.g., Wermuth, 1950), the detailed morphometrical characteristics of species and their hybrids with regard to sexual dimorphism are not satisfactorily elaborated. Our data have clearly shown the larger size of A. colchica for both sexes. Most authors agree that A. colchica reaches a larger average size than A. fragilis (e.g., Sura, 2018), however, the largest reported individual representing the genus Anguis is a male of A. fragilis (Zadravec & Galub, 2018). Recent studies have also proven a larger size of A. colchica than A. fragilis, but only for females (Benkovský et al., 2021). Size (and possibly shape) of the head determines the power of jaws in lizards, a trait under positive sexual selection (Verwaijen, VanDamme & Herrel, 2002). It was reported for A. fragilis that males with larger heads win combats for females (Capula & Filippi, 1998) and may have a stronger grasp of the female during courtship, as recorded for the common lizard Zootoca vivipara (Gvoždík & Van Damme, 2003). According to the aforementioned concept, A. colchica might monopolize all females available for mating, but there are no signs of outcompeting the second species from its area of distribution.

On the other hand, the multidimensional statistical analyses employed in this study (PCA, DFA) showed that the differences in the shapes of the heads of A. fragilis and A. colchica are not distinctive enough, and the classification algorithms are not efficient, especially for A. colchica (correctness of 65% at minimum). As depicted in this study, weak specificity of morphometric features of the two lizard species and their hybrids corroborates with conservative karyotypes (Altmanová et al., 2024) but contradicts evident genetic differentiation based on genome-wide nuclear DNA and mitogenomes (Gvoždík et al., 2023). Still, it must be noted that in some previous studies, such a correlation was more pronounced (Benkovský et al., 2021). The origin of slightly different results in independent studies might have various causes. First, the local environment can modify morphological variability, e.g., Benkovský et al. (2021) suggested that at least different positions of prefrontal shields between species might be caused by environmental conditions or embryonic development. Second, different analytical approaches could affect the resolution of the results at genetic and phenotypic levels. Moreover, both mentioned reasons could interplay in the definite picture of phenotypic variability (Mayr, 1969).

As the head shape of the two slow worm species is similar, but not the same it can be expected that hybrid specimens will resemble one of the parental species more or will show intermediate head proportions. The correctness of slow worm species cross-validated classification in DFA was low; however, individuals identified as molecular hybrids were closer to A. fragilis. Recently published data on phenotypic differences (combining morphometric and meristic traits) between slow worm species and their hybrids also showed that specimens from the hybrid zone in Czechia and Slovakia resemble A. fragilis more than A. colchica (Benkovský et al., 2021). However, hybrids from the Polish part of GZ represented more A. colchica-like phenotypes exclusively in meristic features. Disproportions of parents’ phenotypes in hybrid offspring is a common phenomenon for vertebrates, including examples from Squamata, e.g., Pituophis catenifer and Pantherophis vulpinus hybrids are more like one parental taxon in head shape, but intermediate in meristic traits (e.g., number of ventral scales) (Leclere et al., 2012). Complex genetic and epigenetic mechanisms control the expression of parental phenotype in hybrids (Bartoš et al., 2019), sometimes manifesting maternal dominance (Wolf & Wade, 2009). A closer similarity of all five molecular hybrids to the maternal parent A. colchica in meristic traits may result from maternal inheritance affecting the ultimate phenotype of hybrids.

It can be predicted that some future studies of the contact zones of Anguis species will drive to detect some new hybrid individuals. They seem to be key specimens for explaining the mechanisms that maintain the dynamics in the hybrid zones and play a crucial role in the Europe-wide distribution of slow worms.