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.

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 1^st and 2^nd codon positions combined and TIM2+F+G4 for the 3^rd 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 1^st codon position, HKY+I for the 2^nd codon position, and GTR+G for the 3^rd codon position. We conducted two simultaneous runs of four Markov Chains Monte Carlo (MCMC) tests, each run consisted of 4*10⁶ 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).

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 10⁷ 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).

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 (www.worldclim.org; 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.

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.

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.

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.

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 (D_H0 is 0.85±0.02 vs. D_H1 = 0.82, and I_H0 is 0.98±0.007 vs. I_H1 = 0.97) occupy different niches, and the Central lineage and Southwestern lineage (D_H0 is 0.90±0.06 vs. D_H1 = 0.80, and I_H0 is 0.98±0.02 vs. I_H1 = 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 (D_H0 is 0.88±0.06 vs. D_H1 = 0.90, and I_H0 is 0.88±0.07 vs. I_H1 = 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 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).

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.