Genetic and Environmentally Induced Scalation Variation in Bisexual and Parthenogenetic Lizards

1. Introduction

A switch between parthenogenetic and bisexual reproduction occurred multiple times in metazoans (Green & Noakes, 1995; Knobil & Neil, 1998). Both ways of reproduction have important advantages and not less important shortcomings. Parthenogenetic reproduction allows for the production of more offspring using fewer resources; hence, parthenogenetic populations can grow faster (Maynard Smith, 1978; Bell, 2008). Rapid reproduction allows the parthenogens to occupy intensively new areas and habitats, especially in higher altitudes with conditions suboptimal for their bisexual relatives (Glesner & Tilman, 1978; Kearney et al., 2004; Vrijenhoek & Parker, 2009). On the other hand, sexual reproduction ensures a higher diversity of genotypes in a population due to the recombination process. In case of a rapid change of environment, this may increase the chances of a population to survive in new ecological conditions.

If this basic theory is correct, one can suppose that an optimal adaptive strategy includes alternations of bisexual and parthenogenetic reproduction modes: the former helps the population to survive in extreme or rapidly changing environments, whereas the latter ensures rapid population growth and expansion when the conditions are close to species-specific optimum (Bell, 1982; Maynard Smith, 1998). Indeed, the alternation of asexual and sexual reproduction is commonplace in plants and some animals (Green & Noakes, 1995; Hempel & DeBoer, 2008; Simon et al., 2010; Pandian, 2016; Subramoniam, 2018; Tarkhnishvili et al., 2023). However, many higher animal taxa lack this adaptive strategy (Mittwoch, 1978). Parthenogenesis in vertebrates is, in general, rare. It exists in more than ten genera of lizards and snakes, some salamanders and fish. No cases are known of regular alteration of parthenogenetic and sexual reproduction in these groups.

In most cases, there are occasional parthenogens produced by bisexual parents (Lampert, 2009), which either remain sterile or return to bisexual reproduction, hence with no evolutionary or adaptive consequences. Some obligatory parthenogenetic lineages are produced mainly by distant hybridization, causing possible dysfunction of some processes during meiosis (Kearney et al., 2009). Return to sexual reproduction and producing asexual or sexually breeding backcrosses is exceptionally rare and proven or hypothesized for some fish (da Barbiano et al., 2013; Warren et al., 2018) and lizards (Cole et al., 2010; Tarkhnishvili et al., 2017, 2020a).

Although there is a lack of theoretical evidence of the unavoidable superiority of sexual breeding (Kondrashov, 1988), its prevalence in most vertebrates and other complex unitary animals is an important argument for the higher fitness of sexual breeders. However, some hybrid parthenogens appear to be quite successful, even in relatively long geological timespan. Caucasian rock lizards (Darevskia), a highly diverse group of species from the Caucasus Mountains (Darevsky, 1967; Tarkhnishvili, 2012; Freitas et al., 2024), have seven parthenogenetic “species” with hybrid origin and the time of their origin is estimated as more than 0.5 Mya (Yanchukov et al., 2022). This estimated time suggests that the existing parthenogenetic lineages survived several cycles of Pleistocene glaciation, i.e., dramatic climatic changes. Currently, they coexist with their patrilineal ancestors and are shown to be strong competitors that may even force their parental species from some habitats (Tarkhnishvili et al., 2010). Probable reasons for their success are high heterozygosity related to the hybrid origin and fast population growth potential, typical for parthenogens in general.

Phenotypic and genotypic diversity of parthenogenetic Darevskia dahli and its patrilineal ancestor D. portschinskii were recently analyzed in six river valleys where these species coexist in nature (Tarkhnishvili et al., 2020b). It was shown that the bisexual species has a higher overall diversity, both genotypic (microsatellite allele-based) and phenotypic, than the parthenogen. It was also established that genetic and phenotypic differences between the metapopulations correlate with the geographic distances in the bisexual species but not in the parthenogen. In that research, the aim was to infer the heritability of individual scalation characters, which could, in the future, help with population analysis with limited use or without the use of genetic markers. However, all analyses were conducted on a metapopulation level, and correlations on the individual level were not studied.

