The mitochondrial phylogeography of the Crimean endemic lizard Darevskia lindholmi (Sauria, Lacertidae): Hidden diversity in an isolated mountain system

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 differ­ entiated 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. Vertebrate Zoology 71, 2021, 559–576 | DOI 10.3897/vz.71.e62729 Copyright Oleg Kukushkin et al.. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Oleg Kukushkin et al.: The phylogeography of the Crimean endemic Darevskia lindholmi 560


Introduction
The genus Darevskia Arribas, 1999 includes more than 30 species of mediumsized lizards with ranges locali zed 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 (Evers mann, 1834), D. brauneri (Méhely, 1909), andD. 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 an cestor, 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 Doronin et al. 2021). Based on morphological data, D. lindholmi showed a close relationship to D. brauneri (Szczerbak 1962;Doro nin et al. 2013), and moderate values for genetic distanc es based on cytochrome b between the Crimean and the Caucasian representatives of the complex (3.4-5.3%; Doronin et al. 2013) indicate a PlioPleistocene time of divergence (see Tarkhnishvili et al. 2016). The genetic distances between most bi sexual Darevskia species are much higher (exceeding 10% ;Tarkhnishvili 2012;Ah madzadeh 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 expect ed 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 km 2 (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 pres ence was found here (Vakhrushev and Amelichev 2001). Because of a limited scale and relatively simple topogra phy, along with longterm connections with the adjacent continent, Crimea acted as a "melting pot" rather than re fugia 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 di versity (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. 2017Psonis et al. , 2018Jablonski et al. 2019). In addi tion to the Crimean member of the genus Darevskia, the only endemic reptile taxon to the Crimean Mountains is Lacerta agilis tauridica Suchow, 1927, whose morpho logical uniqueness was supported also by genetic data (Joger et al. 2007;Andres et al. 2014;Kukushkin et al. 2020). Lacerta viridis magnifica Sobolevssky, 1930, for mer 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 sug gested that this species' ancestor was isolated after the disappearance of the mountainforest bridge that directly connected the Crimean Mountains with the northwest ern extremity of the Greater Caucasus during one of the cool and wet periods of the Early Pleistocene (Szczerbak 1962(Szczerbak , 1966Darevsky 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 hypoth esized that the colonization of Crimea by a D. lindholmi ancestor could be associated with its passive transmarine distribution in the Early Pleistocene. However, the ques tion 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 pro vide a biogeographical and taxonomic context.

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 locali ties 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).

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 us ing BioEdit v. 7.0 software (Hall 1999). The genetic vari ability analysis was based on the calculation of the fol lowing 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), di versity 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 analy ses, the dataset was tested for redundancy and saturation using DAMBE 6.4.101 software (Xia 2018) by calcu lating the entropybased 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 sub stitution saturation and suitability of the data for phylo genetic inference. The IQTREE 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 IQTREE 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). IQTREE was ran with the HKY+F+I+G4 substitution model for the 1 st and 2 nd co don positions combined and TIM2+F+G4 for the 3 rd co don 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 Partition Finder 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 6 gener ations 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-likeli hood, 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 re construction 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 (pdistances) 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, mil lions 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 be tween Gallotia and Psammodromus algirus (emergence of Canary Islands): normal distribution, mean 18.0±2.0 (Cox et al. 2010;Carranza and Arnold 2012); the diver gence 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 separa tion of Crete from the Peloponnese leading to the diver gence 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 analy sis was run for 10 7 generations with a sampling frequency 1000 generations, from which 10% were discarded as burnin. A relaxed uncorrelated lognormal clock model, Yule process of speciation, and random starting tree were applied. All branches involved in the respective calibra tion 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 fol low the results of phylogenetic inference. The analysis was repeated three times and parameter log files and the phylogenetic trees were combined using LogCombiner

