Research Article
Research Article
The mitochondrial phylogeography of the Crimean endemic lizard Darevskia lindholmi (Sauria, Lacertidae): Hidden diversity in an isolated mountain system
expand article infoOleg Kukushkin§, Oleg Ermakov|, Iulian Gherghel#¤, Svetlana Lukonina|, Anton Svinin«, Igor Doronin§, Evgeniy Simonov«, Daniel Jablonski»
‡ A. O. Kovalevsky Institute of Biology of the Southern Seas, Russian Academy of Sciences, Feodosiya, Ukraine
§ Zoological Institute, Russian Academy of Sciences, St. Petersburg, Russia
| Penza State University, Penza, Russia
¶ Case Western Reserve University, Cleveland, United States of America
# Lucian Blaga University, Sibiu, Romania
¤ Alexandru Ioan Cuza University, Iași, Romania
« University of Tyumen, Tyumen, Russia
» Comenius University, Bratislava, Slovakia
Open Access


The Lindholm rock lizard, Darevskia lindholmi, is the only member of the genus Darevskia whose range is restricted solely to Europe, representing a local endemism found only in the Crimean Mountains. In our study, we investigated the cytochrome b gene (mtDNA) of 101 D. lindholmi sequences from 65 Crimean localities, representing its entire range. We found that D. lindholmi is highly genetically structured, and its range is divided into populations belonging to three mitochondrial lineages. The Lindholm rock lizard populations inhabiting the middle part of the Crimean Mountains (further referred to as the Central lineage) are sharply differentiated from the other two lineages (the Common and the Southwestern lineages), which are present in most of the species range. The genetic distance between the Central lineage and the other two taken together is 4.6%, according to our results, suggesting that the divergence occurred during the Early Pleistocene. The narrowly distributed Southwestern lineage and the widespread Common lineage, on the other hand, are differentiated by 1%. Field observations on the representatives of the main evolutionary groups show that their ecology is also different: the Central lineage is a mesophilic and cold-resistant form, while the other two closely related lineages are more xerophilic and thermophilic. Results of the potential ranges modeling and ecological niche analysis confirm that the genetic lineages occupy different niches of the Crimea. Furthermore, the area of inhabitation of the Central lineage splits the western and eastern parts of the Common lineage range, while the Southwestern lineage is restricted along the coast of the southwestern coast of the peninsula. The long-term co-existence of deeply divergent sister mitochondrial lineages in a relatively small (circa 7,000 km2) isolated mountain system serves as a mesocosm for understanding the speciation process. Our data suggest that the Central lineage warrants further taxonomic investigation.


Caucasus, Crimea, cryptic lineage, endemism, Lacertidae, Quaternary, speciation


The genus Darevskia Arribas, 1999 includes more than 30 species of medium-sized lizards with ranges localized mainly in the Caucasus, Anatolia, and neighboring regions, including southeastern Europe and northern Iran (Arnold et al. 2007). The Lindholm rock lizard, D. lindholmi (Szczerbak, 1962), is a Crimean endemic and a member of the D. saxicola complex (Szczerbak 1962; Darevsky 1967; MacCulloch et al. 2000; Doronin et al. 2013). The type locality of this species is located in the city of Yalta on the Crimean Southern Coast (Doronin 2012a). Besides D. lindholmi, the D. saxicola complex includes at least three closely related species distributed in the western and central Caucasus: D. saxicola (Eversmann, 1834), D. brauneri (Méhely, 1909), and D. szczerbaki (Lukina, 1963) (Doronin et al. 2013). However, the species status of D. szczerbaki inhabiting the extreme northwestern part of the Caucasus is questionable and it might be a hybrid form (D. saxicola its matrilineal ancestor, and D. brauneri the patrilineal; Tarkhnishvili et al. 2016). On the other hand, some populations or taxa from southern Caucasus are also questionable with their affiliation to particular Darevskia species complexes (Tarkhnishvili et al. 2016). Currently, the assignment of D. alpina (Darevsky, 1967) to the aforementioned species complex or the D. caucasica complex is still a subject of discussion (Murphy et al. 2000; Doronin et al. 2021). Based on morphological data, D. lindholmi showed a close relationship to D. brauneri (Szczerbak 1962; Doronin et al. 2013), and moderate values for genetic distances based on cytochrome b between the Crimean and the Caucasian representatives of the complex (3.4–5.3%; Doronin et al. 2013) indicate a Plio-Pleistocene time of divergence (see Tarkhnishvili et al. 2016). The genetic distances between most bisexual Darevskia species are much higher (exceeding 10%; Tarkhnishvili 2012; Ahmadzadeh et al. 2013; Kurnaz et al. 2019; Doronin et al. 2021). The highest diversity and uniqueness of satellite DNA discovered in D. lindholmi led to the assumption of a possible hybrid origin of the species (Grechko et al. 2006). All this suggests an unclear evolution and leading to confusing taxonomy in the D. saxicola complex that requires further revisions (Speybroeck et al. 2020).

Darevskia lindholmi currently inhabits an isolated mountain system of the Crimean Peninsula, without any geographical connection to the Caucasus that is expected as the radiation center of the genus (Darevsky 1967; Murtskhvaladze et al. 2020). During most of the Neogene, Crimea existed as an island or an archipelago, although its continental connection with adjacent landmasses was established more than once during relatively short periods (Muratov 1960; Esin et al. 2018; Stovba et al. 2020; Palcu et al. 2021). The Crimean Mountains have a range of ca 7,000 km2 (25% of the total area of the Crimean Peninsula) and consist of three ranges stretching from the southwest to the northeast. The Late Pleistocene cooling in Crimea affected the biota to a more substantial extent than in the Caucasus, although no evidence of mountain glacier presence was found here (Vakhrushev and Amelichev 2001). Because of a limited scale and relatively simple topography, along with long-term connections with the adjacent continent, Crimea acted as a “melting pot” rather than refugia for biota (see Mendel et al. 2008; Dufresnes et al. 2016; Kukushkin et al. 2018; Levin et al. 2019). Thus, the level of endemism here is low, and unique genetic diversity (though only slightly different from the mainland populations) was detected only for few vertebrate species, especially reptiles (Fritz et al. 2009; Zinenko et al. 2015; Psonis et al. 2017, 2018; Jablonski et al. 2019). In addition to the Crimean member of the genus Darevskia, the only endemic reptile taxon to the Crimean Mountains is Lacerta agilis tauridica Suchow, 1927, whose morphological uniqueness was supported also by genetic data (Joger et al. 2007; Andres et al. 2014; Kukushkin et al. 2020). Lacerta viridis magnifica Sobolevssky, 1930, former endemic subspecies of Lacertidae from Crimea, was identified as an introduced, currently extinct population of L. bilineata (Kehlmaier et al. 2020). Thus, D. lindholmi, residing outside the Caucasus provides grounds for questioning its origin and evolution. It has been suggested that this species’ ancestor was isolated after the disappearance of the mountain-forest bridge that directly connected the Crimean Mountains with the north-western extremity of the Greater Caucasus during one of the cool and wet periods of the Early Pleistocene (Szczerbak 1962, 1966; Darevsky 1967). At the same time, studies on satellite DNA repeats (see Ryabinin et al. 1996) suggest that D. lindholmi is much older than the Caucasian taxa of the D. saxicola complex, whose radiation was attributed to the postglacial time. Ryabinin et al. (1996) also hypothesized that the colonization of Crimea by a D. lindholmi ancestor could be associated with its passive transmarine distribution in the Early Pleistocene. However, the question of the species origin and phylogeographic structure inside Crimea remains unanswered.

In our study, we thus investigated the mitochondrial (mtDNA) phylogeography and the niche divergence of D. lindholmi, particularly i) studying intraspecific genetic variability, and ii) assessing the spatial distribution and ecological niches of the detected mtDNA lineages to provide a biogeographical and taxonomic context.

Material and methods

Sampling and DNA extraction

Tissue samples were collected mainly by O.K. between 2011 and 2019 throughout the entire species’ range (Fig. 1). In total, 101 sequences from 65 localities were finally analysed, including 92 sequences from 60 localities that were novel for the study (Table S1). The tissue samples comprised of autotomized lizard tails, or finger phalanges fixed in 96% ethanol further processed in the lab via genomic DNA extraction using a standard salt protocol (Aljanabi and Martinez 1997).

Figure 1. 

Distribution of mitochondrial DNA haplogroups in Darevskia lindholmi. Circles indicate the sampling localities: dark-blue – Southwestern lineage, blue – Common lineage; yellow – Central lineage. The type locality of the species is indicated by a star. Geographic position of the localities are given in Table S1. The location of the Crimean Peninsula in the Northern Black Sea region and the position of the Crimean Mountains are shown on the inset maps.

mtDNA sequencing and analyses

The cyt b gene (1131 bp) was amplified using a primer pair LgLu (5’-AAC CRC YGT TGT MTT CAA CTA-3’) and RtHr (5’-GGY TTA CAA GAC CAG YGC TTT-3’; Doronin et al. 2013). PCR products were purified for sequencing by electrophoresis in 6% PAAG. Sequencing was performed using a ABI 3500 automatic sequencer (Applied Biosystems), using BigDye Terminator 3.1 kits (Applied Biosystems), and the same set of primers used for the PCR. After primer trimming, cyt b sequences used for analysis had a length of 1131 bp and were deposited in GenBank (Table S1).