The system, including closely related parthenogens and their parental bisexual species, is potentially ideal for exploring another important evolutionary question: the relation between inherited (genetically determined) phenotypic variation and phenotypic variation that is environmentally induced. Despite the different biological nature of these types of variation, they are difficult to separate – due to the presence of phenocopies (Goldschmidt, 1935; Lenz, 1973) and the lack of straightforward links between gene and character due to polygenic genetic architecture of most quantitative characters, as well as pleiotropy and epistasis (Wright, 1984). Some empirical evidence suggests that selection for best-fit genetically determined phenotypes and phenotypic modification due to the direct environmental impact may have the same direction (Goldschmidt, 1949; Falconer, 1990; West-Eberhard, 2003; Burggren & Mendez-Sanchez, 2023). However, it is not clear whether phenotypic modification complements the existing genotypes-based selection, as it was supposed by Tarkhnishvili et al. (2020b), or whether there is any positive or negative association between them. To answer this question, we extended the analysis of the existing data to reveal the impact of genetic isolation and climatic conditions on phenotypic diversity, focusing on the individual-level study.

2. Materials and Methods

2.1. Studied System

The research area covered the central part of the Lesser Caucasus east and south of the riv. Mtkvari (Kura) valley and southernmost foothills of the Greater Caucasus in Georgia, where parthenogenetic D. dahli and its patrilineal, sexually reproducing ancestor D. portschinskii are commonly found in the same locality (Tarkhnishvili et al., 2010, 2020b). We sampled D. dahli and D. portschinskii from 75 localities and six metapopulation areas (Figure 1), corresponding to valleys of smaller rivers, assuming that gene flow within each area is more intensive than among the areas. 186 D. dahli and 54 females of D. portschinskii were photographed from four perspectives: pileus, right side of the head, dorsum, and anal area. Ten nominal and nine numeric characters were scored based on the collected images. Microsatellite genotypes were scored from 52 D. dahli and 27 females of D. portschinskii (Table S1).

2.2. Scoring Microsatellite Genotypes

The genotypes of the studied individuals are shown in Table S1. We scored individual lizards at five microsatellite loci: Du418, Du47, Du218, Du323, and Du215 (Korchagin et al., 2007). The maximum number of missing loci per individual was 3, and the overall proportion of missing data was 12.39%. All individuals with 2 or 3 missing loci were excluded from further analysis. Microsatellite PCR reactions were performed using the Qiagen Multiplex PCR kit with primer concentrations of 0.11–0.14 µM (approximately 0.1 mM). Thermal cycling was performed at 95 °C for 15 min, followed by touchdown of 15 cycles of 94 °C for 30 s, 56 °C for 45 s incorporating a stepwise decrease of 0.5 °C at each cycle, 72 °C for 1 min, then 22 cycles of 94 °C for 1 min, 48 °C for 90 s, and 72 °C for 1 min, followed by a final extension at 60 °C for 30 min. PCR multiplexes were developed to reduce the time and cost of genetic analyses. We amplified the five loci in two multiplex PCR reactions using the following primer combinations: Multiplex I - Du47, Du281, and Du418; and Multiplex II - Du323 and Du215. Amplicons were run on a 3130xl Genetic Analyser (Applied Biosystems), using deionized formamide and Genescan size standard LIZ 500 (Applied Biosystems Inc., Foster City, CA, USA). We amplified and scored all loci at least three times and calculated genotyping errors (allelic dropouts, false alleles, and genotype reliability), following the guidelines described in Waits et al. (2001) and Miller et al. (2002). We set the reliable change index (RCI) threshold value as 95%. We verified every individual twice for heterozygosity and three times for homozygosity since allelic dropouts are less likely to be detected at homozygous loci (Pompanon et al., 2005).

