Species Distribution Models and Niche Partitioning among Unisexual Darevskia dahli and Its Parental Bisexual ( D. portschinskii , D. mixta ) Rock Lizards in the Caucasus

: Among vertebrates, true parthenogenesis is known only in reptiles. Parthenogenetic lizards of the genus Darevskia emerged as a result of the hybridization of bisexual parental species. However, uncertainty remains about the mechanisms of the co-existence of these forms. The geographical parthenogenesis hypothesis suggests that unisexual forms can co-exist with their parental species in the “marginal” habitats. Our goal is to investigate the inﬂuence of environmental factors on the formation of ecological niches and the distribution of lizards. For this reason, we created models of species distribution and ecological niches to predict the potential geographical distribution of the parthenogenetic and its parental species. We also estimated the realized niches breadth, their overlap, similarities, and shifts in the entire space of predictor variables. We found that the centroids of the niches of the three studied lizards were located in the mountain forests. The “maternal” species D. mixta prefers forest habitats located at high elevations, “paternal” species D. portschinskii commonly occurs in arid and shrub habitats of the lower belt of mountain forests, and D. dahli occupies substantially an intermediate or “marginal” position along environmental gradients relative to that of its parental species. Our results evidence that geographical parthenogenesis partially explains the co-existence of the lizards. and D. mixta ( t = 7.13, p << 0.001). To conclude, the settings by default resulted in less accurate SDMs.