The sequences were aligned and edited manually using BioEdit v. 7.0 software (Hall 1999). The genetic variability analysis was based on the calculation of the following parameters using DnaSP v.5.10.01 (Librado and Rozas 2009): total number of polymorphic positions (S), number of haplotypes (H), diversity of haplotypes (h), diversity of nucleotides per site (π), the average number of nucleotide substitutions (k). Neutrality tests Fu’s Fs (Fu 1997) and Tajima’s D (Tajima 1989) were performed in DnaSP v5.10.01. Before running the phylogenetic analyses, the dataset was tested for redundancy and saturation using DAMBE 6.4.101 software (Xia 2018) by calculating the entropy-based index of substitution saturation (Iss) and its critical value (Iss.c) (Xia et al. 2003; Xia and Lemey 2009). Each set of cyt b codon positions were tested separately. When the Iss values were significantly lower than the Iss.c (p < 0.001), indicates little to no substitution saturation and suitability of the data for phylogenetic inference. The IQ-TREE 2.0.5 software (Nguyen et al. 2015) was used for the maximum likelihood (ML) reconstruction of the phylogenetic relationships. The ModelFinder (Kalyaanamoorthy et al. 2017) macros in IQ-TREE 2.0.5 (Nguyen et al. 2015) was used to choose the best fit models of nucleotide substitution and optimal partitioning scheme (by codon position) using Bayesian Information Criterion (BIC). IQ-TREE was ran with the HKY+F+I+G4 substitution model for the 1st and 2nd codon positions combined and TIM2+F+G4 for the 3rd codon position and replicated using 1000 bootstraps. The Bayesian phylogenetic Inference (BI) was done using MrBayes 3.2.6 (Ronquist and Huelsenbeck 2003) with the optimal partitioning scheme selected in PartitionFinder 2.1.1 (Lanfear et al. 2012) by BIC: K80+I+G for the 1st codon position, HKY+I for the 2nd codon position, and GTR+G for the 3rd codon position. We conducted two simultaneous runs of four Markov Chains Monte Carlo (MCMC) tests, each run consisted of 4*106 generations with a sampling frequency every 500 generations. The first 25% of generations were discarded as burn-in. The convergence of runs was assessed by examination of the average standard deviation of split frequencies and the potential scale reduction factor. Stationarity was confirmed by examining posterior probability, log-likelihood, and all model parameters by the effective sample sizes (ESSs > 200) and trace plots of MCMC output in the program Tracer 1.7 (Rambaut et al. 2018). For outgroups, we used D. parvula, D. daghestanica, D. caucasica, and L. agilis (Table S1). The phylogenetic tree was visualized and edited using FigTree 1.4.4 (Rambaut 2018). The reconstruction of the haplotype network was carried out by the parsimony network algorithm of TCS (Clement et al. 2000), with 95% connection limit in PopART software (Leigh and Bryant 2015). Genetic distances (p-distances) were calculated in MEGA7 (Kumar et al. 2016).

Divergence time estimations

Various lizard groups were used as “external” calibration age constraints in our analysis (Tables S1, S2). We used the previously published six calibration points (Mya, millions of years ago ± standard deviation) widely applied in Lacertidae phylogenetic studies (e.g. Carranza and Arnold 2012; Tamar et al. 2016; Psonis et al. 2018): the split between Gallotia and Psammodromus algirus (emergence of Canary Islands): normal distribution, mean 18.0±2.0 (Cox et al. 2010; Carranza and Arnold 2012); the divergence between the genera Lacerta and Timon (Čerňanský 2010): 17.50±0.30; the split between G. galloti and G. caesaris (emergence of La Gomera Island; Cox et al. 2010; Carranza and Arnold 2012): 6.0±0.30; the separation of Crete from the Peloponnese leading to the divergence of Podarcis peloponnesiacus from P. cretensis and P. levendis (Meulenkamp 1985; Schule 1993; Poulakakis et al. 2005): 5.30±0.10; the divergence of P. lilfordi from P. pityusensis (Terrasa et al. 2004; Brown et al. 2008): 5.25±0.03 Mya; colonization of El Hierro Island by G. c. caesaris from the La Gomera Island (Cox et al. 2010; Guillou et al. 1996): 1.05±0.20 Mya. Only one sequence per species/lineage was retained. The divergence times were estimated using BEAST 1.10.4 (Drummond and Rambaut 2007). The TrN+I+G substitution model was selected in jModelTest2 (Darriba et al. 2012). The analysis was run for 107 generations with a sampling frequency 1000 generations, from which 10% were discarded as burn-in. A relaxed uncorrelated lognormal clock model, Yule process of speciation, and random starting tree were applied. All branches involved in the respective calibration points at the tree were forced to be monophyletic, as well as forced monophyly for P. levendis and P. cretensis, and D. alpina and D. lindholmi were allowed to follow the results of phylogenetic inference. The analysis was repeated three times and parameter log files and the phylogenetic trees were combined using LogCombiner 1.10.4. To assess the convergence and effective sample sizes (for ESSs > 200) for all parameters, we used Tracer 1.7.1 (Rambaut et al. 2018). The final phylogenetic tree was calculated in TreeAnnotator 1.10.4. The phylogenetic trees were visualized using FigTree 1.4.4 software (Rambaut 2018).

Modeling of the potential ranges and ecological niche analysis

The potential distribution estimates for the main evolutionary lineages of D. lindholmi were calculated using eight algorithms (Classification Tree Analysis, Multivariate Adaptive Regression Splines, General Additive Models, Generalized Linear Models, Artificial Neural Networks, MaxEnt, Generalized Boosted Models, and Random Forests) by performing 10 replicates using a bootstrap approach in the “SSDM” R package (Schmitt et al. 2017). By using multiple algorithms to estimate the distributions of each lineage, we ensured that multiple approaches were used to accurately estimate the distribution of species, and to reduce error in the models, even when using a low number of occurrences (Araujo et al. 2019). The distribution estimates were carried out independently for the three main genetic lineages using 44 localities for the Common lineage, nine for the Southwestern lineage, and 12 localities for the newly detected, deeply diverged Central lineage from the middle of the Crimean Mountains (Table S1). We also used data surrounding the pixels of known distribution to ensure the necessary fit for the models, a strategy used by Gherghel et al. (2020). The inclusion of each D. lindholmi population in any given phylogenetic lineage using this analysis was based on the results obtained from the molecular analysis (see Fig. 1 and Table S1). The distribution estimates generated by each algorithm were selected and stacked using an AUC (Area Under the Curve) cut-off of 0.75, where only the model estimates above this value were used to produce the final stacked model for each lineage distribution estimate. This method was chosen due to a low number of occurrences in the training dataset (comprising 75% of the occurrences) and to ensure that the estimates were robust, as it was previously recommended and used in similar applications of species distribution models (Schmitt et al. 2017; Papes and Gaubert 2007). To test the accuracy of the final stacked models, we used an independent test dataset comprising of 25% of occurrences to calculate the omission error, AUC, and Kappa metrics for each distribution estimate of the three lineages. Both the AUC and Kappa are metrices ranging from 0 (no model fit) to 1 (perfect model fit) (Hand and Till 2001) and are used to test species distribution models concerning the environments in which the species lives (Mouton et al. 2010, Jiménez-Valverde 2012). The omission error is used to test the goodness of fit for the models to the known geographic occurrence data. Models with low omission error are better concerning the known distribution of the species (Gherghel et al. 2020). Combining the three approaches to test the models is preferable to ensure the models were properly tested in both geographic and environmental spaces (following the assumptions of the Hutchinson duality paradigm of species distributions; Soberon and Peterson 2005).

To model the potential distribution of each detected lineage as well as the species, we used the WorldClim 2 database at 1 km resolution (; Fick et al. 2017). The WorldClim 2 database comprises 19 climatic variables that are a combination of temperature and precipitation factors previously shown to be highly important in shaping species distributions. The first step in selecting the predictors for our models was to ensure that the variables were not highly correlated and hence to eliminate the collinearity among them by calculating Pearsons’s correlation coefficients for all pairs of variables using SDMToolbox in ArcGIS 10 (Brown 2014). We excluded the variables of a correlated pair with │r│> 0.7 (see Table S3) which are considered to be the less biologically important. The resulting dataset contained seven bioclimatic variables (Table 3): bio3 (isothermality; %), bio7 (temperature annual range; °C × 10), bio8 (mean temperature of the wettest quarter of the year; °C × 10), bio9 (mean temperature of the driest quarter of the year; °C × 10), bio11 (mean temperature of the coldest quarter of the year; °C × 10), bio14 (precipitation of the driest month; mm), and bio 15 (precipitation seasonality; coefficient of variation). We also assessed the importance of the environmental predictors used in the analysis; and to show the relationship between the probability of the occurrence and the environmental variables, we performed Pearson correlation analyses (Schmitt et al. 2017) (Table S1). To predict the range dynamics of the species during the Last Glacial Maximum (LGM – 21 Kya ago) we have projected the present-day model for the species on three paleoclimatic global circulation models (CCSM4, MIROC, and MPI) based on CMIP5 (Coupled Model Intercomparison Project Phase 5). The potential distribution range during the LGM was then reported as the model ensemble of the three paleoclimatic global circulation models. The final model outputs consisted of suitabilities ranging from zero (no suitability) to 1 (maximum suitability).