2.3. Describing Scalation

Thirty-one scalation characters described in Tarkhnishvili et al. (2020b) were analyzed from the digital images of 241 lizards from 75 localities. After removing invariant characters and those showing fixed differences between D. portschinskii and D. dahli, 19 informative traits remained (Figure 2, Table S1); all these characters vary in both taxa. No underlying information on geographic variability was used while selecting the characters.

Nine of the informative traits were numeric (a) the number of scales around the central temporal scale, SC-CT; (b) the number of the upper labials, ULb; (c) the number of the scales between the central temporal and sixth upper labial (CT-UL6); (d) the number of the scales between the tympanal and upper temporal, TM-UT; (e) the number of the scales between tympanal and seventh upper labial, TM-UL7; (f) the number of the scales between the central temporal and the third scale of the second road of postorbitals, CT-PPO3; (g) the number of the scales between the central temporal and the closest upper temporal, CT-UT; (h) the number of the scales in contact with the anal scale, SC-AN; (i) the number of the femoral pores on the left leg (FP). Other ten characters were nominal: (a) the shape of the central temporal scale (round or angular; CT-SHAPE); (b) the orientation of the CT: round, or oval – horizontally or vertically elongated (CT-ORIENT); (c) the structure of tympanal scale (TM-SHAPE; single or divided); (d) post-CT (PCT) (single or divided); (e) PO2-PPO3 – postorbital and post-postorbital (not in contact, in contact, or overlapping); (f) upper labial 3 (UL3; single or divided); (g) CPA-SHAPE – shape of the CPA – absent, small, large, broad; (h) prefrontal (PRF; side smooth, side angled, scale between the prefrontals); (i) interparietal (IP; parallel or angled lines, concaved, broken lines); (j) interparietal thickness (IPT).

2.4. Data Analysis

A comparison of genetic and phenotypic characterizations of the same individuals collected from the sampled localities was performed using a reduced sample of 54 D. dahli and 27 D. portschinskii, for which both phenotypic and genetic data were available. These individuals were also used for the analysis of genetic differentiation among the localities and metapopulations. Individual STR genotypes with missing data at any one locus and phenotypes with up to two missing measurements at each nominal and numeric trait set were included in the data analyses.

We used all studied individuals for estimating the phenotypic distances among the six studied metapopulations: 186 and 54 individual phenotypes of D. dahli and D. portschinskii, respectively.

The structure and relationships between the lizard metapopulations were analyzed based on a proper measurement of dissimilarity between phenotypes and STR genotypes of individuals. A different mode of reproduction of the two studied species allows for assuming a high level of trait and allele associations for the parthenogenetic D. dahli and low or no association for sexually reproducing D. portschinskii. The most suitable metrics and approaches for analyzing variability within and among populations of organisms with high levels of trait or allele association are assignment-based methods (Kosman, 1996; Kosman and Leonard, 2007; Kosman, 2014), whereas average-based methods are more appropriate in cases of no association. Since our objective was to compare phenotypic and genetic variability in populations of D. portschinskii and D. dahli, which exhibit sexual and asexual reproduction, respectively, both average- and assignment-based approaches were employed in this study.

We considered separately the lizard phenotypes for numeric and nominal traits (9 and 10 traits, respectively). Dissimilarity of the numeric phenotypes was calculated with the mean character difference (mcd; pp. 122-123 in Sneath and Sokal 1973) after the min-max data standardization of all trait values (for each trait, $v_{stand}=(v-v_{min})/(v_{max}-v_{min})$, where $v$, $v_{max}$, and $v_{min}$ are original, maximum and minimum its values; ${0\le v}_{stand}\le 1$) to guarantee an equal weight of all traits. Differences between the ‘nominal’ phenotypes were measured with the simple mismatch dissimilarity (m).

