Phenotypic differentiation of the slow worm lizards (Squamata: Anguis) across their contact zone in Central Europe

Material

For this study we analyzed 326 slow worms in total, 88 of which belonged to Anguis fragilis, 156 to A. colchica, and 82 to the group we refer to as “slow worms/individuals from the hybrid zone” or shortened as “HZ slow worms”. All individuals originated from the Central European countries Czech Republic and Slovakia and altogether they were collected at 140 localities (Fig. 1A, Table 1). A large portion of the material we obtained from museum collections (see Suppl. Information), however we also used uncatalogued specimens originating mainly from road-kills. Live individuals were individually marked to avoid pseudo-replications caused by taking data from re-captures (for the method description see Winne et al., 2006) and released at the spot of their capture right after the morphological data was taken. The sex of each individual was detected either by visual confirmation of hemipenes displayed by some live males after a careful palpation, by dissection of the tail base in freshly killed individuals, or endoscopically by probing the anal sacs (significantly shorter in females). Permits to collect the data were provided by the Ministry of Environment of the Slovak Republic (No. 9303/2009-2.1/jam and 4145/2011-2.2) and Ministry of Environment of the Czech Republic (CR: MZP/2018/630/2449).

Taxonomic assignment and genotyping

Our a priori taxonomic assignment was based on the information obtained from an analysis of the distribution of the slow-worm genotypes based on three genes—mitochondrial ND2 gene (mtDNA) and phased gametic haplotypes of nuclear genes RAG1 and PRLR (nDNA). For amplifications we used the same protocol and primers as described (Gvoždík et al., 2010; Gvoždík et al., 2013) and for the RAG1 gene we used the amplification primers R13 and R18 characterized in Groth & Barrowclough (1999; see also Gvoždík et al., 2021). Due to practical limitations, it was not possible to genotype a portion of the material used for the morphological analyses (e.g., long-term fixed material). Therefore, we constructed a detailed map of the slow-worm distribution based on the analysis of the molecular-genetic markers and assigned the material for morphological analyses comparing the distribution of the localities with the distribution of haplotypes of Anguis fragilis, A. colchica, and their hybrids (Fig. 1B). Due to this approach, our group of slow worms from the hybrid zone does not necessarily represent only the true hybrids, but in fact it may include a few individuals of A. fragilis and/or A. colchica of non-hybrid origin as well. Hybrid zone of A. fragilis and A. colchica has not been sufficiently characterized yet, but hybridization in slow worms has recently been described from the Czech Republic and Hungary based on the same genetic markers that we have been using in our study (Szabó & Vörös, 2014; Gvoždík et al., 2015). In this study we consider hybrids either (1) individuals with incongruent mtDNA and nDNA, i.e., combination of mtDNA of one species and nDNA of another species, and/or (2) individuals with interspecific heterozygotic combination within a nuclear marker, i.e., with alleles of both species in a particular nuclear locus. All hybrids were found to have A. fragilis mtDNA, but different hybrids presented five different combinations of the nuclear DNA allelic composition (Table 1).

Phenotypic characters

In the studied slow worms, we evaluated 10 metric, 11 meristic (scale numbers) and six categorical characters (for a complete list, definitions and abbreviations see Table S1). All metric variables were taken to the nearest 0.1 mm, with the exception of the snout-vent length (SVL) and the tail length (TL), both of which were taken to the nearest 1 mm. The categorical characters include the type of the prefrontal scales position, distinctiveness of the ear opening and four characters of the pattern and coloration (Figs. 2 and 3, Table S1). The coloration characters were subjectively scored based on the extent or intensity of the trait, where the lowest score meant the state of the character most resembling the character state in juveniles—thus showing the level of ontogenetic shift of the character. Due to the structure of our dataset (missing data for some traits/individuals due to preservation, injuries or autotomy) and in order to avoid decreasing the sample sizes and the analyses robustness, we decided to analyze all different types of traits, i.e., measures, scale numbers, and categorical data, separately. All data were taken by the same person (NB under supervision of DJ).