To further test the unicity of the ecological niches occupied by the main D. lindholmi mtDNA lineages (see Results) we performed two niche overlap metrics to assess the degree of observed ecological overlap between each lineage by performing Schoener’s D (Warren et al. 2008) and standardized Hellinger’s-based I (Schoener 1968) niche overlaps in ENMTools v. 1.4.4 (Warren et al. 2010). Schoener’s D calculates the suitable range for a given evolutionary unit based on probability distributions for inhabiting particular regions (cells), calculating niche overlap based upon species abundance in those locations. Hellinger’s-based I is based purely on probability distributions without the assumptions of Schoener’s D (Warren et al. 2010). The value of each niche overlap index was statistically measured from zero, in the case of complete niche divergence/absence of the overlap, to 1 if niches are identical and they completely overlap (Rödder and Engler 2011). To further understand whether the similarity or dissimilarity of the ecological niches was due to random chance, we performed a background identity test where we built 100 pseudoreplicates, using the eight algorithms mentioned above, to produce a null model distribution of each lineage based on their respective range of distribution in Crimea to compare with the observed value of niche overlap of the two metrics used (Schoener’s D and Hellinger’s I; Warren et al. 2008). The respective range of distribution for each lineage was computed using a convex-hull-based approach in the ArcGIS Spatial Tools toolkit that was used to select the background from which the pseudopresences were randomly generated. The observed values of niche similarity between lineages were compared with the null distribution to test their statistical significance (Warren et al. 2008). The final distribution maps were created in ArcGIS 10.


Phylogenetic relationships and genetic variability

Among 101 sequences of D. lindholmi, three different mitochondrial lineages were identified: The Common, inhabiting most of the species range; the Southwestern, localized in the extreme south of the Crimean Peninsula, and the strongly differentiated Central lineage, confined to the middle part of the Crimean Mountains. The phylogenetic reconstructions, based on the ML and BI, had similar topologies.

When comparing D. lindholmi with other members of the D. saxicola complex (Table 1), the minimum values of the uncorrected p-distances were between the Common and the Southwestern lineages and D. szczerbaki (4.0±0.6% in both cases), while the highest value was between the Central lineage and D. brauneri (5.5±0.6%). Simultaneously, the latter value of p-distance was established as a maximum in the pairwise comparison among all members of D. saxicola complex. The minimum p-distance (2.7±0.5%) within this species complex was found between D. szczerbaki and D. saxicola.

Table 1.

The average genetic difference (p-distantce, %) between cyt b sequences per site (below the diagonal) and standard deviation (above the diagonal) within the Darevskia saxicola complex.

Species/lineage D. lindholmi/ Common D. lindholmi/ Southwestern D. lindholmi/ Central D. saxicola D. brauneri D. szczerbaki
D. lindholmi/ Common 0.3 0.6 0.6 0.6 0.6
D. lindholmi/ Southwestern 0.9 0.6 0.6 0.6 0.6
D. lindholmi/ Central 4.4 4.1 0.7 0.6 0.6
D. saxicola 4.9 4.7 5.2 0.7 0.5
D. brauneri 4.7 4.4 5.5 5.2 0.6
D. szczerbaki 4.0 4.0 4.7 2.7 4.9

The average value of uncorrected pairwise genetic distances in the analyzed dataset is 1.6±0.2%. The genetic distance of 4.4±0.6% was found between the Central and the Common lineages, 4.1±0.6% between the Central and Southwestern lineages, and 0.9±0.3% between the Southwestern and Common lineages (Table 1). Maximum and mean values of the intra-lineage divergence are 0.3% (0.1±0.0) for the Southwestern lineage, 0.4% (0.1±0.0) for the Central, and 1% (0.4±0.1) for the Common lineage.

Due to the length of the different sequences of the final dataset, we analyzed the haplotype networks based on 92 sequences of the species where we found 35 unique haplotypes. The Common lineage (67 samples from 42 localities) comprises most of the analyzed samples and represents the highest number of haplotypes (27 out of 35; Table 2). Some samples inside this lineage present further substructuring, namely those from the small areas at the western and eastern parts of the Main Range of the Crimean Mountains (respectively, D and F in Fig. 2, and Fig. 3). The remaining four subdivisions correspond to large physical-geographical units differing in their landscape and climatic characteristics: the extreme eastern and northern parts of the Crimean Mountains including coastal, low-mountain, and foothill regions (A), the eastern (B) and western (E) parts of the mountain-forest Crimea, and the northwestern forest-steppe foothills (C). The haplotype network supports the grouping into the lineages and subdivisions mentioned above (Fig. 3). The Southwestern (nine samples from six localities) and the Central lineages (16 samples from 12 localities) are not deeply structured: only four haplotypes were identified in each (Fig. 3). The lowest genetic diversity was found in the Central lineage (Table 2). Negative values of the Fs-test and Tajima’s D test in the Common lineage provide evidence for an increase in growth of the populations and past expansion events (Table 2).

Figure 2. 

Bayesian tree reconstructed with 1131 bp cytochrome b gene sequences. Numbers show posterior probabilities/ bootstrap support (values above 50% are given); **indicates 1.0 or 100% support, while *indicates 0.95–0.99 or 80–99% support. Haplogroup colors reflect the identified mtDNA lineages of Darevskia lindholmi. The numbers on the terminal branches of the tree correspond to the localities on Fig. 1. The capital letters (A–F) indicate the main subdivisions inside the Common lineage. The individuals pictured represent Darevskia lindholmi of both sexes: the Central lineage, Ayudag Mount (above), the Common lineage, Sevastopol city (photos by Oleg Kukushkin and Vitaly Giragosov, respectively).

Figure 3. 

The haplotype network for Darevskia lindholmi based on cytochrome b. Haplogroup colors correspond with the identified mtDNA lineages. The numbers of transverse strokes on the branches (in bold) correspond to the number of nucleotide substitutions. The numbers on network elements correspond to the localities on Fig. 1. The capital letters (A–F) represent substructure within the Common lineage (see Figs 1, 2). Expected ranges of D. lindholmi lineages and haplogroups of the Common lineage are shown on the map.

Table 2.

Genetic diversity in Darevskia lindholmi: n – sample size; H – number of haplotypes; h – haplotype diversity; π – nucleotide diversity (per site); S – total number of polymorphic positions; k – average number of nucleotide substitutions; Fu’s Fs is the Fs test value; Tajima D.

Lineage n H h±SD π±SD (%) S k Fu’s Fs Tajima D
Common 67 27 0.915±0.021 0.40±0.02 33 4.53 –11.61* –1.11
Southwestern 9 4 0.583±0.183 0.08±0.03 4 0.88 –1.28 –1.49
Central 16 4 0.442±0.145 0.07±0.03 5 0.83 –0.68 –1.61
D. lindholmi in the whole 92 35 0.935±0.013 1.60±0.19 81 18.14
*P support

Time of the divergence

Darevskia lindholmi divergence from other members of the genus is estimated at 2.1 Mya (95% HPD: 1.6–2.7 Mya). The data suggests that the Central lineage became isolated 1.8 Mya (95% HPD: 1.2–2.3 Mya). The split of the Southwestern and the Common lineages occurred only 0.5 Mya (95% HPD: 0.2–0.8 Mya). The mean (as well as the 95% HPD interval) estimated divergence times inside the D. saxicola complex and main lineages within D. lindholmi are presented in Fig. 4.

Figure 4. 

Maximum credibility chronogram (MCC) inferred in BEAST v. 1.10.4 with bars, showing the 95% HPDI of divergence time estimates. The vertical gray line marks the Messinian Salinity Crisis (MSC; 5.96–5.33 Mya). Red circles at the nodes indicate used calibration points, white circles indicate posterior probabilities above 0.95, and black circles denote groups of taxa forced to be monophyletic in the analysis.

Distribution modeling and niche differentiation