The STR genotypes were compared assuming two models of STR evolution: the Infinite Alleles Model (IAM) and the Stepwise Mutation Model with a variable mutation rate (SMMv). Dissimilarities between the STR genotypes were calculated using the corresponding metrics: IAM dissimilarity (eq. 3 in Kosman and Leonard. 2005; eq. 6 in Kosman and Jokela, 2019) and SMMv dissimilarity (eq. 4 in Kosman and Jokela, 2019); all mentioned equations were applied for $q=2$ as for diploids.

2.4.1. Variability Within the Metapopulations.

Dispersion of individuals within each metapopulation was estimated using the assignment-based and average-based metrics ${KW=KW}_{\rho }$ and ${ADW=ADW}_{\rho }$, respectively, with regard to the relevant dissimilarity ρ between individual profiles (Kosman 1996; Kosman and Leonard 2007; Kosman 2014), where $\rho =mcd,m$for the numeric and nominal phenotypes, respectively, and $\rho =IAM,SMMv$ for the IAM and SMMv models of STR evolution.

The effective number of different individuals (ENDI) within a metapopulation of size N was estimated with the metric of functional trait dispersion (eqs. 5, 6 in Scheiner et al. 2017): $ENDI={}_{1}D\left(T,M\right)$ for the corresponding dispersions $M={KW=KW}_{\rho }$ and $M={ADW=ADW}_{\rho }$ in the case of assignment- and average-based approaches, respectively; $ENDI$ ranges from 1 to N. Since the effective number depends on an actual number of individuals sampled in a given metapopulation, the normalized version of this indicator $nENDI={}_{1}nD\left(T,M\right)$ was considered (corrected eq. 5 in Kosman et al. 2019b; eq. 3 in Sun et al. 2020). Note, values of $nENDI$ belong to the $\left[\mathrm{0,1}\right]$ interval, $nENDI$ estimates are independent of a metapopulation size and, therefore, allows for comparison of variability within metapopulations of different sizes.

2.4.2. Uniqueness of Individuals and Metapopulations.

Uniqueness of each individual phenotype (numeric and nominal) and STR genotype (for the IAM and SMMv models of STR evolution) was estimated separately for D. portschinskii and D. dahli with the singularity metric based on the corresponding dissimilarity between the phenotypes/genotypes (eqs 1–3 in Kosman et al. 2019a). Then the singularity of a metapopulation was calculated as the average singularity of all phenotypes that belong to it (eqs 7–8 in Kosman et al. 2019a).

2.4.3. Variability Among the Metapopulations.

Distance between the metapopulations was calculated using the assignment-based and average-based metrics ${KB=KB}_{\rho }$ and ${DAD=DAD}_{\rho }$, respectively, with regard to the relevant dissimilarity ρ between individual profiles (Kosman 1996; Kosman and Leonard 2007; Kosman 2014).

Differentiation among the metapopulations was estimated with the permutation test (1000 random partitions) for differentiation statistics ${dif}_{KW}$ and ${dif}_{ADW}$ (eq. 1 in Gultyaeva et al. 2020; Czajowski et al. 2021) for the assignment- and average-based dispersions $KW$ and $ADW$, respectively [see also Kosman (2014; eq. 13) and Kosman et al. (2014; p. 565)].

The effective number of different metapopulations (ENDP) for the assignment-based (KB) and average-based (DAD) distances between them and the corresponding estimates of the extent of structural differentiation of metapopulations were obtained according to eqs. 7–9 and eq. 10 in Kosman et al. (2024).

2.4.4. Association of Dissimilarities Between Phenotypes, Genotypes and Environmental Conditions.

Association of relationships between numeric phenotypes, nominal phenotypes, STR genotypes for IAM model and STR genotypes for SMMv model was analyzed separately for D. portschinskii and D. dahli using the Mantel test with the corresponding matrices of dissimilarity between individuals. Differences between individuals for all kinds of phenotypes and genotypes were also compared (Mantel test) with the corresponding changes in two environmental conditions: average summer temperature and annual rainfall.