Statistical analyses

We used multivariate as well as univariate statistics to test for the differences among the three groups (i.e., Anguis fragilis, A. colchica, and HZ slow worms), with males and females treated separately. Beside the tests comparing all three groups, we also ran all analyses with only A. fragilis and A. colchica to prevent bias caused by the imprecisely defined group of slow worms from the hybrid zone. In cases when the data transformation did not improve normality or homogeneity of variance in the dataset (tested according to Keating & Hensley, 1983), we used the original input data offering better interpretability of the results (Games, 1983; Tabachnick & Fidell, 2007) and we used non-parametric tests for datasets not conforming to the assumptions of parametric tests. Exploratory analyses did not show differences between the data obtained from the live or fixed material therefore these were pooled in all subsequent analyses.

We compared SVL with analysis of variance (ANOVA). The tail lengths of the slow worms with intact tails were tested with analysis of covariance (ANCOVA) with SVL as co-variate. Head dimensions and scale numbers were analyzed using multivariate analysis of variance (MANOVA) or covariance (MANCOVA), with the exception of the subcaudal scale numbers, which were compared using the Kruskal-Wallis non-parametric test due to low sample numbers and non-normal data distribution.

The individual characters were subsequently compared using post-hoc tests, with the Bonferroni-corrected level of significance. For the analyses using large numbers of individuals and homoscedastic data we used REGWQ, Tukey post-hoc tests, Gabriel’s or Hochsberg’s GT2 post-hoc tests for uneven number of individuals, while the Games-Howell post-hoc test was always employed for heteroscedastic data (Rafter, Abell & Braselton, 2002). To reduce the multivariate data and identify which variables contribute to the observed variation the most, we also performed correlation matrix-based Principal Component Analysis (PCA), with metric and meristic characters treated separately. Multicollinearity in the metric dataset was reduced using residuals of the measures from their regressions on SVL instead of the original measures (Adnan, Ahmad & Adnan, 2006a; Adnan, Ahmad & Adnan, 2006b; Aguilera & Escabias, 2008; D’Ambra & Sarnacchiaro, 2010). To maintain potentially underlying relationship between the newly obtained components as well as to simplify their interpretation, we applied an oblique rotation of the components (Everit & Dunn, 2001; Jennrich, 2002; Costello & Osborne, 2005; Tabachnick & Fidell, 2007). Factor loadings were interpreted based on both pattern and structure matrices (Graham, Guthrie & Thompson, 2003; Henson & Roberts, 2006). In addition to the PCA with non-restricted component number, we also ran PCAs with the number of components reduced to one and plotted the obtained scores against the geographical longitude (since the HZ has north-to-south orientation in the studied region) to better visualize the geographical pattern within morphological variation. We also used Discriminant Function Analysis (DFA) to predict membership of individuals within the clusters and membership of slow worms from the hybrid zone to one of the parental species. As omitting cross-validation in the DFA model mainly reflects complexity of the dataset, we used leave-one-out cross-validation to prevent model overfitting and to improve its predictability (Lance, Kennedy & Leberg, 2000; Hawkins, 2004; McLachlan, 2004). Similar to PCA, we removed the effect of the overall body-size of analyzed individuals by using residuals from regressions of an actual measure on SVL in all DFAs. Due to the relatively high level of the tail autotomy and regeneration in slow worms (67.7% of individuals in our entire dataset), the tail length and subcaudal scale number were excluded from the multivariate analyses to avoid reducing the number of individuals and biasing against older/younger individuals and males/females (which show slightly different frequency of the tail loss; see the Supplemental Information). The frequency of the occurrence of all categorical variables was tested with Pearson χ²-tests. In analyses of coloration characters, however, we also applied the log-linear analysis which controls for the correlation between the analyzed variables (Tabachnick & Fidel, 2007; Parizanganeh et al., 2011). All analyses were performed using SPSS 17.0 (Chicago, IL).