Our models performed very well in estimating the potential distributions and ecological niches of D. lindholmi, as demonstrated by both the AUC and Kappa metrics [Central: AUC= 0.96 (±0.04), Kappa= 0.65 (±0.24); Common: AUC= 0.97 (±0.02), Kappa= 0.7 (±0.26); Southwestern: AUC= 0.99 (±0.01), Kappa= 0.8 (±0.2)]. This suggests that the ecological niches estimated are highly reliable for all genetic lineages modeled. The potential distribution of the lineages is also well estimated as shown by the low omission scores [Central lineage omission error = 0.05 (±0.06), Common omission error = 0.01 (±0.01), Southwestern omission error = 0.01 (±0.01)]. The intermodel variability also shows a high degree of stability in our ensamble among all eight algorithms used (Fig. S1), which also suggests that the models are highly reliable. The predictors with the highest power for predicting the ranges and ecological niches for all of the main evolutionary lineages of the Crimean Darevskia are the mean temperature of the wettest quarter of the year (bio8), temperature annual range (bio7), and the amount of precipitation during the driest month (bio14). The response curves of the three lineages to these three environmental predictors are very strong as well, and cumulatively are capable of explaining 57.2% of the variation from the potential distribution estimates across all lineages. To disentangle the direction of the relationship between the probability of distributions of each lineage and environemntal variables, we found positive correlations between the probability of presence of the lineages and isothermality (Common lineage vs. bio3: r=0.274; Central lineage vs. bio3: r=0.406; Southwestern lineage vs. bio3: r=0.209), the mean temperature of driest quarter of the year (Common lineage vs. bio9: r=0.232; Central lineage vs. bio9: r=0.29; Southwestern lineage vs. bio9: r=0.31), mean temperature of coldest quarter of the year (Common lineage vs. bio11: r=0.588; Central lineage vs. bio11: r=0.582; Southwestern lineage vs. bio11: r=0.625) and the amount of precipitation during the driest month (Common lineage vs. bio14: r=0.468; Central lineage vs. bio14: r=0.719; Southwestern lineage vs. bio14: r=0.537), while we found negative correlations between the temperature annual range (Common lineage vs. bio7: r= –0.743; Central lineage vs. bio7: r= –0.82; Southwestern lineage vs. bio7: r= –0.77), the mean temperature of wettest quarter of the year (Common lineage vs. bio8: r= –0.847; Central lineage vs. bio8: r= –0.913; Southwestern lineage vs. bio8: r= –0.869) and the precipitation’s seasonality (Common lineage vs. bio15: r= –0.196; Central lineage vs. bio15: r= –0.128; Southwestern lineage vs. bio15: r= –0.018). The potential ranges based on the fundamental niches of the Common, Southwestern, and Central lineages suggests the possibility of some potential geographic overlap in areas where the lineages enter in contact, however, this was not observed based on the omission error of the models which suggests that likely the lineages do not cohabit the same localities (Fig. 1, Fig. 5). The optimal fundamental niche for the Central lineage is significantly smaller and limited primarily to the central region of the Crimean Mountains (Fig. 5) as also observed in the occurrence data (Fig. 1). During the Last Glacial Maximum, our models suggest that the species had its refugia exclusively in the Crimean Mountains (Fig. 6).

Figure 5. 

Potential distribution of Darevskia lindholmi evolutionary lineages: Central, Common, and the Southwestern. The intensity of the colour is related with the degree of suitability for the species presence. The black polygon represents the minimum convex polygon presence localities for each lineage.

Figure 6. 

Potential distribution of Darevskia lindholmi in present-days and during the Last Glacial Maximum (21 Kya). The intensity of the colour is related with the degree of suitability for the species presence.

To test the unicity of the niches modeled in relation to the available environmental conditions we performed a background analysis where we found significant differences between the niches occupied by the three lineages. The observed niche overlap between the lineages was very high (Central vs. Common: Schoener’s D=0.82, Hellinger’s-based I=0.97; Central vs. Southwestern: Schoener’s D=0.80, Hellinger’s-based I=0.97; Common vs. Southwestern: Schoener’s D=0.90, Hellinger’s-based I=0.98). However, due to small gradients, these overlaps were still shown to represent different niches based on the identity background analysis (t-test, df =99, P<0.05). The identity test indicated that the null-hypothesis regarding niche overlap can be rejected and that Central lineage and Common lineage (DH0 is 0.85±0.02 vs. DH1 = 0.82, and IH0 is 0.98±0.007 vs. IH1 = 0.97) occupy different niches, and the Central lineage and Southwestern lineage (DH0 is 0.90±0.06 vs. DH1 = 0.80, and IH0 is 0.98±0.02 vs. IH1 = 0.97) are marginally, but significantly, differentiated in their respective niches. The identity background test also shows that the Common lineage and the Southwestern lineage do not occupy different niches and a high degree of overlap exists (DH0 is 0.88±0.06 vs. DH1 = 0.90, and IH0 is 0.88±0.07 vs. IH1 = 0.98), which is in congruence with the observed genetic relationships described above (Fig. 3). As such, we can conclude that the ecological niches occupied by the considered groups are not identical, and the Central and Common evolutionary lineages occupy distinct ecological niches.



To reconstruct the origin of D. lindholmi in Crimea is challenging due to the controversial paleogeography of the North Black Sea region as well as the existence of some gaps on the local orogeny, especially from later geological periods (Shalimov 1966; Shnyukov et al. 1997; Esin et al. 2018; Popov et al. 2010, 2019; Stovba et al. 2020; Palcu et al. 2021), and limited fossil evidence (Vremir and Ridush 2005; Ratnikov 2015; Kovalchuk et al. 2020). Current paleogeographic reconstruction suggest that during the Plio-Pleistocene transition (the Akchagylian epoch, 3.6–1.8 Mya), the Black Sea and the Caspian-Azov marine basins were separated by land located approximately in the Kerch – Taman region or even south-lying areas (Lavrischev et al. 2011; Krijgsman et al. 2019). The existence of this hypothetical corridor might thus be the reason for the high similarity of biota between Crimea and the northwestern Caucasus (Novosad 1992). Based on this, we assume that the isolation of D. lindholmi might be related to the time of disappearance of this land-bridge, which corresponds with our estimated time of divergence of D. lindholmi from other species of the D. saxicola complex (1.6–2.7 Mya; Fig. 4).

Our data also shows intraspecific genetic divergence in D. lindholmi, although it is less diverged in comparison with some Caucasian populations of the D. saxicola complex (Doronin et al. 2013; Tarkhnishvili et al. 2016). The Central lineage is sister to the Common and the Southwestern lineages with uncorrected p-distance 4.6±0.6%, suggesting a long-term independent evolution related to topographical and climatic changes in the region, probably during the Early Pleistocene (see Babak 1959; Pisareva et al. 2019; Sirenko 2019).

Present-day deep-water parts of the Black Sea were formed at the end of the Pliocene and beginning of the Pleistocene (Nikishin et al. 2003). It means that during the Akchagylian epoch the powerful, rapid tectonic subsidence of the Black Sea hollow occurred, which was accompanied by the sinking of peripheral parts of surrounding mountain systems. The subsidence of the southern “wing” of the Crimean orogen below the sea level led to the emergence of isolated plots of land, like off-shore islets and peninsulas with complex coastlines (Shalimov, 1966; Stovba et al. 2020). The influence of sea-level fluctuations as well as the creation of land-bridges between different land masses is known to have determined the genetic structuring in other lizard species (Podnar et al. 2004; Salvi et al. 2014; Senczuk et al. 2017; Bernardo et al. 2019). Anyway, it can be assumed that the Central lineage was separated about 2 Mya or even later (1.2–2.3 Mya; Fig. 4) apparently during temporary disintegration of the Crimean Uplift into several isolated areas or due to other landscape changes.

The level (less than 1%) and time (0.2–0.8 Mya; Fig. 4) of difference between the Southwestern and the Common lineages suggest that their divergence could be later associated with substantial cooling during the Middle-Late Pleistocene (Sirenko 2019). The closest ancestral haplotypes of D. lindholmi represented by the Southwestern lineage probably persisted in the extreme south of Crimea (potential refugium), where the influence of climatic oscillations was less pronounced (Gerasimenko 2007, 2011; Pisareva et al. 2019). A similar pattern is expected for the Crimean montane subspecies of L. agilis (Kukushkin et al. 2020). The subsequent population divergence of D. lindholmi could thus correspond with landscape and climatic changes or niche differentiation (see below). In general, the aforementioned events support our data of the divergence in the detected main lineages of D. lindholmi (Fig. 4). Alternatively, the Common and the Southwestern lineages might simply represent geographically structured diversity within a single population that was never split into allopatric areas of vicariance but was subjected to spatially restricted gene flow due to isolation by distance and habitat fragmentation.

The range dynamics of D. lindholmi distribution pattern could be dependent on the dynamics of the forest massif boundaries. Our models on the range dynamics of D. lindholmi based on climatic niche modeling suggest that during the Last Glacial Maximum (21 Kya) the species inhabited the Crimean Mountains (Fig. 6). As D. lindholmi prefers rocky microhabitats, using niche conservatism theory we expect that the species would have inhabited similar habitats during the Last Glacial Maximum in the warmest and wettest areas, such as large rock massifs and warm river canyons of the Crimean Mountains. The most favorable conditions during this period developed probably for the Central lineage (Fig. S4). When the climate got warmer and wetter in the Holocene (Cordova et al. 2011), lizards recolonized the previously unsuitable regions of the Crimean Mountains. The genetic uniformity of the Central lineage allows us to assume that it persisted within the geographically restricted area of the Crimean Mountains or alternatively it lost their genetic diversity due to extinciton events of local populations. The observed structure of the Common lineage (Fig. 3) might be the result of a rapid post-glacial recolonization of the mountain areas from several refugia located on the southern macroscope of the Crimean Mountains. The star-like pattern in the haplotype network, as well as neutrality tests, indicate past dispersal events of the Common lineage, possibly related to sea-level changes during the Mid-Holocene, later 8.5 Kya (see Lericolais et al. 2010).