2.4.5. Association Analysis of Individual and Climate Differences, and Comparison of Phenotypic, Genotypic and Geographic Distances Between Metapopulations.

Comparison of individual phenotype and genotype dissimilarities on the individuals and on the locality levels with differences in temperature and rainfall in exact sampling localities of those individuals was done using the Mantel test. This test was also applied to analyze association of the assignment-based (${KB}_{\rho }$) and average-based (${DAD}_{\rho }$) distances between the metapopulations of D. portschinskii and D. dahli and the distances between localities of those metapopulations (dissimilarity $\rho$ between individuals corresponded to each data type).

2.5. Software Used

The MxComp program of the NTSYSpc package (version 2.2; Exeter Software, Setauket, NY) were used for performance of the Mantel test. Dissimilarities between STR profiles were calculated using the LOCUS software. Other above-mentioned calculations were performed with the FUNCTIONAL DIVERSITY ANALYSIS TOOLS (FDAT) software. The LOCUS and FDAT packages are available at https://en-lifesci.tau.ac.il/profile/kosman. Standard statistical analyses and visualization of data were performed with Excel 16.0 (Microsoft corp., 2021) and SPSS v. 29 (IBM corp. 2023).

3. Results

It is very challenging to compare population variability (genetic or functional) of organisms with different mode of reproduction because proper approaches and metrics for the corresponding analyses are not generally the same: the assignment-based methods are suitable for asexually reproducing organisms (e.g., D. dahli), whereas in the case of sexual propagation (e.g., D. portschinskii), the average-based methods can be used for populations with completely random matting (panmixia) and no association among genes or traits. Since there is no evidence that the studied D. portschinskii populations met the latter two conditions, a comparative study of D. dahli and D. portschinskii with the assignment-based methods seems more reasonable, and, therefore, the corresponding results are reported here. Note, however, that the average-based approaches were also applied in this study and provided rather qualitatively similar outcomes (not shown here).

3.1. Variability Within each Species and Metapopulations

Genetic variability within D. portschinskii was higher than for D. dahli: estimates of the KW dispersion and the normalized effective number of different individuals (nENDI) were 0.349 vs. 0.083 and 0.339 vs. 0.072, respectively, for the SMMv dissimilarity between SSR genotypes (Table 1). Phenotypic variability based on numeric traits was also significantly higher in D. portschinskii, with the nENDI estimate 0.298 vs. 0.188. For the nominal traits, estimates of variability were nearly equal (the corresponding nENDI values 0.841 for D. portschinskii vs. 0.836 for D. dahli; Table 1). In individual metapopulations, the genotypes and numeric phenotypes within D. portschinskii were also much more variable than for D. dahli, whereas the direction of the differences between the two species based on nominal characters varied (in three out of six metapopulations nENDI estimates were higher for one of the species; Table 1).

3.2. Differentiation Among Species and Metapopulations

Differentiation among the six metapopulations of each species was also significant for all types of data (Table 2). The extent of differentiation (ED) among the metapopulations was higher in D. portschinskii for the SSR genotypes as well as for both the numeric and nominal phenotypes (e.g., ED = 0.180 vs. 0.045 for the SSR genotypes assuming the SMMv model, and ED = 0.368 vs. 0.282 for the nominal data; Table 2). For each species, differentiation among its metapopulations was more substantial for nominal characters than the numeric ones.

The relative range of variation, $\left(max-min\right)/max$, of the metapopulation-specific nENDI values based on the numeric phenotypes was higher in D. portschinskii than in D. dahli: 0.230 versus 0.183. For the nominal phenotypes, the situation was reverse: 0.169 in D. dahli versus 0.056 for D. portschinskii (follows from Table 1). For the SSR genotypes (SMMv model), the relative range of variation of the nENDI values was much larger for D. dahli: 0.68 versus 0.38.