Modeling of the potential ranges and ecological niche analysis
The potential distribution estimates for the main evolu tionary lineages of D. lindholmi were calculated using eight algorithms (Classification Tree Analysis, Multi variate 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 ap proaches 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 Moun tains (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 in clusion of each D. lindholmi population in any given phy logenetic lineage using this analysis was based on the re sults 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 occurrenc es) 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;Pa pes and Gaubert 2007). To test the accuracy of the final stacked models, we used an independent test dataset com prising of 25% of occurrences to calculate the omission error, AUC, and Kappa metrics for each distribution esti mate 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 dis tribution models concerning the environments in which the species lives (Mouton et al. 2010, JiménezValverde 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 concern ing 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 vari ables 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 bio logically 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; coef ficient 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 oc currence 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 presentday model for the species on three paleoclimatic global circulation models (CCSM4, MI ROC, and MPI) based on CMIP5 (Coupled Model Inter comparison 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 suit ability).
To further test the unicity of the ecological niches oc cupied by the main D. lindholmi mtDNA lineages (see Results) we performed two niche overlap metrics to as sess the degree of observed ecological overlap between each lineage by performing Schoener's D (Warren et al. 2008) and standardized Hellinger'sbased 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'sbased I is based purely on probability distri butions 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 distri bution 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 convexhullbased approach in the ArcGIS Spatial Tools toolkit that was used to select the background from which the pseudopresences were randomly generated. The ob served 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 phy logenetic 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 pdistances were between the Com Simultaneously, the latter value of pdistance was estab lished as a maximum in the pairwise comparison among all members of D. saxicola complex. The minimum pdistance (2.7±0.5%) within this species complex was found between D. szczerbaki and D. saxicola. The average value of uncorrected pairwise genetic dis tances 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 South western and Common lineages (Table 1) 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 east ern and northern parts of the Crimean Mountains includ ing coastal, lowmountain, and foothill regions (A), the eastern (B) and western (E) parts of the mountainforest Crimea, and the northwestern foreststeppe 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).    The inter model 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 evolu tionary 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 pre cipitation during the driest month (bio14). The response curves of the three lineages to these three environmen tal 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.

Time of the divergence
To disentangle the direction of the relationship between the probability of distributions of each lineage and en vironemntal variables, we found positive correlations between the probability of presence of the lineages and isothermality ( The potential ranges based on the fundamental niches of the Common, Southwestern, and Central lineages suggests the possibility of some poten tial 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 spe cies had its refugia exclusively in the Crimean Mountains (Fig. 6).
To test the unicity of the niches modeled in relation to the available environmental conditions we performed a background analysis where we found significant dif ferences between the niches occupied by the three line ages. The observed niche overlap between the lineages was very high (Central vs. Common: Schoener's D=0.82,Hellinger'sbased I=0.97;Central vs. Southwestern: Schoener's D=0.80,Hellinger'sbased I=0.97;Common vs. Southwestern: Schoener's D=0.90,Hellinger'sbased I=0.98). However, due to small gradients, these overlaps were still shown to represent different niches based on the identity background analysis (ttest, df =99, P<0.05). The  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.

Phylogeography
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 la ter geological periods (Shalimov 1966;Shnyukov et al. 1997;Esin et al. 2018;Popov et al. 2010Popov et al. , 2019Stovba et al. 2020;Palcu et al. 2021), and limited fossil evidence (Vremir and Ridush 2005;Ratnikov 2015;Kovalchuk et al. 2020). Current paleogeographic reconstruction sug gest that during the PlioPleistocene transition (the Ak chagylian epoch, 3.6-1.8 Mya), the Black Sea and the CaspianAzov marine basins were separated by land lo cated approximately in the Kerch -Taman region or even southlying 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 landbridge, 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 com plex (Doronin et al. 2013;Tarkhnishvili et al. 2016). The Central lineage is sister to the Common and the South western lineages with uncorrected pdistance 4.6±0.6%, suggesting a longterm independent evolution related to topographical and climatic changes in the region, proba bly during the Early Pleistocene (see Babak 1959;Pisare va et al. 2019;Sirenko 2019).
Presentday deepwater 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 sub sidence of the Black Sea hollow occurred, which was accompanied by the sinking of peripheral parts of sur rounding mountain systems. The subsidence of the south ern "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 fluc tuations as well as the creation of landbridges 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 as sociated with substantial cooling during the MiddleLate Pleistocene (Sirenko 2019). The closest ancestral haplo types of D. lindholmi represented by the Southwestern lineage probably persisted in the extreme south of Crimea (potential refugium), where the influence of climatic os cillations was less pronounced (Gerasimenko 2007(Gerasimenko , 2011Pisareva 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 South western lineages might simply represent geographically structured diversity within a single population that was never split into allopatric areas of vicariance but was sub jected to spatially restricted gene flow due to isolation by distance and habitat fragmentation.
The range dynamics of D. lindholmi distribution pat tern 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 sug gest that during the Last Glacial Maximum (21 Kya) the species inhabited the Crimean Mountains (Fig. 6). As D. lindholmi prefers rocky microhabitats, using niche con servatism theory we expect that the species would have inhabited similar habitats during the Last Glacial Maxi mum in the warmest and wettest areas, such as large rock massifs and warm river canyons of the Crimean Moun tains. 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 genet ic diversity due to extinciton events of local populations. The observed structure of the Common lineage (Fig. 3) might be the result of a rapid postglacial recolonization of the mountain areas from several refugia located on the southern macroscope of the Crimean Mountains. The starlike pattern in the haplotype network, as well as neu trality tests, indicate past dispersal events of the Common lineage, possibly related to sealevel changes during the MidHolocene, 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 northsouth axis and lies eastward from Yalta and the Simferopol meridian and westward to the vast KarabiYaila Plateau (Fig. 1). Thus, the Central lineage distribution is confined to the highest and wet test 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 nar row strip of the precipitous shore below 750 m a. s. l. (i.e. area with the most pronounced subMediterranean character), whereas the Common lineage occupies the rest of the species range (Fig. 1). The altitudinal diapason of the most diverged evolutionary lineages' habitats over laps significantly. Nevertheless, the average elevation of their representative's observations has an almost twofold 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 coldtolerant spe cies. Its distribution at low altitudes is limited by insuf ficient 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 ar eas 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 physicalgeographical areas of the Crimean Mountains and dwell in a wide range of habi tats 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 per sonal data). These observations are congruent with our results obtained from the ecological niche modeling (Fig. 5) which confirms the assumptions based on in situ ob servations.
Our models show significant variation in the response of each lineage to the environmental variables used to model their distributions. In general, three variables ac counted for almost twothirds 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 pre dictors 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 cli matic 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 particu larly high and low temperatures. The ecological niches occupied by the considered groups are not identical, al though show a high degree of overlap. The divergence of ecological niches between the studied populations might occur because of their different attitudes towards the hu midification factor (Fig. S2), though there is a wide dia pason 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 pre vious work, closely related species tend to occupy sim ilar ecological niches, while divergence can occur when significant shifts in habitat use are followed by morpho logical 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 line ages 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 ex amples 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 envi ronmental conditions (CaeiroDias 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 pres ence 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 line ages. 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. dobrogi cus range in two with relatively strong levels of admix ture 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 Papeş 2015, Vörös et al. 2016). Similarly, in the case of the D. lindholmi lineages, we suspect a similar process of track ing suitable ecological niches of the lineages has occurred throughout the Quaternary glacial oscillations which re sulted in the observed contemporary patterns of genetic structuring in the species (Figs 1 and 3).

Taxonomy
The uncorrected pdistance 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;Pso nis 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. bithyni ca (Méhely, 1909). On the other hand, the taxo nomic position of this species belonging to the D. rudis complex is supported by differences in biparental mark ers 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 sub species rank. In addition to mtDNA data, the observed parapatry of the lineages, differences between the nich es, and the existence of certain morphological differenc es between them (O. Kukushkin and A. Svinin personal data;Figs S5 and S6), together with wellknown agree ment between mtDNA data and taxonomy in the west ern 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 his tory inside D. lindholmi implies that ranges between lin eages could be in contact, which predicts potential gene flow (see hybridization between closely related Lacer tidae species in areas of sympatry; Tarkhnishvili et al. 2013;CaeiroDias et al. 2021). Moreover, further inves tigation 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 dis covered 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, particu larly 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.