Distribution and ecology

Distribution ranges of the lizard genetic lineages do not show any significant overlap (Fig. 1; although it could be due to sampling bias). Based on the current knowledge, the phylogeographical pattern of D. lindholmi lineages looks exceptional since the range of the Central lineage splits the area inhabited by the Common lineage. The range of the Central lineage is shaped like a wide strip oriented along the north-south axis and lies eastward from Yalta and the Simferopol meridian and westward to the vast Karabi-Yaila Plateau (Fig. 1). Thus, the Central lineage distribution is confined to the highest and wettest part of the Main Range of the Crimean Mountains as well as the most humid forested part of the Southern Coast. In the foothills, its range is limited to the valleys of the small rivers Zuya and Burulcha (right tributaries of the Salgir) where only two isolated populations were recorded, virtually surrounded by a flat landscape. The range of the Southwestern lineage is restricted to a narrow strip of the precipitous shore below 750 m a. s. l. (i.e. area with the most pronounced sub-Mediterranean character), whereas the Common lineage occupies the rest of the species range (Fig. 1). The altitudinal diapason of the most diverged evolutionary lineages’ habitats overlaps significantly. Nevertheless, the average elevation of their representative’s observations has an almost two-fold difference: the Common with the Southwestern, 1–1200 m a. s. l. (350±35.9 m), the Central lineage, 50–1320 m (666±72.2 m).

Based on observations made in the field, D. lindholmi can be characterized as a mesophilic cold-tolerant species. Its distribution at low altitudes is limited by insufficient humidity, while in the highlands by the average temperature of the summer months (Kukushkin 2009). For these reasons, this lizard is absent in arid rocky areas in the eastern extremity of the Main Range, as well as at the summits above 1300–1500 m a. s. l. in different parts of highlands, where the average temperature of the hottest months is below 15 °C. Analysis of climatic maps (Ved’ 2000) revealed that the Central lineage is confined to humid areas with annual precipitation of about 500–1100 mm. The Common and Southwestern lineages are more resistant to the arid climate and inhabit coastal areas with annual precipitation near 320–350 mm, as well as cold and wet highlands where precipitation reaches 600–1000 mm. Overall, the Common and Central lineages are documented in all physical-geographical areas of the Crimean Mountains and dwell in a wide range of habitats with diverse climatic characteristics and vegetation (Fig. S3). However, our field observations suggest that the moderately humid, warm climate of forested terrains at an elevation below 1000 m a. s. l. is the most favorable for representatives of both lineages (O. Kukushkin personal data). These observations are congruent with our results obtained from the ecological niche modeling (Fig. 5) which confirms the assumptions based on in situ observations.

Our models show significant variation in the response of each lineage to the environmental variables used to model their distributions. In general, three variables accounted for almost two-thirds of the variation (Table 3) such as the mean temperature of the wettest quarter of the year (bio8), temperature annual range (bio7), and the amount of precipitation during the driest month (bio14). Previous work on the potential distribution of D. lindholmi by Doronin (2012b) found similar environmental predictors to be of the highest importance in shaping species distributions as our models (Table 3). At a closer look, we found that each genetic lineage favors a suite of climatic conditions that have little overlap (Fig. S2); from which we can conclude that these lineages have diverged in climatic tolerances, especially regarding their ability to survive droughts and prolonged periods of particularly high and low temperatures. The ecological niches occupied by the considered groups are not identical, although show a high degree of overlap. The divergence of ecological niches between the studied populations might occur because of their different attitudes towards the humidification factor (Fig. S2), though there is a wide diapason of habitats that have intermediate suitability for all groups. Ecological niches for the closely related species are expected to be more similar compared with deeply diverged species, generally following a pattern of niche conservatism (Wiens and Graham 2005). Evidence of both niche conservatism and niche divergence has been reported in other representatives of the Darevskia genus (Ahmadzadeh et al. 2013; Kurnaz et al. 2019; Kurnaz and Yousefkhani 2020; Rato et al. 2020). As shown in previous work, closely related species tend to occupy similar ecological niches, while divergence can occur when significant shifts in habitat use are followed by morphological changes of functional traits (Schulte et al. 2004; Collar et al. 2010). In the case of D. lindholmi, we detect a pattern of niche conservatism, where niches of the lineages follow a pattern of niche overlap that is mirroring the lineage divergence, with lineages that are less divergent having more similar niches and lineages that are more divergent having dissimilar niches. However, some examples of the significant climatic niche overlap in deeply divergent sister evolutionary lineages of Lacertidae have been previously documented, which is considered to be a result of ancient allopatric speciation under similar environmental conditions (Caeiro-Dias et al. 2018). Thus, the current isolation of the Central lineage along the border of nearly its entire range is substantially determined by the geomorphological features, primarily by the presence of bedrock outcrops to which the populations of this petrophilous lizard are usually confined. As such, based on our field observations, the Central lineage distribution show no geographic overlap in with the rest of the lineages. Intriguingly, the Central group splits the Common group distribution in half, and our analyses show limited to no level of mitochondrial admixture between the two lineages. This situation resembles the Triturus cristatus case that splits the range of the closely related T. dobrogicus range in two with relatively strong levels of admixture in areas of contact between the two species (Wielstra et al. 2013). This is the result of past range shifts where the species track their fundamental niches with climate change throughout the Quaternary oscillations by using corridors that might be hidden or undetectable in the contemporary species distributions (Gherghel and Papeş 2015, Vörös et al. 2016). Similarly, in the case of the D. lindholmi lineages, we suspect a similar process of tracking suitable ecological niches of the lineages has occurred throughout the Quaternary glacial oscillations which resulted in the observed contemporary patterns of genetic structuring in the species (Figs 1 and 3).

Table 3.

The percentage contribution (%) of selected bioclimatic predictors to estimating the potential distributions and ecological niches of Darevskia lindholmi lineages.

Bioclimatic variables Code Central Common Southwestern
Isothermality bio3 8.87 10.25 12.27
Temperature annual range bio7 23.65 27.77 16.45
Mean temperature of wettest quarter bio8 16.26 18.67 19.67
Mean temperature of driest quarter bio9 9.33 10.07 12.52
Mean temperature of coldest quarter bio11 11.55 7.20 10.86
Precipitation of driest month bio14 21.14 16.06 13.17
Precipitation seasonality bio15 9.21 9.98 15.06


The uncorrected p-distance and time of divergence of the Central lineage (Tab. 1, Fig. 4) are comparable with the data of some taxa and/or mitochondrial lineages within lacertid lizards, e.g. Podarcis (Podnar et al. 2004; Psonis et al. 2017; Senczuk et al. 2017; Spilani et al. 2019) or Lacerta (Andres et al. 2014; Marzahn et al. 2016). In Darevskia, differences between 3.8–4.8% in cyt b correspond to species level detected for example in D. bithynica (Méhely, 1909). On the other hand, the taxonomic position of this species belonging to the D. rudis complex is supported by differences in biparental markers and ecological niches (Rato et al. 2020). Our data, however, allows the hypothesis that the Central lineage may represent a distinct taxon, most likely of the subspecies rank. In addition to mtDNA data, the observed parapatry of the lineages, differences between the niches, and the existence of certain morphological differences between them (O. Kukushkin and A. Svinin personal data; Figs S5 and S6), together with well-known agreement between mtDNA data and taxonomy in the western Palearctic reptiles (Joger et al. 2007), convinces us for such a conclusion. However, an integrative approach using ecological, morphological, and biparental genetic data is needed to resolve the taxonomy of the species (see Kindler and Fritz 2018; Abreu et al. 2020; Jauss et al. 2021). Future work on clarifying the taxonomic status of D. lindholmi is particularly needed, especially to better understanding the genetic lineages’ spatial isolation from each other, and to understanding the level of the lineages’ admixture. Moreover, ecological niches and the habitats used by the individuals corresponding to the D. lindholmi lineages are quite different (Fig. S3) and future studies should also focus on functional traits adaptations to these habitats as well, as they could provide important insight into the initial steps of speciation in the Darevskia group. The estimated time of the independent evolutionary history inside D. lindholmi implies that ranges between lineages could be in contact, which predicts potential gene flow (see hybridization between closely related Lacertidae species in areas of sympatry; Tarkhnishvili et al. 2013; Caeiro-Dias et al. 2021). Moreover, further investigation of populations close to the type locality (Yalta city) is needed. At present, the haplogroup of the Central lineage was found on the seaside slope of the Main Range northwards and eastwards from the city of Yalta (Fig. 1). In turn, the haplogroup of the Common lineage was discovered not only westwards from Yalta, but also within city limits (O. Kukushkin and O. Ermakov personal data, 2021). Particularly, future work to solve the taxonomic situation within D. lindholmi is recommended, particularly to elevate the Central lineage to a distinct taxon in the light of the current results presented here. Ultimately, a taxonomic revision of the whole D. saxicola complex remains a challenging, and yet warranted task.