3.3. Uniqueness of Individuals and Metapopulations

Singularity analysis of D. portschinskii did not reveal individuals of this species with an extreme genotype or phenotype except one numeric phenotype: for the SSR genotypes, the spectrum of normalized singularity estimates was ‘continuous’ without significant gaps at the top of the list of individuals ranged in decreasing order of their singularity (Table S2). Although the average and median singularity measures based on the numeric phenotypes were higher in D. portschinskii, the relative variability or relative range of the spectrum of the corresponding singularity values (e. g., in terms of coefficient of variation) was higher for D. dahli (Table 3Tables 3, S2; Figure 3). Singularities based on the nominal traits showed nearly equal values for D. portschinskii and D. dahli (Table 3).

, where $max$ and $min$ are maximum and minimum singularity estimates for the corresponding type of data; ^f $V=SD/M$ .

Both the genetic and numeric phenotypic singularities showed presence of outliers in D. dahli but to less extent in D. portschinskii (Figure 3; Table S2). Genotypic outliers in D. dahli are explained by identical or nearly identical genotypes in most of the studied individuals; even a few de novo mutations can make an individual unique and highly dissimilar from the others. However, the phenotypic singularity indicates the presence of deviant phenotypes specific to individual local habitats. Note, a large number of D. dahli individuals with very low estimates of genotypic and phenotypic (numeric) singularities (bottom of the lists in Table S2) means a high level of genotypic and phenotypic redundancy in the parthenogen.

3.4. Association Between Individual Phenotypic and Genotypic Differences

At the individual level, there was no significant association between genetic and phenotypic dissimilarities in either of the two studied species for the SMMv dissimilarity estimates between SSR genotypes (using the SMMv model of SSR evolution) and for two types of phenotypes based on numeric and nominal traits. For both D. portschinskii and D. dahli, the relationships between individuals based on numeric and nominal phenotypes did not concur (p > 0.05).

3.5. Association of Genotypes and Phenotypes with Geography and Environmental Variables

A significant association between genetic differences and geography was established in D. portschinskii, but not in D. dahli. In the former species, genotype SMMv dissimilarities of individual lizards correlated significantly with geographic distances between individual localities (the coefficient of Mantel correlation was 0.38, p < 0.001). In contrast, in D. dahli, the correlation coefficient was close to zero and insignificant.

The outcome of correlations between the phenotypic differences among the individuals versus differences in their exact localities and the environmental variables is shown in Table 4. In D. dahli, phenotypic differences calculated from the numeric data weakly albeit significantly correlated with the geographic distances among the individuals and with disparities of the average rainfall levels at the corresponding localities. On the level of individuals, correlation of the geographic distance with phenotypic dissimilarity, based on both nominal and numeric traits, was significant for both species. Simultaneously, in D. dahli, it was a weak but significant correlation between the distance based on the numeric traits and rainfall at a locality; the correlation remained significant on the locality level. In contrast, dissimilarities between the nominal phenotypes of the D. portschinskii individuals significantly correlated with the differences in the temperature (Table 4).

4. Discussion

Not surprisingly, the overall genetic variation was higher in the recombinant population of the bisexual species than in its daughter parthenogen. In the bisexual D. portschinskii, but not in the parthenogenetic D. dahli, the pairwise genetic distances among metapopulations correlated with the corresponding distances among the studied areas, exhibiting a typical isolation-by-distance pattern. Tarkhnishvili et al. (2020b) showed that the fixation index (Rst) between pairs of D. portschinskii metapopulations from eastern Georgia correlates with geography. In the present study, we revisited the hypothesis of isolation by distance, based on genetic distances between metapopulations, and also considered the dissimilarities between STR genotypes of individuals versus the distances between the exact localities of the corresponding individuals. Such segregation is typical of a population structure that has been established for a sufficiently long time (Wright, 1943; Slatkin, 1991; Weir, 1990). Considering the significant correlation between geographic distance and genetic dissimilarity among individuals in relatively close localities (a road distance of up to approximately 130 km), we conclude that migration rates between populations are low and decrease rapidly with distance.