To further investigate morphological differences across the ca. 655-km west-east transect of A. fragilis and A. colchica contact/hybrid zone, we fitted five cline models including trait interval fixed to the observed values and five combinations of fitting tail (none fitted, both tails, mirror tails, left only, right only) for selected individual phenotypic traits and PC scores of measures and scales (see above) using R package HZAR (Derryberry et al., 2014). Since the analysis requires localities with multiple individuals, we modified our dataset by decreasing the geographic resolution of localities to 0.1 degree of latitude and longitude. Convergence of the models was tested using three independent runs of each model keeping the original seeds while switching to the new seed channel and default settings of chain length, burn-in, and thinning. The best-fitting model for each trait was selected based on the lowest AICc (i.e., AIC score corrected for small sample size) implemented in HZAR package. Cline center and width were extracted from the best-fitting model for each analyzed trait. All analysis were performed in Rstudio version 1.1.453 (R Studio Team, 2015) using R software (R-Core-Team, 2020).

Anguis fragilis and A. colchica are morphologically differentiated

The results of our study show that Anguis fragilis and A. colchica are morphologically differentiated. This is in accordance with the differentiation on the genetic level, which clearly shows both species as two separate phylogenetic lineages (Gvoždík et al., 2010; Gvoždík et al., 2013). Within the genus Anguis, A. cephallonica appears to be the most morphologically divergent lineage (Grillitsch & Cabela, 1990), whereas the differentiation between A. fragilis and A. colchica is less conspicuous. However, the species differ more from each other than A. fragilis does from A. veronensis (Gvoždík et al., 2013). Although in our previous work (Gvoždík et al., 2010), we employed the concept of the genetic species (Baker & Bradley, 2006), recently discovered morphologically differentiated phylogenetic lineage and species A. veronensis was justified as an evolutionary species following the definitions of Simpson (1951) and Wiley (1978). Given that we supplement our molecular-genetic differentiation data with clear morphological differences between A. fragilis and A. colchica, both these taxa represent evolutionary species as well.

Morphological differentiation between these two taxa, considered either forms or subspecies of a single species in the past, has been debated for several decades (e.g., Wermuth, 1950; Voipio, 1962; Lác, 1967; Dely, 1974; Džukić, 1987). Previous studies suffered from the lack of evidence allowing to group individuals based on independent characters, such as genetic markers. Therefore, there was a potential risk of misidentification, mixing up the taxa, or including unrecognized hybrids into analyses. Our comparison taking advantage of the analyses of independently grouped individuals based on the genetic information found differences in all types of the studied characters including measures, scale numbers and arrangement, and coloration. However, we also show that there is no single morphological or coloration trait that is exclusively exhibited by either species and could be used as a reliable discriminatory identification character.

Patterns of morphological differentiation of Anguis fragilis and A. colchica

Slow worms from the hybrid zone are more phenotypically similar to Anguis fragilis than to A. colchica

We also analyzed morphology of the slow worms originating from the hybrid zone of Anguis fragilis and A. colchica (Fig. 1). Its course roughly corresponds to the region where the contact zone was suggested by the previous studies (Wermuth, 1950; Lác, 1967; Dely, 1981; Szabó & Vörös, 2014; Gvoždík et al., 2015).

The analyses showed that slow worms from the hybrid zone are morphologically intermediate between A. fragilis and A. colchica but showing stronger resemblance to A. fragilis than to A. colchica (Figs. 2–7; Tables 3, 4, S3–S8). This is well illustrated by a huge overlap in virtually all types of characters—measures, scalation, and coloration. Although we cannot exclude the possibility that this pattern represents an artifact and that the similarity is caused by a higher percentage of individuals with a higher portion of the A. fragilis genome or possibly A. fragilis of non-hybrid origin, it could also have arisen through selective pressures favoring the fragilis over the colchica phenotypes in the hybrid zone. A detailed genetic study linking the phenotypes with genotypes on an individual level would allow for testing these hypotheses.

Phenotypic variation shows abrupt, not gradual clines