We would like to express sincere gratitude to the following persons for the assistance during the field works, the providing of photos, maps, and other materials (alphabetically): Vitaly Giragosov, Roman Gorelov, Marina Khramova, Marina Khrisanova, Yuliya Krasylenko, Spartak Litvinchuk, Konstantin Milto, Anton Nadolnyi, Pavel Oksinenko, Ekaterina Polyakova, Elena Sviridenko, Oksana Tishchenko, Sergey Tokarev, Aleksandr Trofimov, Ilya Turbanov, and Larisa Znamenskaya. We also thank to Riley Tedrow, for valuable feedback and improving the English of the manuscript at its final stages. We thanks the anonymous reviewers for their attention to our work and criticism, which contributed to the improvement of the manuscript quality. The work of O.K. and I.D. was carried out within the framework of the State assignments Nos 121032300023-7 and АААА-А19-119020590095-9, respectively. O.K. and O.E. conducted research at personal expenses without involving additional sources of funding. The participation of I.D. was funded from the research grant No. 19-04-00514 of the Russian Foundation of Basic Research. The work of D.J. was supported by the Slovak Research and Development Agency under contract No. APVV-19-0076.


  • Abreu JM, Salvi D, Perera A, Harris JD (2020) Mitochondrial lineages or discrete species? Assessing diversity within Timon tangitanus (Lacertidae) using mitochondrial and nuclear DNA sequences. Amphibia-Reptilia 41: 123–132.
  • Ahmadzadeh F, Flecks M, Carretero MA, Mozaffari O, Böhme W, Harris DJ, Freitas S, Rodder D (2013) Cryptic speciation patterns in Iranian rock lizards uncovered by integrative taxonomy. PLoS ONE 8(12): e80563.
  • Aljanabi SM, Martinez I (1997) Universal and rapid salt-extraction of high genomic DNA for PCR-based techniques. Nucleic Acids Research 25(22): 4692–4693.
  • Araújo MB, Anderson RP, Barbosa AM, Beale CM, Dormann CF, Early R, Garcia RA, Guisan A, Maiorano L, Naimi B, O’Hara RB (2019) Standards for distribution models in biodiversity assessments. Science Advances: 5(1): eaat4858.
  • Arnold EN, Arribas O, Carranza S (2007) Systematics of the Palaearctic and Oriental lizard tribe Lacertini (Squamata: Lacertidae: Lacertinae), with descriptions of eight new genera. Zootaxa 1430: 1–86.
  • Babak VI (1959) A brief account of the neotectonics of the Crimea. Byulleten Moskovskogo obshchestva ispytateleyi prirody 34(4): 51–65. (In Russian).
  • Bernardo PH, Sánchez-Ramírez S, Sánchez-Pacheco SJ, Alvarez-Castaneda ST, Aguilera-Miller EF, Mendez-de la Cruz FR, Murphy RW (2019) Extreme mito-nuclear discordance in a peninsular lizard: the role of drift, selection, and climate. Heredity 123: 359–370.
  • Brown JL (2014) SDMtoolbox: a python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. Methods in Ecology and Evolution 5(7): 694–700.
  • Brown RP, Terrasa B, Pιrez-Mellado V, Castro JA, Hoskisson PA, Picornell A, Ramon MM (2008) Bayesian estimation of post-Messinian divergence times in Balearic Island lizards. Molecular Phylogenetics and Evolution 48: 350–358.
  • Caeiro-Dias G, Brelsford A, Kaliontzopoulou A, Meneses-Ribeiro M, Crochet P-A, Pinho C (2021) Variable levels of introgression between the endangered Podarcis carbonelli and highly divergent congeneric species. Heredity 126: 463–476.
  • Caeiro-Dias G, Luís C, Pinho C, Crochet P-A, Sillero N, Kaliontzopoulou A (2018) Lack of congruence of genetic and niche divergence in Podarcis hispanicus complex. Journal of Zoological Systematics and Evolutionary Research 56(4): 1–14.
  • Carranza S, Arnold EN (2012) A review of the geckos of the genus Hemidactylus (Squamata: Gekkonidae) from Oman based on morphology, mitochondrial and nuclear data, with descriptions of eight new species. Zootaxa 3378: 1–95.
  • Cordova CE, Gerasimenko NP, Lehman PH, Kliukin AA (2011) Late Pleistocene and Holocene paleoenvironments of Crimea: Pollen, soils, geomorphology, and geoarchaeology. Geological Society of America Special Papers 473: 133–164.
  • Darevsky IS (1967) Rock lizards of the Caucasus (systematics, ecology and phylogenesis of the polymorphic groups of Caucasian rock lizards of the subgenus Archaeolacerta). Nauka, Leningrad, 214 pp. (In Russian).
  • Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: More models, new heuristics and parallel computing. Nature Methods 9: 772.
  • Doronin IV (2012a) Review of type specimens of rock lizards of Darevskia (saxicola) complex (Sauria: Lacertidae). Proceedings of the Zoological Institute RAS 316(1): 22–39. (In Russian).
  • Doronin IV (2012b) The use of GIS for the analysis of the distribution of rock lizards Darevskia (saxicola) complex (Sauria: Lacertidae). Current Studies in Herpetology 12(3/4): 91–122. (In Russian).
  • Doronin IV, Dzhelali PA, Lotiev KY, Mazanaeva LF, Mustafaeva GA, Bunyatova SN (2021) Phylogeography of a Darevskia (caucasica) complex (Lacertidae: Sauria) based on the cytochrome b mitochondrial gene analysis. Proceedings of the Zoological Institute RAS 325(1): 49–66. (In Russian).
  • Doronin IV, Tuniyev BS, Kukushkin OV (2013) Differentiation and taxonomy of the rock lizards Darevskia (saxicola) complex (Sauria: Lacertidae) according to morphological and molecular analyses. Proceedings of the Zoological Institute RAS 317(1): 54–84. (In Russian).
  • Dufresnes C, Litvinchuk SN, Leuenberger J, Ghali K, Zinenko O, Stöck M, Perrin N (2016) Evolutionary melting pots: a biodiversity hotspot shaped by ring diversifications around the Black Sea in the Eastern tree frog (Hyla orientalis). Molecular Ecology 25(17): 4285–4300.
  • Esin NV, Yanko-Hombach VV, Esin NI (2018) Evolutionary mechanisms of the Paratethys Sea and its separation into the Black Sea and Caspian Sea. Quaternary International 465: 46–53.
  • Fick SE, Hijmans RJ (2017) WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37(12): 4302–4315.
  • Fritz U, Ayaz D, Hundsdörfer AK, Kotenko T, Guicking D, Wink M, Tok CV, Çiçek K, Buschbom J (2009) Mitochondrial diversity of European pond turtles (Emys orbicularis) in Anatolia and the Ponto-Caspian Region: multiple old refuges, hotspot of extant diversification and critically endangered endemics. Organisms, Diversity and Evolution 9: 100–114.
  • Fu Y (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147: 915–925.
  • Gerasimenko N (2007) Environmental changes in the Crimean mountains during the Last Interglacial – middle pleniglacial as recorded by pollen and lithopedology. Quaternary International 164–165: 207–220.
  • Gerasimenko NP (2011) Climatic and environmental oscillations in southeastern Ukraine from 30 to 10 ka, inferred from pollen and lithopedology. Geological Society of America Special Papers 473: 117–132.
  • Gherghel I, Brischoux F, Nyári AS, Papeş M (2020) A simple framework for estimating potential distributions of amphibious marine species and implications for conservation. Coral Reefs 39(4): 1081–1090.
  • Grechko VV, Ciobanu DG, Darevsky IS, Kosushkin SA, Kramerov DA (2006) Molecular evolution of satellite DNA repeats and speciation of lizards of the genus Darevskia (Sauria: Lacertiae). Genome 49(10): 1297–1307.
  • Guillou H, Carracedo JC, Torrado FP, Badiola ER (1996) K-Ar ages and magnetic stratigraphy of a hotspot-induced, fast grown oceanic island: El Hierro, Canary Islands. Journal of Volcanology and Geothermal Research 73: 141–155.
  • Hall TA (1999) BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symposium Series 41: 95–98.
  • Hand DJ, Till RJ (2001) A simple generalisation of the area under the ROC curve for multiple class classification problems. Machine Learning: 45(2): 171–186.
  • Jablonski D, Nagy ZT, Avci A, Olgun K, Kukushkin OV, Safaei-Mahroo B, Jandzik D (2019) Cryptic diversity in the smooth snake (Coronella austriaca). Amphibia-Reptilia 40: 179–192.
  • Jauss R-T, Solf N, Kolora SRR, Schaffer S, Wolf R, Henle K, Fritz U, Schlegel M (2021) Mitogenome evolution in the Lacerta viridis complex (Lacertidae, Squamata) reveals phylogeny of diverging clades. Systematics and Biodiversity: in press.
  • Jiménez-Valverde A (2012) Insights into the area under the receiver operating characteristic curve (AUC) as a discrimination measure in species distribution modelling. Global Ecology and Biogeography 21: 498–507.
  • Joger U, Fritz U, Guicking D, Kalyabina-Hauf S, Nagy ZT, Wink M (2007) Phylogeography of western Palaearctic reptiles – Spatial and temporal speciation patterns. Zoologischer Anzeiger 246: 293–313.
  • Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS (2017) ModelFinder: Fast model selection for accurate phylogenetic estimates. Nature Methods 14: 587–589.
  • Kehlmaier C, Zinenko O, Fritz U (2020) The enigmatic Crimean green lizard (Lacerta viridis magnifica) is extinct but not valid: Mitogenomics of a 120-year-old museum specimen reveals historical introduction. Journal of Zoological Systematics and Evolutionary Research 58: 303–307.
  • Kindler C, Fritz U (2018) Phylogeography and taxonomy of the barred grass snake (Natrix helvetica), with a discussion of the subspecies category in zoology. Vertebrate Zoology 68(3): 269–281.
  • Kovalchuk O, Rekovets L, Tsvelykh A, Yanenko V, Manko V, Tajkova S (2020) Living in a time of change: Late Pleistocene/Holocene transitional vertebrate fauna of Grot Skeliastyi (Crimea, Ukraine). Historical Biology.
  • Krijgsman W, Tesakov A, Yanina T, Lazarev S, Danukalova G, Van Baak CGC, Agustí J, Alçiçek MC, Aliyeva E, Bista D, Bruch A, Büyükmeriç Y, Bukhsianidze M, Flecker R, Frolov P, Hoyle TM, Jorissen EL, Kirscher U, Koriche SA, Kroonenberg SB, Lordkipanidze D, Oms O, Rausch L, Singarayer J, Stoica M, van de Velde S, Titov VV, Wesselingh FP (2019) Quaternary time scales for the Pontocaspian domain: Interbasinal connectivity and faunal evolution. Earth-Science Reviews 188: 1–40.
  • Kukushkin OV (2009) On patterns of spatial distribution of Lindholm’s rock lizard Darevskia lindholmi (Sauria: Lacertidae) in the South-Eastern Coast of the Crimea. Byulleten Samarskaya Luka: Problemy regionalnoyi i globalnoyi ekologii 18(1): 68–75. (In Russian).
  • Kukushkin OV, Ermakov OA, Ivanov AY, Doronin IV, Sviridenko EY, Simonov EP, Gorelov RA, Khramova MA, Blokhin IG (2020) Cytochrome b mitochondrial gene analysis-based phylogeography of a Sand lizard in the Crimea: ancient refugium at the peninsula, late expansion from the North, and first evidence of Lacerta agilis tauridica and L. a. exigua (Lacertidae: Sauria) hybridization. Proceedings of the Zoological Institute RAS 324(1): 56–99. (In Russian).
  • Kukushkin OV, Ivanov AY, Ermakov OA (2018) Genetic geterogenity of Marsh frog (Pelophylax (ridibundus) complex; Anura, Ranidae) population in the Crimea revealed by mitochondrial amd nuclear DNA analyses. University proceedings. Volga region. Natural sciences 3(23): 33–55. (In Russian).
  • Kumar S, Stecher G, Tamura K (2016) MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Molecular Biology and Evolution 33(7): 1870–1874.
  • Kurnaz M, Kutrup B, Yousefkhani SSH, Koҫ H, Bülbül U, Eroğlu AI (2019) Phylogeography of the red-bellied lizard, Darevskia parvula in Turkey. Mitochondrial DNA Part A 30(3): 556–566.
  • Lanfear R, Calcott B, Ho SYW, Guindon S (2012) PartitionFinder: Combined Selection of Partitioning Schemes and Substitution Models for Phylogenetic Analyses. Molecular Biology and Evolution 29(6): 1695–1701.
  • Lavrischev VA, Sheikov AA, Andreev VM, Semenov VM, Ermakov VA, Grekov II, Shelting SK, Shishov VS, Navodnichenko SP (2011) State Geological Map of the Russian Federation. Scale 1: 1 000 000 (third generation). Series Scythian. VSEGEI, Saint Petersburg, 431 pp. (In Russian).
  • Lericolais G, Guichard F, Morigi C, Minereau A, Popescu I, Radan S (2010) A post Younger Dryas Black Sea regression identified from sequence stratigraphy correlated to core analysis and dating. Quaternary International 225: 199–209.
  • Levin BA, Gandlin AA, Simonov ES, Levina MA, Anna E., Barmintseva AE, Japoshvili B, Mugue NS, Mumladze L, Mustafayev NJ, Pashkov AN, Roubenyan HR, Shapovalov MI, Doadrio I (2019) Phylogeny, phylogeography and hybridization of Caucasian barbels of the genus Barbus (Actinopterygii, Cyprinidae). Molecular Phylogenetics and Evolution 135: 31–44.
  • MacCulloch RD, Fu J, Darevsky IS, Murphy RW (2000) Genetic evidence for species status of some Caucasian rock lizards in the Darevskia saxicola group. Amphibia-Reptilia 21: 169–176.
  • Marzahn E, Mayer W, Joger U, Ilgaz Ç, Jablonski D, Kindler C, Kumlutaş U, Nistri A, Schneeweiss N, Vamberger M, Žagar A, Fritz U (2016) Phylogeography of the Lacerta viridis complex: mitochondrial and nuclear markers provide taxonomic insights. Journal of Zoological Systematics and Evolutionary Research 54(2): 85–105.
  • Mendel J, Lusk S, Vasil’eva ED, Vasil’ev VP, Lusková V, Ekmekci FG, Erk’akan F, Ruchin A, Koščo J, Vetešník L, Halačka K, Šanda R, Pashkov AN, Reshetnikov SI (2008) Molecular phylogeny of the genus Gobio Cuvier, 1816 (Teleostei: Cyprinidae) and its contribution to taxonomy. Molecular Phylogenetics and Evolution 47: 1061–1075.
  • Meulenkamp JE (1985) Aspects of the Late Cenozoic evolution of the Aegean Region. In: Stanley DJ, Wezel FC (Eds.) Geological Evolution of the Mediterranean Basin. Springer, New York, 307–321.
  • Muratov MV (1960) A brief outline of the geological structure of the Crimean Peninsula. Gosgeoltechizdat, Moscow, 207 pp. (In Russian).
  • Murphy RW, Fu J, MacCulloch RD, Darevsky IS, Kupriyanova LA (2000) A fine line between sex and unisexuality: the phylogenetic constraints on parthenogenesis in lacertid lizards. Zoological Journal of the Linnean Society 130: 527–549.
  • Murtskhvaladze M, Tarkhnishvili D, Anderson CL, Kotorashvili A (2020) Phylogeny of Caucasian rock lizards (Darevskia) and other true lizards based on mitogenome analysis: Optimisation of the algorithms and gene selection. PLoS ONE 15(6): e0233680.
  • Nguyen L-T, Schmidt HA, von Haeseler A, Minh BQ (2015) IQ-TREE: A fast and effective stochastic algorithm for estimating maximum likelihood phylogenies. Molecular Biology and Evolution 32: 268–274.
  • Nikishin AM, Korotaev MV, Ershov AV, Brunet M-F (2003) The Black Sea Basin: tectonic history and Neogene-Quaternary rapid subsidence modelling. Sedimentary Geology 156: 149–168.
  • Novosad VV (1992) Flora of Kerch-Taman region (comparative analysis of the structure, ecoflorotopological differentiation, genesis, prospects of the rational use and conservation). Naukova dumka, Kiev, 277 pp. (In Russian).
  • Papeş M, Gaubert P (2007) Modelling ecological niches from low numbers of occurrences: assessment of the conservation status of poorly known viverrids (Mammalia, Carnivora) across two continents. Diversity and Distributions 13(6): 890–902.
  • Pisareva VV, Faustova MA, Zyuganova IS, Karpukhina NV, Zakharov AL, Konstantinov EA, Semenov VV, Kurbanov RN (2019) Landscape and climate changes in Eastern Europe in the Early Pleistocene. Stratigrafia, Geologicheskaya korrelyatsiya 27(4): 93–116. (In Russian).
  • Podnar M, Mayer W, Tvrtković N (2004) Mitochondrial phylogeography of the Dalmatian wall lizard, Podarcis melisellensis (Lacertidae). Organisms, Diversity and Evolution 4: 307–317.
  • Popov SV, Antipov MP, Zastrozhnov AS, Kurina EE, Pinchuk TN (2010) Sea-level fluctuations on the northern shelf of the Eastern Paratethys in the Oligocene – Neogene. Stratigraphy and Geological Correlation 18(2): 200–224.
  • Popov SV, Rostovtseva YV, Pinchuk TN, Patina IS, Goncharova IA (2019) Oligocene to Neogene paleogeography and depositional environments of the Euxinian part of Paratethys in Crimean – Caucasian junction. Marine and Petroleum Geology 103: 163–175.
  • Poulakakis N, Lymberakis P, Valakos E, Pafilis P, Zouros E, Mylonas M (2005) Phylogeography of Balkan wall lizard (Podarcis taurica) and its relatives inferred from mitochondrial DNA sequences. Molecular Ecology 14: 2433–2443.
  • Psonis N, Antoniou A, Kukushkin O, Jablonski D, Petrov B, Crnobrnja-Isailović J, Sotiropoulos K, Gherghel I, Lymberakis P, Poulakakis N (2017) Hidden diversity in the Podarcis tauricus (Sauria, Lacertidae) species subgroup in the light of multilocus phylogeny and species delimitation. Molecular Phylogenetics and Evolution 106: 6–17.
  • Psonis N, Antoniou A, Karameta E, Leaché AD, Kotsakiozi P, Darriba D, Kozlov A, Stamatakis A, Poursanidis D, Kukushkin O, Jablonski D, Crnobrnja–Isailović J, Gherghel I, Lymberakis P, Poulakakis N (2018) Resolving complex phylogeographic patterns in the Balkan Peninsula using closely related wall-lizard species as a model system. Molecular Phylogenetics and Evolution 125: 100–115.
  • Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA (2018) Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. Systematic Biology 67(5): 901–904.
  • Ratnikov VY (2015) Squamate reptiles from Upper Pleistocene deposits of Emine-Bair-Khoar Cave in the Crimea. In: Turbanov IS, Marin IN, Gongalsky KB (Eds) Biospeleology of the Caucasus and other regions of Russia. Proceedings of All-Russian conference of young scientists, Moscow (Russia), December 2015. Kostromskoy pechatnyi dom, Kostroma, 64–65. (In Russian).
  • Rato C, Stratakis M, Sousa-Guedes D, Sillero N, Corti C, Freitas S, Harris DJ, Carretero MA (2020) The more you search, the more you find: Cryptic diversity and admixture within the Anatolian rock lizards (Squamata, Darevskia). Zoologica Scripta 50: 193–209.
  • Ryabinin DM, Grechko VV, Darevky IS, Ryskov AP, Semenova SK (1996) Comparative study of DNA repetitive sequences by means of restriction endonucleases among populations and subspecies of some lacertid lizard spcies. Russian Journal of Herpetology 3(2): 178–185.
  • Salvi D, Catarina Pinho C, Harris DJ (2017) Digging up the roots of an insular hotspot of genetic diversity: decoupled mito-nuclear histories in the evolution of the Corsican-Sardinian endemic lizard Podarcis tiliguerta. BMC Evolutionary Biology 17:63.
  • Schmitt S, Pouteau R, Justeau D, de Boissieu F, Birnbaum P (2017) SSDM: An R package to predict distribution of species richness and composition based on stacked species distribution models. Methods in Ecology and Evolution 8: 1795–1803.
  • Schulte JA, Losos JB, Cruz FB, Núñez H. (2004) The relationship between morphology, escape behaviour and microhabitat occupation in the lizard clade Liolaemus (Iguanidae: Tropidurinae: Liolaemini). Journal of Evolutionary Biology 17(2): 408–20.
  • Schule W (1993) Mammals, vegetation and the initial human settlement of the Mediterranean islands: a palaeoecological approach. Journal of Biogeography 20: 399–412.
  • Senczuk G, Colangelo P, De Simone E, Aloise G, Castiglia R (2017) A combination of long term fragmentation and glacial persistence drove the evolutionary history of the Italian wall lizard Podarcis siculus. BMC Evolutionary Biology 17(6): 1–15.
  • Shalimov AI (1966) Novel tectonic scheme of Crimea and connection of folded structures of the Mountainous Crimea and North-Western Caucasus. In: The structure of the Black Sea depression: Collection of the papers. Nauka, Moscow, 49–58. (In Russian).
  • Shnyukov EF, Shcherbakov IB, Shnyukova EE (1997) Paleoisland arс in the northern part of the Black Sea. National Academy of Sciences of Ukraine, Kiev, 288 pp. (In Russian).
  • Sirenko OA (2019) Changes in Pleistocene vegatation and climate of Ukraine in the range of 1.8–0.4 million years. Journal of Geology, Geography and Geoecology 28(2): 355–366.
  • Soberón J, Peterson AT (2005) Interpretation of models of fundamental ecological niches and species’ distributional areas. Biodiversity Informatics 2: 1–10
  • Speybroeck J, Beukema W, Dufresnes C, Fritz U, Jablonski D, Lymberakis P, Martínez-Solano I, Razzetti E, Vamberger M, Vences M, Vörös J, Crochet P-A (2020) Species list of the European herpetofauna – 2020 update by the Taxonomic Committee of the Societas Europaea Herpetologica. Amphibia-Reptilia 41(2): 139–189.
  • Spilani L, Bougiouri K, Antoniou A, Psonis N, Poursanidis D, Lymberakis P, Poulakakis N (2019) Multigene phylogeny, phylogeography and population structure of Podarcis cretensis species group in south Balkans. Molecular Phylogenetics and Evolution 138: 193–204.
  • Szczerbak NN (1962) On the systematics of Lacerta saxicola Eversmann of the Crimea and North Caucasus. Zoologicheskyi zhurnal 41(9): 1374–1385. (In Russian).
  • Szczerbak NN (1966) Amphibians and reptiles of Crimea (Herpetologia Taurica). Naukova dumka, Kiev, 240 pp. (In Russian).
  • Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123: 585–595.
  • Tamar K, Carranza S, Sindaco R, Moravec J, Trape JF, Meiri S. (2016) Out of Africa: Phylogeny and biogeography of the wide-spread genus Acanthodactylus (Reptilia: Lacertidae). Molecular Phylogenetics and Evolution 103: 6–18.
  • Tarkhnishvili DM (2012) Evolutionary history, habitats, diversification, and speciation in Caucasian rock lizards. In: Jenkins OP (Ed.) Advances in Zoology Research. Volume 2. Nova Science Publishers, Hauppauge, 79–120.
  • Tarkhnishvili D, Gabelaia M, Mumladze L, Murtskhvaladze M (2016) Mitochondrial phylogeny of the Darevskia saxicola complex: two highly deviant evolutionary lineages from the easternmost part of the range. Herpetological Journal 26: 175–182.
  • Tarkhnishvili D, Murtskhvaladze M, Gavashelishvili A (2013) Speciation in Caucasian lizards: climatic dissimilarity of the habitats is more important than isolation time. Biological Journal of the Linnean Society 109: 876–892.
  • Terrasa B, Picornell A, Castro JA, Ramon MM (2004) Genetic variation within endemic Podarcis lizards from the Balearic Islands inferred from partial Cytochrome b sequences. Amphibia-Reptilia 25: 407–414.
  • Vakhrushev BA, Amelichev GN (2001) To an issue of the possibility of glaciation of the Crimean mountains. Fizicheskaya geografiya i geomorfologiya 40(1): 1–11. (In Russian).
  • Ved’ IP (2000) Climatic atlas of Crimea. Tavriya-Plus, Simferopol, 120 pp. (In Russian).
  • Vremir M, Ridush B (2005) The Emine-Bair-Khosar “Mega-Trap” (Ukraine). Mitteilungen der Kommission für Quartärforschung Österreichischen Akademie Wissenshaften 14: 235–239.
  • Vörös J, Mikulíček P, Major Á, Recuero E, Arntzen JW. (2016) Phylogeographic analysis reveals northerly refugia for the riverine amphibian Triturus dobrogicus (Caudata: Salamandridae). Biological Journal of the Linnean Society 119(4): 974–991.
  • Wielstra B, Crnobrnja-Isailović J, Litvinchuk SN, Reijnen BT, Skidmore AK, Sotiropoulos K, Toxopeus AG, Tzankov N, Vukov T, Arntzen JW. (2013) Tracing glacial refugia of Triturus newts based on mitochondrial DNA phylogeography and species distribution modeling. Frontiers in Zoology 10(1): 13.
  • Xia X (2018) DAMBE7: New and improved tools for data analysis in molecular biology and evolution. Molecular Biology and Evolution 35: 1550–1552.
  • Xia X, Lemey P (2009) Assessing substitution saturation with DAMBE. In: Lemey P, Salemi M, Vandamme A-M (Eds) The Phylogenetic Handbook: a Practical Approach to Phylogenetic Analysis and Hypothesis Testing. Cambridge University Press, 615–630.
  • Zinenko O, Stümpel N, Mazanaeva L, Bakiev A, Shyryaev K, Pavlov A, Kotenko T, Kukushkin O, Chikin Y, Duisebajeva T, Nilson G, Orlov NL, Tuniyev S, Ananjeva NB, Murphy RW, Joger U (2015) Mitochondrial phylogeny shows multiply independent ecological transition and northern dispersion despite of Pleistocene glaciations in meadow and steppe vipers (Vipera ursinii and Vipera renardi). Molecular Phylogenetics and Evolution 84: 85–100.

Supplementary materials

Supplementary material 1 

Figures S1–S6

Kukushkin O, Ermakov O, Gherghel I, Lukonina S, Svinin A, Doronin I, Simonov E, Jablonski D (2021))

Data type: .pdf

Explanation note: Supplementary illustrative materials.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (1.54 MB)
Supplementary material 2 

Tables S1–S3

Kukushkin O, Ermakov O, Gherghel I, Lukonina S, Svinin A, Doronin I, Simonov E, Jablonski D (2021))

Data type: .xls

Explanation note: Table S1. A list of Darevskia tissue samples and cytochrome b sequences used in analysis. — Table S2. GenBank accession numbers of cytochrome b sequences used for the divergence time estimates and chronogram construction. — Table S3. Pearson correlation table shows high correlated variables (│r│> 0.70)

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (91.00 kb)
login to comment