Genetic variation in the parthenogenetic D. dahli was low. No correlation was found between geographic distances and genetic differences, neither for metapopulations nor for separate individuals. The complete absence of association of this genetic variability pattern of D. dahli with geographic separation may result from the recent expansion of this form throughout the studied area from a single center (Slatkin, 1993). This explanation aligns with considering parthenogens as opportunistic “weed” species (Wright & Lowe, 1968; Moritz, 1991; Barateli et al., 2022). Females of D. dahli lay more but smaller eggs than D. portschinskii, which is also typical for the “weed” reproduction strategy (Barateli et al., 2022).

The overall variation based on the numeric phenotypes was also higher in D. portschinskii than in D. dahli, reflecting substantial genetic component of this variation (see also Tarkhnishvili et al., 2020b). This trend was observed for the species both for the overall samples and for individual metapopulations; D. portschinskii also showed higher differentiation among the metapopulations than its parthenogenetic daughter species. The analysis based on the nominal characters did not show such strong differences, although differences among metapopulations were again significantly higher in D. portschinskii. In contrast, singularity analysis showed that the number of outliers was higher in D. dahli for both numeric and nominal phenotypes.

Different from the genetic distances, phenotypic dissimilarities among the individuals correlated with the geographic distances among the individuals’ exact localities in both the sexually reproducing and the parthenogenetic species, and even stronger in the latter. This geography-dependent variation in D. dahli is hardly determined genetically and rather reflects environmental differences among the distant localities. A supportive argument for this suggestion is significant correlation between differences in numeric phenotypes of D. dahli and average rainfall level at the corresponding localities, both on the individual and locality level. In D. portschinskii, correlation with the annual rainfall did not exist (but a weak correlation with differences in the average summer temperatures at the corresponding localities was revealed).

These results lead us to the conclusion that the environmentally induced variation in D. portschinskii is at least not higher than in the parthenogenetic form, and even conversely, phenotypes of the latter may even stronger change under the effect of some environmental variables, e.g., humidity.

Considering the minute genetic variation in D. dahli, the correlation of numeric phenotypes with local climates should result from phenotypic plasticity. Therefore, it appears that the humidity affects an individual developmental pathway (e.g., egg development) in the parthenogen more than in its parental sexually reproducing lizard. The question remains whether this high developmental plasticity may increase the adaptability of D. dahli to changing conditions, compensating for the lack of recombinant variation. The adaptive power of a broad norm of reaction is well-known (Jablonka & Lamb, 1995; Pfenning et al., 2010), although it is not universal. For instance, intensive UV radiation increases the intensity of pigmentation in Cladocerans, but this can potentially harm the population, making individuals more vulnerable to predation by fish (Scoville & Pfrender, 2010). High phenotypic plasticity may reflect a misbalance and weakening of developmental homeostasis (Nijhout & Davidowitz, 2003; Fusco & Minelli, 2010). In our case, numeric characters are associated with the number of scales in some parts of the temporal and anal areas. At first glance, this should have, if any, a minute impact on the survival of individuals. On the other hand, Sakich and Tattersal (2021) showed that reduced scalation in bearded dragons (Pogona vitticeps) leads to faster dehydration; this may decrease viability in warmer habitats. Scalation characters may also correlate with more important traits, such as the body size of juveniles. Remarkably, a higher phenotypic plasticity of parthenogens versus their sexually reproducing relatives is not universal, even in lizards. Kearney & Shine (2004), who studied the effect of incubation temperature on the phenotype (including scalation) of sexually reproducing and hybrid parthenogenetic geckoes (Heteronotia), observed the effect opposite to our results: reduced phenotypic plasticity of the hybrid parthenogenetic lineage.