Introduction
True parthenogenesis is very rare in vertebrate animals and less than 0.46% of reptile species exhibits it [1]. In true parthenogenesis, female gametes develop into embryo without fertilization. Some species of reptiles are the only representatives of the chordates which are reproduced by parthenogenesis, and for this reason they can be objects for studying the evolution of true parthenogenesis [2][3][4]. In some reptile species, parthenogenetic lineages are produced by interspecific hybridization between closely related bisexual species [3,5]. Despite the fast rate of unisexual reproduction in general, the growth of parthenogenetic forms is constrained by accumulation of harmful mutations and their poor adaptation to parasites [4,6]. Nevertheless, parthenogenetic lineages can be used as exclusive objects for identifying key mechanisms that ensure long-term preservation of unisexual forms despite their long co-existence of lizards [39][40][41][42][43][44]. However, these studies were mainly based on bioclimatic predictors using maximum entropy (Maxent) [45,46] whereas topographic and landscape variables were rarely applied.
In the present study, we have increased the number of occurrence records and involved more predictor variables to establish the "marginal" habitats of the parthenogenetic form D. dahli. We applied species distribution models (SDMs) and ecological niche models (ENMs), which are commonly used in theoretical and applied studies in ecology and biogeography [47]. The applied models are of particular interest because they are based on two distinct approaches answering the different questions [48]. Species distribution models attempt to estimate species distribution over the geographic space while ENMs estimate fundamental niches of species, and can be applied to predict a potential distribution and overlapping of niches in time and space [49]. We used maximum entropy distribution modelling (SDM Maxent) to predict the potential ranges of the studied species and the ordination method to assess the overlapping, similarities and shifts of niches in the entire space of predictor variables. We chose Maxent as the most appropriate modelling approach for the case of continuous and categorical predictor variables which can be very effective with a relatively small number of occurrence records [39,[50][51][52]. In contrast to ordination methods, Maxent enabled us to select and rank variables in accordance to their importance for detecting species ranges using optimal model parameters based on the Akaike information criterion [53,54]. Noteworthy is the fact that MaxEnt has previously been successfully applied to construct SDMs for some lizard species using large or limited occurrence records [31,[39][40][41][42][43][44]. This approach can be applied not only to any unisexual and bisexual species but also to their complexes.
In this study, we used a set of original, museum and published data on the distribution of D. portschinskii, D. mixta, and D. dahli to: (1) identify bioclimatic, topographic, and landscape variables that determine their potential habitats, (2) construct maps of their potential geographical distribution, (3) estimate the overlapping, similarities, and differences of their ecological niches, and (4) verify whether the hypothesis of GP is true for D. dahli.

Materials and Methods
In this study, we used original data from fieldwork as well as published and museum data of species (D. dahli, D. mixta, and D. portschinskii), verified in 2018-2019. The spatial resolution of the predictor layers was 3 arc s (~100 m). The data analysis was conducted using five consecutive steps: (1) dataset acquisition and preparation of vector and raster layers; (2) spatial thinning of record sites and predictor variables; (3) determination of the parameters of the MaxEnt models using AIC information criterion; (4) construction of species distribution models (SDMs) using the method of maximum entropy and realized ecological niches models (ENMs) in the space of first two axes of the principal component analysis; (5) comparison of ecological niches (breadth, similarity, overlapping, shifts) and quantitative parameters of habitat exploitation. This approach can be applied to the other unisexual taxa.

Dataset Acquisition and Preparation of Vector and Raster Layers
We created vector layers of the locality records of studied lizards using original data, collected during the field surveys in June and July of 2018 and 2019. The new records allowed us to make corrections of the coordinates of lizard occurrence known from published data and museum collections. Further we referred to new records as the "original" field records of species occurrence. Fieldwork (2018-2019) was carried out in the habitats of the species located on rock outcrops, individual blocks of stone, large stones and clay cliffs in mountain-steppe, mountain meadow, mountain forest, subalpine and alpine zones. The total length of the surveyed route was 4800 km. The geographical coordinates and elevations of locality records were determined by Garmin Montana 680 t GPS receiver (Garmin Corp., Olathe, KS, USA). Geographic coordinates were determined with an error of ± 3.5 m. In field studies lizards were captured by noose in each locality and identified using the species guides [2,30]. After identification, all the lizards were immediately released in the same place where they were captured. We took photographs of each lizard's anterolateral and temporal areas of the head and its anal areas with a Nikon Coolpix B500 digital camera for subsequent control of species identification in the laboratory.
Thus, we collected 487 records with the geographical coordinates of the three studied lizards occurrence data in Armenia, Georgia, Nagorno-Karabakh, and Azerbaijan (D. dahli-164; D. portschinskii-258; D. mixta-67). Although many lizard habitats were previously known [2,25,29,31,34,[55][56][57], we collected new records during our field surveys in 2018-2019 to enlarge dataset of lizard occurrence and to specify geographical coordinates of lizard habitats which were described in published data and museum collections without coordinates. Using original field records and data from museums and published articles, we built vector layers of lizards occurrences in ArcGIS Desktop 10.6.1.
We analysed available published data to select the most important variables that determine the lizards distribution [2,25,29,31,[42][43][44]. All predictor variables including bioclimatic, topography, and land cover/land use are presented in Supplementary Table S1. Spatial standard bioclimatic variables (19) were taken from the WorldClim 2.1 dataset [58,59], which characterize annual trends in seasonality and temperature and precipitation variations. In addition, we included 2 variables for insolation and wind speed [58]. The raster layers of the relief were taken from the freely distributed dataset of shuttle radar topographic mission (SRTM) [60]. The predictor variables also comprised raster layers of elevations, inclination angles, and exposure aspects derived from SRTM using ArcGIS Desktop 10.6.1. Data on land use, roads and railways, human settlements, rivers, vegetation types, soils in Armenia, Georgia and Azerbaijan was obtained from the Open Street Map resource [61,62].

Spatial Thinning of Records and Predictor Variables
We used a two-step procedure to verify the spatial autocorrelation of species occurrence points. First, the occurrence points were separated by a distance less than 1 km between each other using the subsample algorithm available in the spThin R package [63]. Then the dataset was tested by cluster analysis using the average nearest neighbor index in ArcGis 10.6.1 [64]. After applying this procedure, we excluded autocorrelated points in the datasets for D. dahli, D. mixta, and D. portschinskii (Supplementary Table S2 and Figure S1).
We selected predictor variables for the model approach conducting two step analysis. First, using the MaxEnt 3.4.1 package [45,46], we generated 10,000 randomly distributed background points for each lizard species individually. Then with the help of corSelect functions of the fuzzySim package [65], we reduced the set of variables between which the absolute value of Pearson pair correlation coefficient was >0.75. Multicollinearity was assessed using the variation inflation factor (VIF) in package fuzzySim [65] with a threshold value of VIF > 10 [66]. The reduced list of variables is presented in Table 1.

Determination of the Maxent Models Parameters
Although in the models of spatial distribution, the parameters' settings selected by default in Maxent were based on extensive empirical material [46], recent studies have shown that they may be inefficient [67]. For this reason, we determined the optimal parameters of the Maxent models for each species using the AIC information criterion in the ENMeval R package [54]. The default Maxent parameters were different from those defined using the AIC criterion. Although we did not find any general trends in selection of Feature classes (FC), the FCs selected using the AICc models had higher regularization parameters (RM) than the default value of 1.0 (Supplementary Figure S2).

Construction of Species Distribution (SDM) and Ecological Niche (ENM) Models
Species distribution models for D. dahli, D. portschinskii, and D. mixta were constructed by the maximum entropy method using Maxent 3.4.1 in the Dismo environment [68]. These models were created as a result of 10 Maxent runs to randomly select test and training samples. In all Maxent runs, 80% of the occurrence records were used as training samples while 20% were served as test samples. We used the Boyce index (B ind ) to assess model performance [69][70][71] with the help of the EcoSpat R package [71]. Boyce index lacks those drawbacks which has AUC index [72,73]. It requires only data on lizard occurrence and measures how much the predictive models differ from random distribution. We calculated B ind for all the 10 models for each lizard; afterwards we averaged the obtained values to get the final estimates. The importance of each predictor variable of SDMs was assessed by analysis of Maxents variable contribution table using the jackknife method [45,46]. Variables with essential impact on the model, i.e., having high values of permutation importance (PI > 5) and/or high values of the percent of contribution (PC > 5) were considered as the most important. Niche models (ENMs) of lizards were constructed in the frameworks of the general concept [73,74] which states that ecological niche is represented in a grid of predictor variables of the habitats, depicted by the first two components of the principal component analysis (PCA). In particular, we tied environmental space of the lizards to a grid, and further the lizards' occurrence data converted into densities using the kernel function to smooth out the density distribution. The environmental data were also converted into densities. Thus, habitat preferences of the lizards were estimated based on their occurrence density and environmental densities [71]. We assessed the niche overlap by the method described elsewhere [74,75] using the Schoener's D index. This index measures niche overlap in the ecological space. It varies in the range from 0 (no overlap) to 1 (complete overlap). We used the methodology described elsewhere [71,74] to perform a niche similarity test. Niche similarity was tested using EcoSpat package [71].

Comparison of Ecological Niches and Parameters of Habitat Exploitation
The niches breadth of lizard was estimated by the Mahalanobis distance (MD) [76]. First, we obtained the distance from each locality to the centroid of the lizards using the identical set of predictor variables. Further we found the coefficient of variation of MD (CV MD) for each lizard. The CV MD between species-specific centroids and individual niches of the species was interpreted as a multidimensional quantitative estimate of the niches breadth (N b ). The difference in the niches breadth was estimated in two steps. First, the uniformity of the niches breadth was tested using multiple comparisons [77]. Then, if the null hypothesis was rejected, niches were compared pairwise [78].
Ecological niches of bisexual and parthenogenetic lizards in multidimensional predictor space were compared by distinguishing three types of overlapping: (1) stable areas where lizards occur in both ranges; (2) unfilling areas that are present only in the range of the parental species, and (3) new areas that are present only in the range of the parthenogenetic lizards [31]. We used three indices to quantify these types of overlap with the help of the EcoSpat package [71]. Index of stability (S) is the part of the range of D. dahli that overlaps with that of its parental species. Index of unfilling (U) is the part of the range of a parental species that is absent in the range of D. dahli. Index of expansion (E) is the part of the range of D. dahli which is located in environmental conditions that are different from those of parental species.
In the EcoSpat package, the average positions (centroids) of ecological niches and their shifts along environmental gradients were determined based on the data obtained using the smoothing functions of the densities of the species occurrence and the accessible habitats. Therefore, we had to additionally test the significance of these shifts with the help of the Generalized Linear Model (GLM) procedure. Centroids were compared only on the basis of species occurrence data using one-way analysis of variance (GLM ANOVA) with equal and unequal numbers of replicates in the cells. In all the cases, we used GLM ANOVA models of type I, i.e., models with fixed factors. Pairwise comparison of niche centroids was a factor with three levels. If the GLM ANOVA with fixed effects showed a significant difference between factor (species) levels, a post-hoc Tukey multiple comparison HSD test was used to identify which factor levels were different. If the sample sizes were unequal, i.e., in the case of an unbalanced model of the 1st type, we used a Tukey-Kramer test for a multiple comparison. For a multiple comparison with unequal variances based on the Levene criterion, we used the Tukey-Kramer test with the Welch modification [77]. All data were log10-transformed prior the analyses to achieve normality of residuals and improve homoscedasticity of variance.
Graphic representations of centroid shifts of ecological niches along environmental gradients were obtained with the EcoSpat package [71]. Significance of shifts tested using multcomp package [79] in R 3.3 [80]. The RStudio Desktop version 1.1.463 used as an IDE for R language [81].

Results
After applying the spThin subsampling procedure and the sequential removal of clustering records, 45, 84 and 41 location records remained for D. dahli, D. portschinskii and D. mixta, respectively (Supplementary Table S2 and Figure S1). Further reduced datasets without autocorrelated points were used to build species distribution models and compare ecological niches.
Contribution of each predictor variable in SDMs and estimates of B ind at optimal parameters of Maxent are given in Table 2. These results show that three variables, namely precipitation of the warm quarter (C_PWarmQ), solar irradiation (C_SRad) and elevation (T_EL), are significant for all species. The remaining predictor variables determining the suitable habitats varied among the species. The habitat preferences of D. dahli determined by isothermality (C_ISOT) and distance to road (L_DHW); while the "paternal" species D. portschinskii depends on isothermality (C_ISOT), seasonal coefficient of humidity variation (C_PCoefVar), annual temperature range (C_TAnR), air temperature in dry quarter (C_MTDrQ), distance to road (L_DHW), and vegetation type (L_VEG). The habitat preferences of D. mixta, mainly depend on the seasonal coefficient of humidity variation (C_PCoefVar) and vegetation type (L_VEG). There are three common predictor variables for D. dahli and D. mixta ( Table 2). All the variables responsible for D. dahli preferences are relevant for D. portschinskii except for vegetation type. Supplementary Figures S3-S5 indicated that for all the rock lizards, the curves of dependencies of probabilities of habitat preferences on climatic and topographic variables were bell-shaped with p-values > 0.5. Distance to roads is the most important for the distribution of D. dahli and D. portschinskii (Table 2), while vegetation type is significant only for parental species (Table 2).

Potential Range of Studied Lizards
Models predicted that lizards' ranges mainly cover north-eastern Armenia, western Azerbaijan, and southern, central, and northern Georgia ( Figure 1). The occurrence records and preferable habitats of D. dahli are divided into seven expansive isolated regions located in the highland forest, meadow and steppe zones in the north-eastern part of Armenia (Figure 1a). The range of D. portschinskii is located in southern Georgia, north-eastern Armenia, the western parts of Nagorno-Karabakh, and Azerbaijan ( Figure 1b). The range of D. mixta is found in the south-to-north gorges from the Lesser to the Greater Caucasus in Georgia (Figure 1c).     Estimates of niche overlap between parthenogenetic and its parental species are given in Figure 2. These results showed that the first axis correlated with precipitation in the warm quarter of year while the second axis correlated with elevation. The first and second components accounted for 81.9% of the total variation of predictor variables (Supplementary Figure S6). We did not include more axes because they were associated only with a small part of the total variation. Figure 2 illustrates that the realized niche of D. dahli is shifted upwards along elevation gradient whereas in regard to precipitation of the warm quarter of year, it is slightly lower relative to D. portschinskii. The marginal habitats of D. dahli are characterized by high elevation and by decrease or increase of precipitation of the warm quarter (Figure 2a, Niche expansion). The realized niche of D. dahli is shifted downwards in elevation and precipitation of warm quarter relative to the niche of the "maternal" species D. mixta. All marginal habitats of D. dahli are distinguished from those of D. mixta by decrease of precipitation in the warm quarter (Figure 2c, Niche expansion). Niche overlap data estimated using the Schoener's D index is presented in Table 3. Table 3 shows that D. dahli exploits about 88% of the habitats used by D. mixta. Only 57% of its habitats are located outside the ecological niche of D. mixta. Estimates of similarity between the niches of D. dahli and D. mixta show that there is a statistically significant difference between ecological niches of these species (p = 0.09). However, D. dahli uses 98% habitats occupied by its "paternal" species D. portschinskii. In addition, the range of D. dahli exceeds by 2% that of its "paternal" species. The similarity of ecological niches of D. dahli and D. portschinskii is statistically confirmed (p = 0.009).    Table 3 shows that D. dahli exploits about 88% of the habitats used by D. mixta. Only 57% of its habitats are located outside the ecological niche of D. mixta. Estimates of similarity between the niches of D. dahli and D. mixta show that there is a statistically significant difference between ecological niches of these species (p = 0.09). However, D. dahli uses 98% habitats occupied by its "paternal" species D. portschinskii. In addition, the range of D. dahli exceeds by 2% that of its "paternal" species. The similarity of ecological niches of D. dahli and D. portschinskii is statistically confirmed (p = 0.009).

Shifts of Ecological Niches Centroids along Environmental Gradients
Supplementary Figures S7 and S8 demonstrate the centroid shifts of the ecological niche of the parthenogenetic lizard relative to those of its parental species along environmental variables. The X axis represents a predictor variable whereas the Y axis depicts the occurrence density (defined by EcoSpat) in the multidimensional environmental space. These figures show that the niche centroids of D. dahli are shifted along the axes of many environmental variables. Along the axes (gradients) of isothermality (C_ISOT-BIO3) (Supplementary Figure S7a) and annual temperature (C_TAnR-BIO7) (Figure S7b), the niche centroids of D. dahli are shifted upwards, i.e., "marginal" habitats are characterized by the largest ranges of changes in mean annual and diurnal temperatures. "Marginal" habitats are determined by low air temperature in dry quarter of year (C_MTDrQ-BIO9) (Supplementary Figure S7c), i.e., centroid of the niche along this variable is shifted downwards. Habitats of D. dahli are "marginal" for D. mixta and determined by high seasonal variation of precipitation (C_PcoefVar-BIO15) (Supplementary Figure  S7d). The centroid for this variable is significantly shifted upwards. This means that D. dahli prefers habitats where humidity variation range is greater than of D. mixta. However, the humidity of the habitats during the warmest season (C_PWarmQ-BIO18) (Supplementary Figure S7e) of D. dahli is significantly less than of D. mixta. Supplementary Figure S7f shows that solar irradiation (C_SRad) in habitats preferred by D. dahli is similar or higher than in habitats of D. mixta. The realized niches of D. dahli evidenced that the preferred habitats of the parthenogenetic lizard are shifted along the elevation variable (T_EL) downwards relative to that of D. mixta (Supplementary Figure S7g) due to a smaller number of the marginal habitats occupied by D. mixta than by D. dahli. Preferred habitats of D. dahli are located significantly closer to roads (L_DHW) than in D. mixta (Supplementary Figure S7h). Darevskia mixta mainly inhabits forests while D. dahli dwells in forests, steppes, meadows, semi-deserts, and human settlements (Supplementary Figure S7i).
In regard to isothermality (C_ISOT-BIO3) and annual temperature (C_TAnR-BIO7), D. dahli occupies both "marginal" and regular habitats of the "paternal" species D. portschinskii. The centroid of D. dahli is shifted upwards along isothermal variable (Supplementary Figure S8a Figure S8d). In regard to humidity of warm quarter of year (C_PWarmQ-BIO18)(Supplementary Figure S8e), the centroid of D. dahli is shifted upwards. Nevertheless, D. dahli rarely uses regular habitats of "paternal" species and habitats with lower and higher values of this environmental variable. Supplementary Figure S8f,g suggests that the niche centroids of D. dahli are shifted upwards along variables of solar irradiation (C_Srad) and elevation (T_El). In contrast to its "paternal" species, D. dahli occurs close to roads (L_DHW) (Supplementary Figure S8h), occupies forest, steppe, meadow, and semi-desert habitats, as well as human settlements (Supplementary Figure S8i). Thus, D. dahli uses more diverse habitats than its parental species.

Discussion
Establishment of lizard ecological niches using modelling approach is a powerful tool for studying separation of lizard niches and species-specific requirements for the environmental variables. Moreover, while SDMs predict potential preferable habitats in geographic space (in Gspace), the ordination method (ENMs) enables to perform comparative niche analysis in space of environmental predictor variables (in E-space) [48] for parthenogenetic and parental species to measure their overlap and separation. We used this approach to reliably establish ecological niches and to assess the extent of separation between niches of competing and/or co-existing unisexualbisexual forms of the genus Darevskia [31]. Our data indicate that differences in ecology between the parthenogenetic form and its parental species facilitate the co-existence of these species. Based on our results, we can suggest a set of environmental factors that determine the distribution of the

Discussion
Establishment of lizard ecological niches using modelling approach is a powerful tool for studying separation of lizard niches and species-specific requirements for the environmental variables. Moreover, while SDMs predict potential preferable habitats in geographic space (in G-space), the ordination method (ENMs) enables to perform comparative niche analysis in space of environmental predictor variables (in E-space) [48] for parthenogenetic and parental species to measure their overlap and separation. We used this approach to reliably establish ecological niches and to assess the extent of separation between niches of competing and/or co-existing unisexual-bisexual forms of the genus Darevskia [31]. Our data indicate that differences in ecology between the parthenogenetic form and its parental species facilitate the co-existence of these species. Based on our results, we can suggest a set of environmental factors that determine the distribution of the parthenogenetic form D. dahli and its parental bisexual species (D. portschinskii, D. mixta) in the Caucasus.

Predicted Distribution Range with Optimal Model Parameters
The available occurrence data on the distribution of the studied species in the Caucasus are still not sufficient. Especially small number of records is reported in regard to D. dahli in the northern and western parts of Georgia (Figure 1). The SDM of D. dahli predicted that the northern border of the main range in the west expands along the southern foothills of the Trialeti Range from the village of Tsalki  (Figure 1a). The SDM also predicted the existence of small isolated locations of lizard populations far from the main range in Georgia. Occurrence records and SDM demonstrated that potential range of D. dahli is located in the two north-eastern provinces Lori and Tavush in Armenia, and in western Azerbaijan (Figure 1a). Darevskia dahli distribution in the south-western part of the range predicted by SDM is in a good conformity with the field data. Rare occurrence of the species in the vicinity of Gyumri (40.7942 • N, 43.84528 • E) of the Shirak province is largely related to the lack of suitable habitats at elevations of 1500-1600 m above sea level. Most places where D. dahli was previously recorded were predicted by SDM [2,25,29,30]. The SDMs predicted a wider area of distribution beyond already known range in the Caucasus which has been established based on lizard records. Thus, this map indicates the potential boundaries of isolated lizard populations located in Armenia which requires performing further field surveys there.
Occurrence records and SDM showed that the range of D. portschinskii covers the valleys of the middle reaches of the Kura River in the central and southern Georgia, northern Armenia, and north-western Azerbaijan (Figure 1b). The north-eastern border of the range is constrained by the Kartalinsky Ridge and the valley of the middle reaches of the Iori River, near the village Bochormy in Azerbaijan expand to the border with Armenia, specifically, to the valley of the middle reaches of the Gandzhachaya River. The SDM indicated that suitable habitats for this lizard are available in the Shahumyan province of Nagorno-Karabakh. The occurrence of D. portschinskii in Nagorno-Karabakh was confirmed during our field survey in 2018. Thus, SDM predicted a wider distribution range of D. portschinskii than was known before [2]. In particular, it expands outside the previously known suitable habitats in Georgia north of the city of Gori (42.208 • N, 43.988 • E) and up to the border of Armenia with Azerbaijan.
Darevskia mixta predicted distribution was given elsewhere [31]. This species was used in our previous study of geographical distribution of the parthenogenetic lizard D. armeniaca and its parental species (D. mixta and D. valentini). However, in that study, SDM was constructed using Maxent parameters by default. The new version of SDM was obtained using the optimal values of the Maxent parameters (Figure 1c). Although these two versions of SDMs are different, the main predictions of the lizard distribution in north-western Turkey and in South Ossetia are similar [31,82,83]. Predictions of the adjusted version of SDM and occurrence records in central and southern Georgia show that they are in a good agreement [29].
Finally, although SDMs predicted most habitats where species had been previously recorded, they additionally indicated new distribution areas, e.g., beyond the already known ranges in the Caucasus. This may be attributed to the limited number of reduced occurrence points which were used for the two species: D. dahli (45) and D. mixta (41). Therefore, it requires to carry out further field surveys to clarify the boundaries of ranges, and to get additional occurrence points to improve SDMs. The accuracy of these models can be enhanced by increasing the dataset excluding autocorrelating points following the condition that the sample size (the number of occurrence points) should be ten times larger than the number of predictors used for modelling [84].

Niches of the Unisexual and Bisexual Lizards: Breadth, Overlap, Similarity and Shifts
Ecological niche of D. dahli is different from that of its parental species. However, the niches of D. dahli and its "maternal" species D. mixta differ greater than niches of D. dahli and D. portschinskii. The low overlap index of D. dahli and D. mixta niches and insignificant similarity of their niches can be attributed to the separation of geographical locations of the parthenogenetic and its "maternal" species. Since the distance between the niche centroids of D. dahli and D. mixta was significantly greater than 2 × SD MD of D. mixta, the "maternal" species could not suppress D. dahli. However, the converse statement is not true. This conclusion is also consistent with the fact that the niche breath of D. dahli is significantly greater than in D. mixta. Despite the fact that D. dahli and D. mixta often co-exist in different geographical areas, the index of stability (S) is quite high. This means that D. dahli uses a lot of habitats suitable for the "maternal" species D. mixta besides marginal habitats.
The high indices of stability (S) and their overlap between D. dahli and its "paternal" species D. portschinskii as well as a small expansion of the range of the parthenogenetic form suggest availability of a large number of sympatric zones (SZ). Indeed, 30 SZ of D. dahli and its "paternal" species were identified based on our field data on lizard occurrence in Armenia (2018-2019) and published reports in Georgia [29]. The distance between D. dahli and D. portschinskii niche centroids and the 2 × SD MD for D. dahli suggest that D. dahli can suppress the population development of the "paternal" species. Darevskia portschinskii also has an inhibitory effect on D. dahli. Niche comparison evidences that niche overlap of the "paternal" and parthenogenetic form is asymmetrical, in particular, D. portschinskii occupies more habitats at low elevations than D. dahli.
Estimates of niche shifts of D. dahli relative to its "maternal" species D. mixta are in good agreement with ecological characteristics of these lizards. High values of occurrence density of D. mixta in forest habitats completely agrees with field data, thus supporting our hypothesis that this lizard mainly distributes in forest areas. Darevskia portschinskii prefers habitats located in the lower belt of mountain forests and/or shrub biotopes. Estimates of the occurrence density are in a good agreement with field data.
In contrast to its parental species, vegetation type hardly has impact on distribution of D. dahli because its contribution to SDM accounts for less than 2%. This finding is supported by evidences that D. dahli occasionally occurs in habitats of mountain steppes, mainly on the borders of forests and highways (roads) [25,30]. This means that D. dahli is able to live in habitats with various types of vegetation, specifically in mountain forests, meadows, steppes, and human-transformed habitats.
The ecological plasticity of D. dahli was demonstrated during intentional introduction of D. armeniaca from Armenia to Ukraine in 1963 [85,86]. Among 126 adult females, there was at least one female of parthenogenetic lizard D. dahli [87]. It was enough to create two isolated populations in the Zhytomyr region (Ukraine) [87]. It was shown that morphological traits of the specimens of D. dahli from Armenia and Ukraine were similar. However, there are environmental factors which constrain its further expansion into the range of D. mixta in Georgia. In particular, although D. dahli can dwell in forests, D. mixta prefers forest habitats with higher elevations, colder climate with lower solar irradiation and higher precipitation than D. dahli.
Darevskia armeniaca did not also expand into the range of D. mixta because it cannot tolerate high humidity at low elevations [31]. It is likely that these two parthenogenetic forms (D. dahli, D. armeniaca) do not withstand competition with D. mixta in wet habitats. We hypothesize that competition with D. dahli and D. armeniaca caused the shift of the range border of D. mixta to the west from its native range by displacement of both parthenogenetic forms from their habitats. Thus, as a result of competitive interactions, the contact between their parental species D. mixta and D. portschinskii is impossible which excludes their further hybridization.
Our results show that D. dahli and D. portschinskii rarely co-exist with D. mixta, therefore, sexual interactions between them are unlikely. However, more complex competitive interactions exist between D. dahli and D. portschinskii. The extensive material analysis indicates that triploid hybrids between D. dahli and D. portschinskii were occasionally recorded in Armenia [88]. Unfortunately, due to the rare hybridization event between D. dahli and D. portschinskii as well as the lack of sufficient material, it is impossible to assess the impact of hybridization on the genetic pool of D. dahli.

Mechanisms of Coexistence of Unisexual and Bisexual Forms
There are several hypotheses which suggest mechanisms of co-existence of bisexual and unisexual forms in the frameworks of the GP [11]. The "General Purpose Genotype" hypothesis (GPG) and the "Frozen Niche Variation" hypothesis (FNV) are the most popular. Although the GPG and FNV models are often regarded as mutually exclusive hypotheses, they have much in common. In addition, they are focused on environmental changes over time and space [9].
Hypothesis of GPG admits that the following conditions are fulfilled: (1) the individual features of the bisexual population vary from narrow to wide tolerance ranges so that their genotypes have specific structures to survive under extreme environmental conditions; (2) the range of potential genotypes is frozen among the clones produced by bisexual parental species; (3) natural selection will act in favor of clones with a wider tolerance. Therefore, GPG model assumes the existence of a multiclonality of parthenogenetic forms because natural selection can act more efficiently towards polyphyletic clones [89]. Such a feature is very important for long-living clonal forms because the spontaneous mutagenesis of a single parthenogenetic lineage can have negative effects in the long run due to mutational meltdown [90].
Our original and available published data evidenced that GPG model showed a reliable prediction in case of D. dahli. In particular, (1) variation of three microsatellite loci for 111 specimens of D. dahli collected from five different populations in Armenia suggest that the clonal diversity and creation of its clones are associated with multiple hybridization events [23]; (2) Darevskia dahli has intermediate requirements for habitat variables relative to those of its parental species, but some variables including isothermality, solar irradiation, seasonal variation in humidity and average temperature in the dry quarter of year strongly differ from its parental species; (3) ecological plasticity of the unisexual form has been proved experimentally [87,91], e.g., D. dahli successfully naturalized in the invasive part of the range in Ukraine.
The hypothesis of FNV suggests that the competitive interactions between unisexual and bisexual forms can result in: (1) stable co-existence of bisexual species with genetically heterogeneous clones while clones with large niche overlap with bisexual species are eliminated by natural selection; (2) competitive exclusion of a parental species by clones with high reproductive rate [9]. This model also indicates that parental species have a wider ecological niche than clones. Therefore, a bisexual population can avoid displacements by some clones because, for example, they are able to use resources which are not accessible to clones. Their co-existence can last as long as the niche of a parental form is wider than that of a clone [92].
Thus, the FNV model can be a useful tool to describe competitive interactions between lizards. The niches breadth of parthenogenetic forms D. dahli is significantly larger than its "maternal" species D. mixta. For this reason, FNV model suggests that bisexual species are not able to co-exist with a large number of parthenogenetic lineages which can suppress them [92][93][94]. This finding is supported by the fact that D. mixta was excluded from its native range by a large number of clones of D. armeniaca and D. dahli [23,95]. It is noteworthy that D. mixta is the only "maternal" species of the genus Darevskia that competed with multiple parethenogenetic clones of D. armeniaca and D. dahli simultaneously, since niche elevation centroid of D. mixta is located at the intermediate position between centroids of D. armeniaca (top) [31] and D. dahli (bottom).
Nevertheless, GPG and FNV hypotheses do not involve other mechanisms of evolutionary change in phenotypic diversity of parthenogenetic forms [96,97]. The polyphyletic structure of D. dahli provides high phenotypic diversity, thus facilitating the survival of the parthenogenetic form under the effects of environment changes. It is likely that D. dahli co-exists with the D. portschinskii partially due to its phenotypic plasticity.
Overall, GP model involves mainly polyphyletic clones consisting of several unisexual lineages which emerge as a result of one or multiple events of hybridization. Thus, GP only partially explains co-existence of long-living clonal forms suggesting that there can be other post-mutational and environmental mechanisms which ensure lizard co-existence. Phenotypic diversity of the lizard clones largely depends on post-mutation events, environmental conditions, age of lizards, and size of the range [23,32,[96][97][98][99].

Conclusions
This study increased our knowledge about the parthenogenetic rock lizard and its parental species distribution in the Caucasus. In particular, (1) we supplemented new field data to the available lizards occurrence records and (2) constructed two types of models, specifically SDMs for D. dahli and its parental species (D. portschinskii and D. mixta) and resized ENMs using the ordination method to estimate the breadth of lizard niches, their overlap, similarities, and shifts in the entire space of predictor variables. High indexes of model performance (B ind ) for all studied lizards calculated with optimal Maxent parameters using AICc criterion evidenced that we were able to select the most important environmental predictor variables that determine habitat suitability for lizards. The narrow niche breadth of the "maternal" species D. mixta relative to that of D. dahli, and separation of their habitats confirm that the initial assumption of GP model was fulfilled. However, the fact of the displacement of D. mixta from its native range by polyphyletic clones of D. dahli has not been completely proven and requires further studies. On the other hand, the availability of polyphyletic clones, a significant superiority in the niche breadth of the "paternal" species D. portschinskii, and significant shifts of the niche centroids of clonal forms facilitate co-existence of these species in the Caucasus. We regard the differentiation of unisexual and bisexual lizards' niches as a mechanism of their survival. Finally, we developed a new methodological approach based on the SDMs and ENMs which can be further applied for studying the niche partitioning in unisexual and its parental bisexual forms. These results can be helpful for conducting future field surveys and can be used by environmental agencies and/or decision makers to preserve natural habitats for rock lizards.
Supplementary Materials: The following are available online at http://www.mdpi.com/2227-7390/8/8/1329/s1, Figure S1: Locations of the study areas based on the available data sets, where (a), (c), (e) are initial clustered data sets; (b), (d), (f) are reduced non-autocorrelated data sets. Dotted areas represent masks used to fit the potential distribution models of Darevskia spp. Figure S2: Evaluation metrics for D. dahli, D. portschinskii иD. mixta resulting from MaxEnt models made across a range of feature-class combinations and regularization multipliers. AICc is the Akaike Information Criterion corrected for small sample sizes, delta AICc (DAIC) is the difference between the AICc of a given model and the AICc of the model with the lowest AICc. Dotted horizontal line represents delta AICc = 2, which delimits models that are generally considered to have substantial support relative to others examined -that is those below the line. Default settings and settings that yielded minimum AICc are indicated with arrows. Legends denote feature classes allowed (L = linear, Q = quadratic, H = hinge, p = product and T = threshold). Note that for these lizards, AICc consistently selected regularization multipliers higher than the default value. Figure S3: Relationships between each of the most important environmental predictors (see Table 2) and the likelihood of D. dahli occurrence. Figure S4: Relationships between each of the most important environmental predictors (see Table 2) and the likelihood of D. portschinskii occurrence. Vegetation type: 1-Mountain forest zone, 2-Mountain meadows, 3-Mountain steppe, 4-Arid mountain steppe, 5-Nival zone, 6-Semi-desert, 7-Cultivated areas, 8-Wetland areas, 9-Alpine zone, 10-Urban areas. Figure S5: Relationships between each of the most important environmental predictors (see Table 2) and the likelihood of D. mixta occurrence. Vegetation types are presented in Figure S4. Figure S6: Correlation between predictor variables and the first two components of the principal component analysis calibrated on the environmental conditions in parental and "daughter" lizards. First and second components explain 81.1% of the total variation. The abbreviations of variables are described in Table S1. Figure S7: Graphic representation of the shift of the niche centroid of the parthenogenetic lizard D. dahli relative to the "maternal" species D. mixta along the most important environmental gradients. The red arrow indicates the direction of niches shift. Vegetation types are presented in Figure S4. Figure S8: Graphic representation of the shift of the niche centroid of the parthenogenetic lizard D. dahli relative to the "paternal" species D. portschinskii along the most important environmental gradients. The red arrow indicates the direction of niches shift. Vegetation types are presented in Figure S4. Figure S9: Comparing the means positions (centroids) of the ecological niches of the studied rock lizards along the gradients of the environment using Post Hoc Tukey HSD test. Differences between each pair of means are presented with 95% family-wise confidence level. The names of the predictor variables (see Table S1) for the panels (a)-(j) are the same as in Figure 3. Table S1: Habitat variables considered in the species distribution models. Table S2: Average Nearest Neighbor Index (ANNI) of the species occurrence data, where n is the number of sampling sites, Z score is the statistic value showing validity of the null hypothesis of a random distribution of points.