For the variation of most of the discussed phenotypic traits, clinal pattern was suggested by several authors (Wermuth, 1950; Voipio, 1962; Lác, 1967; Dely, 1974; Dely, 1981). According to this view, an East to West cline could be observed in decreasing numbers of scales around the body and frequencies of dorsal spotting and external ear opening. In the same direction the relative frequencies of the type C vs. type A (no contact vs. broad contact) of the prefrontal scales position was supposed to change within the entire range of the genus Anguis. However, relatively sharp differences in these characters revealed by our analyses from the material originating from the contact zone indicate rather abrupt than gradual variation (Figs. 6 and 7). In support of this the HZAR analysis revealed that particularly the scale numbers around the mid-body and anterior to the tail (SCR2, 3), both strongly significant interspecifically different traits, show very steep clines with extremely narrow widths (0.7 and 1.4 km, respectively). The clines of the other two types of scales surrounding the body (SCR1 and 4) are more gradual, reflecting the general pattern observed when analyzing PC scores of overall scalation and measures. The projected widths and centers of all clines are consistent and with relatively small variation, always well under 100 km and with less than 45 km shift of the approximated cline centers. This corresponds to the revealed genetic structure and supports the hypothesis that both taxa form unique entities not only genetic (Gvoždík et al., 2010), but also morphological with reproductive-isolation mechanisms acting between them and preserving the species integrity (Barton & Hewitt, 1985; Arnold, 1997).

Divergent evolutionary histories of Anguis fragilis and A. colchica explain their phenotypic differentiation

Traditional hypothesis on the evolutionary history of the genus Anguis (or previously the species Anguis fragilis sensu lato) based on the analysis of phenotype suggested the importance of separate Pleistocene glacial refugia for differentiation between A. fragilis and A. colchica (Voipio, 1962; Lác, 1967; Dely, 1974). Anguis fragilis was believed to disperse northwards from a refugium in the Iberian and/or Italian Peninsula, while the refugium of A. colchica was supposed to be in the Balkan Peninsula and/or Anatolia and the Caucasus-Caspian region. Previously we showed that molecular-genetic variation does not fully correspond with this view (Gvoždík et al., 2010; Gvoždík et al., 2013; Gvoždík et al., 2021; Jablonski et al., 2016; Jablonski et al., 2017). It seems that A. colchica really had multiple glacial refugia within the Balkan Peninsula as speculated before, but the distribution of mtDNA haplotypes suggests that also A. fragilis might have survived the Pleistocene glacial periods in multiple refugia in the northern Balkan Peninsula and spread later to central, northern, and at least partly to western Europe (Jablonski et al., 2016; Jablonski et al., 2017; Gvoždík et al., 2021). Furthermore, the Italian Peninsula is inhabited by a different species, A. veronensis (Gvoždík et al., 2013). More importantly, the history of the separate species is older than the Pleistocene glaciations, and it is more likely that the lineages originated during the Late Miocene or Early Pliocene (Gvoždík et al., 2010; Lavin & Girman, 2019). The Pleistocene climatic oscillations are presumably only responsible for shaping their intraspecific variation (Gvoždík et al., 2013). As it remains unclear what drove the genetic divergence and speciation in slow worms, we do not know if the divergent morphologies arose during the process of speciation or resulted from the subsequent processes during the Pliocene and Pleistocene. In this scenario, the Pleistocene climatic oscillations causing population fragmentations could have played an important role in shaping the phenotypic divergence and variation. The multiple refugia of both species, despite a possibility that both species could have been surviving in the Balkan Peninsula (up to eight hypothetical refugia were identified within the northern Balkans and Carpathians; Jablonski et al., 2016), were presumably characterized by different environments. As a consequence, the populations surviving the glaciations in restricted ranges might have been exposed to variable climatic conditions. The slow worm phenotypes initially characterized by a variation largely overlapping between the lineages and/or by developmental plasticity, could have been canalized by adaptation to these conditions followed by genetic fixation of the acquired traits. In the case of lack of selection acting on a particular trait, the variation could have been lost due to genetic drift in the populations of the reduced size. Both these scenarios could have resulted in the observed morphological differentiation of the species.