In any event, the higher phenotypic plasticity of a parthenogen requires an explanation or, at the very least, a hypothesis underlying this fact. The first, mechanistic explanation is developmental instability caused by additive genetic variation (Leary et al., 1985; Markow, 1995; Janick, 1999). All parthenogenetic Darevskia originated from hybrids of evolutionarily distant bisexual species (Murphy, 2000; Tarkhnishvili et al., 2020a; Freitas et al., 2022). The estimated divergence time of the two parental species of D. dahli, D. portschinskii and D. mixta exceeds 10 Myr (Garcia-Porta et al., 2019; Murtskhvaladze et al., 2020). We should expect multiple differential changes in genomic architecture aggregated during the period of divergence between D. portschinskii and D. mixta. These genome transformations may cause dysfunction in multiple genes (Burton et al., 2013), including those responsible for the normal meiosis process, as well as other features, such as the control of embryonic development in the hybrid. It has been previously demonstrated that certain developmental malfunctions occur in these forms, including abnormal phenotypes (Darevsky, 1960) and twinning during egg development (Barateli et al., 2025). In contrast to our results, the low phenotypic variation in hybrid parthenogenetic Heterodontia can be explained by a much more recent split among their parental species compared to the parthenogenetic Darevskia (Fujita et al., 2010).

An alternative (or, rather, supplementary) explanation is that the phenotypic plasticity, when it does not lead to the development of abnormal phenotypes, may increase the adaptive potential of D. dahli. Multiple studies are exploring the phenotypic plasticity of reptiles, which is dependent on incubation temperature (Deeming, 2004; Noble et al., 2018; Booth, 2018). Most of these studies show little effect of the incubation temperature on morphological traits such as body size or relative limb length. In contrast, the impact of temperature on egg development, growth rates, performance, and survival is substantial (Noble et al., 2018). None of the 175 studies summarized in Noble et al. (2018) described the phenotypic plasticity of scalation characters. Variation in scalation characters in D. portschinskii aligns with this general finding, suggesting that phenotypic plasticity (dependence of the scalation traits on local climate) is closely related to the genetically determined variation expressed in the differences associated with isolation by distance. However, the situation in parthenogenetic D. dahli is the opposite, due to high phenotypic plasticity, which is supposedly related to its hybrid origin and additive genetic variation.

The variability of individuals based on the nominal traits was higher in D. dahli than in D. portschinskii, and the former species had more deviant phenotypes. Since the genetic-based variability within and differentiation among metapopulations associated with isolation were higher in the bisexual species, phenotypic plasticity makes a critical contribution to the phenotypic variability of D. dahli. A significant association of variation of numeric traits in this species with the temperature of a locality suggests that the phenotypic plasticity of this parthenogenetic form is directed by a temperature significantly stronger than that of D. portschinskii.

Pfenning et al. (2010) suggested that phenotypic plasticity can significantly contribute to adaptive evolution. The increased fitness can be achieved by canalizing a particular developmental pathway (genetic assimilation or genetic accommodation) and, conversely, by widening the reaction norm. If parthenogens can evolve a broader relative reaction norm, and their phenotypic plasticity leads to the development of the most singular outlier phenotypes, this could compensate for the lack of rapid genotypic variation and, hence, the population’s inability to adapt rapidly to changing environmental conditions through natural selection. In this case, the evolution of broadening reaction norms does not require substantial genetic change and can be driven by simple dominant alleles that can arise during the mutation process.

The mechanism of adaptation through phenotypic plasticity in parthenogens remains unclear and requires detailed genomic studies in the future. However, the survival of the parthenogenetic species throughout multiple glacial cycles (Freitas et al., 2022; Yanchukov et al., 2022), as well as their ability to outcompete their sexually breeding relatives, suggests that the broadening of reaction norms can be a highly efficient adaptation mechanism.