Genetic Variability and the Potential Range of Darevskia rostombekowi in Transcaucasia

The results of the analysis of the genetic variability of parthenogenetic Darevskia rostombekowi (Darevsky, 1957) species using four microsatellite-containing loci are presented. Based on 118 records with geographical coordinates of the presence of this species in Transcaucasia, the maps of potential range were created. The analysis of the genetic structure of populations demonstrated that despite the established multiclonality (seven clonal lines in four populations), D. rostombekowi was formed as a result of a single act of hybridization between closely related bisexual species. The predicted distribution of D. rostombekowi using the modelling of potential range revealed new suitable habitats, where the presence of the species has not been reported previously. The results of this study and the absence of multiple acts of hybridization during the formation of these clones may indicate a regression of population size of the species. Consequently, the estimation of the conservation status of this parthenogenetic species seems to be justified.


INTRODUCTION
Over the past 50 years, the results of cytogenetic, molecular genetic analyses of unisexual and bisexual complexes of the species of the genus Darevskia were obtained (Darevsky, 1967;Uzzell and Darevsky, 1975;Martirosyan et al., 2002;Malysheva et al., 2002;Malysheva et al., 2006;Osipov et al., 2016;Ryskov et al., 2017). These results indicate the concept of hybridogenic speciation in parthenogenetic species of rock lizards, which consists in the formation of new unisexual hybrid species carrying the genomes inherited from two closely related parental species (Uzzell and Darevsky, 1975;Ryskov et al., 2017). The complex analysis demonstrated that the parthenogenetic D. rostombekowi species was formed as a result of natural hybridization between closely related bisexual species D. portschinskii (paternal species) and D. raddei raddei (maternal species) (Uzzell and Darevsky, 1975). Despite some success in the study of D. rostombekowi (MacCulloch et al., 1997;Martirosyan et al., 2002;Osipov et al., 2016;Ryskov et al., 2017), questions about the boundaries of its range and speciesspecific parameters of abiotic environmental factors (determining the area of distribution of isolated populations) remain unclear. The analysis of allozyme loci in individuals of different D. rostombekowi populations (except for the isolated Tsovak population) revealed no variability, which presumably indicated a monoclonal genetic structure (MacCulloch et al., 1997).
There are fragmentary data on the distribution of isolated populations of D. rostombekowi in Armenia. It is known that the D. rostombekowi lizard (one of seven parthenogeneticspecies of Darevskia rock lizards) occupies a relatively small range consisting of several isolated populations of different sizes within Northern Armenia, adjacent regions of Northwest Azerbaijan, and a small alpine relict (~12000 years) population (detached from the main range) on the southeastern coast of Lake Sevan (MacCulloch et al., 1997;Martirosyan et al., 2002;Arakelyan et al., 2011, Petrosyan et al., 2020a. The boundaries of the range of the species , as well as of isolated populations, are absent in the literature. The aim of this study was (1) to determine the genetic structure of the populations inhabiting different environmental conditions on the territory of Armenia; (2) to construct a map of the spatial distribution of D. rostombekowi based on species occurrence records; and (3) to identify bioclimatic, topographic, and landscape variables that determine the potential range of the species.

Molecular Genetic Analysis
DNA samples of parthenogenetic D. rostombekowi species (n = 42) from four populations on the territory of Armenia were used for this study (Table 1). The lizards were caught in the period from 1997 to 2006 in their natural habitats. All field studies were carried out until D. rostombekowi species was included in the Red Book of Armenia as an "endangered" species according to IUCN criteria (Endangered, EN B1ab) (Agasyan et al., 2010). In this work, the analysis was carried out on the basis of three previously described microsatellite-containing loci Du215, Du281, and Du323 (Ryskov et al., 2017) and the new Du47G locus. These loci were established to study other parthenogenetic species of this genus: D. unisexualis (Badaeva et al., 2008) and D. dahli (Vergun et al., 2014).
A locus-specific polymerase chain reaction was performed using previously selected primer pairs and amplification conditions (Korchagin et al., 2007;Vergun et al., 2014).
The analysis and comparison of the sequences obtained were carried out in the MEGA v.7.0.21 software. The number of alleles, allelic richness as a measure of the number of alleles (adjusted for the sample size), and expected heterozygosity as a measure of gene diversity were calculated for each locus and population using Fstat v.2.9.3.2., GenePop v.4.2, and the online version of POPTREEW. AMOVA analysis was performed using the Poppr package (Kamvar et al., 2014). Realization of multiple alignments of sequences of microsatellite alleles for three species was carried out using a MUSCLE algorithm (Edgar, R.C. et al., 2004). These matrices were used to construct a network by means of the nearest neighbor method using MEGA v.7.0.21 (Tamura et al., 2013). The composi-tion of genotypes was carried out by comparing and analyzing the combinations of allelic variants in the populations using the LabConverter V1 program (Omelchenko et al., 2016).

Creating a Map of the Potential Range of Species
Collecting field data. Data was collected by the authors during field studies in the period 1997-2019 and collection materials of the Department of Zoology of Yerevan State University (YSU), the Zoological Museum of Moscow University (ZMMU), the Museum of the Zoological Institute of Russian Academy of Sciences (ZISP), and the Canadian Royal Ontario Museum (ROM) were used. The potential habitats on routes with a total length of 3000 km in 2018 and 1800 km in 2019 were studied. Field surveys were carried out in the sunny windless hours in the morning. For each locality, the geographic coordinates and heights were determined using a Garmin Montana 680t GPS navigator (Garmin Corp., Olathe, KS, United States) with an error of ±3.5 m. The lizards from each location were captured and identified in the field conditions using a guide (Darevsky, 1967;Arakelyan et al., 2011).
Raster layers for constructing potential range maps. A review of the literature was conducted to select the most important variables determining the distribution of rock lizards (Darevsky, 1967;Uzzell and Darevsky, 1975;Kaliontzopoulou et al., 2008;Tarkhnishvili et al., 2010;Arakelyan et al., 2011;Doronin, 2015;Freitas et al., 2016, Petrosyan et al., 2019, 2020b. The selected environmental variables characterize the climate, topography, and land cover/land use. The spatial climatic variables were taken from a WorldClim 2.0 data set (http://worldclim.org/version2) with a resolution 30 s including 21 bioclimatic variables. In this work, the analysis was carried out for three different resolutions (3, 9, and 30 arc second, that is, ~90, 270, and ~900 m, respectively) of the raster relief layers, climatic parameters, and landscape. All these layers were prepared using the appropriate Arc GIS Desktop 10.4.1 software module using a cubic function of approximation. As a result of these transformations, all raster layers were presented with a resolution of ~90, ~270, and ~900 m, respectively.
To verify the spatial autocorrelation of the occurrence records from VDB, a two-step procedure described in the works (ESRI, 2017;Aiello-Lammens et al., 2015;Petrosyan et al., 2020b) was used.
The potential distribution model of D. rostombekowi was constructed by the maximum entropy method using MaxEnt 3.4.1 (Phillips et al., 2006) (Di Cola et al., 2017). The average Boyce index was determined by averaging the indices of ten models for each resolution separately. The importance of each predictor variable was determined using a table of contribution (AVC) MaxEnt (Phillips et al., 2006). Those variables that make a significant contribution to the model (that is, have high values of permutation importance (PI) or high values of the percentage of contribution (PC)) were considered as important variables (Phillips et al., 2006). Similarity of ecological niches. A comparative analysis of the ecological niches of D. rostombekowi populations obtained for different data sets was carried out using a general concept presented in the following works (Petitpierre et al., 2012;Di Cola et al., 2017). The estimation of the niche overlap obtained based on different sets of occurrence records was conducted using Schoener's D index (Di Cola et al., 2017).
A comparative analysis of the centroids (average values of variables) of niches determining the spatial distribution of D. rostombekowi populations was carried out using the general linear model (GLM) procedure. In this model, the analysis of variance (ANOVA) was used. If ANOVA demonstrated a significant difference in the centroids, then it was determined using the Post hoc Tukey HSD multiple comparison methods precisely which centroids were different from each other.
Statistical analysis and construction of diagrams were performed in the RStudio 1.1.463 software.

Genetic Characteristics of Populations
The analysis of four microsatellite-containing loci Du215, Du281, Du323, and Du47G detected that all D. rostombekowi individuals studied are heterozygous for these loci. From two to five alleles (depending on the locus) were detected in the populations studied. The allelic variants of each locus are presented in Table 2.
The number of genotypes was established by the analysis of the combination of allelic variants of each individual (Table 3, Fig. 1). The analysis of allelic combinations for the four loci detected seven genotypes unevenly distributed in the studied populations. Individuals with identical genotypes formed clonal lines. The major clone (genotype 1) was established in 24 individuals (57% of all individuals studied) and presented in three populations, including Gosh, Papanino (vicinity of Dilijan town), and Spitak. The clones with the genotype 2 were found in seven individuals (17%). The clones with the genotypes 3 (five individuals, 12%) and 4 (three individuals, 7%) were found only in the isolated Sevan population of Tsovak. Other genotypes (5, 6, and 7) were represented by only one individual and were found in the populations Gosh, Papanino, and Spitak.
Genotypes 1-7 differ from each other in the microsatellite cluster sequence. Rare genotypes found in single individuals in the populations are probably mutational derivatives of the main (represented by a large number of individuals and widespread in three out of four populations studied) genotype. The genotypes 3 and 4 are present only in the isolated population of Tsovak and are the «most distant» from the main genotype. The analysis of spatial-frequency distributions of most widespread genotypes 1 and 2 and population-specific genotypes 3 and 4 identified the dependence between the frequencies of clones and the geographical distribution between three northern populations and the Sevan population of Tsovak.
Taking genotype 1 (G1) as the initial, according to Parker's model (Parker et al., 1989), we can assume that the remaining genotypes detected could have formed as a result of microsatellite mutations of the already formed initial clone. This event is reflected in the scheme (Fig. 2), which was constructed based on comparative analysis of the combinations of allelic variants for each locus taking into account the genotype distribution and its frequency in populations. As is seen from the scheme, the genotypes G2, G7 differ from the genotype G1 by a single GATA mutation in the microsatellite cluster by one out of four loci (G1 differs from G7 by a single mutation in the Du215 locus; G1 from G2, by a single mutation in the Du281 locus). The same pattern is observed between the genotypes G2 and G5. Although the genotype G6 is presented in the population Gosh, it is similar to the genotype G2 (present in the populations Spitak and Papanino), except for deletion of one GCAA element in the Du215 locus. Despite the fact that no genotypes G3 and G4 were found in other populations (except for Tsovak), they are associated with the main G1 genotype by common allelic variants; however, they are "distant" from it by multiple changes within the microsatellite cluster.
Data on the calculations of the genetic diversity indices in the populations of parthenogenetic species BIOLOGY BULLETIN Vol. 48 No. 6 2021 OSIPOV et al.

Most Important Predictor Environmental Variables and the Potential Range of Parthenogenetic Species
We collected 118 records with geographical coordinates of D. rostombekowi occurrence data in Armenia and Azerbaijan (49 published accounts and 69 of our own field records). After removing cluster points, there are 56 records left included in the reduced set for further analysis. Subsequently, these records (full n = 118, reduced n = 56) were used to construct species distribution models and to compare ecological niches for three resolutions ~90, ~270, and ~900 m, respectively.
We obtained "high" values of the Boyce index (B ind ±SD) and AUC (±SD) for all species distribution models. The contributions of variables for different SDM are presented in Table 5. It follows from Table 5 that the AUC index (Phillips et al., 2006) has sufficiently high values and varies from 0.978 (±0.007) to 0.985 (±0.007), which are slightly different from each other. The latter means that, regardless of the type of data use (reduced, full) and the spatial resolution of the layers, the AUC index of model suitability changes slightly (Table 5). Unlike the AUC, the Boyce index is sensitive and allows us to select the best model. In our case, the highest accuracy of the species distribution model is observed only for a single case, that is, when full data with the spatial resolution ~90 m are used. In this case, the Boyce index is 0.958 ± 0.004. The worst model according to the Boyce index 0.813 (±0.02) was constructed when using reduced data with the spatial resolution of ~270 m. The maps of the species distribution for these cases are presented in Fig. 3. Du215(2 + 5) + Du281(2 + 4) + Du323(1 + 2) + Du47G(4 + 5) 0 0 6 1 7 3 Du215(4 + 5) + Du281(3 + 4) + Du323(1 + 2) + Du47G(2 + 1) 0 5 0 0 5 4 Du215(4 + 5) + Du281(3 + 4) + Du323(1 + 2) + Du47G(4 + 5) 0 3 0 0 3 5 Du215(2 + 5) + Du281(2 + 4) + Du323(1 + 2) + Du47G (  It follows from Table 5 that only five out of eight predictor variables (MeanTempDrQ, Sradmean, Pre-cipCoefVar, IsoTerm, PrecWarmQ) are common for all 60 models (ten replications × two types of data sets × three resolutions). The importance of the StdDiv-Temp variable is highlighted for almost all models except for the models obtained for reduced data for the resolution of ~90 m. The importance of the HighWay variable is highlighted for all models with different resolutions when using the full set of find point data, as well as for the models with the resolution of ~90 m when using the reduced data. Elevation is the only variable that is highlighted only for a single combination of the data set (full) and resolution (~90 m).
A visual qualitative analysis of the best (Fig. 3a) and the worst (Fig. 3b) maps demonstrates that they coincide in general terms along the boundaries of the range, but there is a great difference in the details. The best map (Fig. 3a) indicates that the most suitable habitats (with a probability >0.6) are mainly located in northern Armenia, including all known habitats of the populations (Sevan, Spitak, Gosh, Dilijan, Ijevan-Noyemberyan), as well as new habitats of the western part of Lake Sevan and in the eastern part of Nagorno-Karabakh, the northern part of Armenia on the border with Georgia. The worst map demonstrates extensive suitable habitats on the western part of Lake Sevan, a vast area in the eastern part of Nagorno-Karabakh, in northern Armenia at the border with Georgia, and plots in the southern part of Georgia.

Similarity Indices of Ecological Niches of Populations of Parthenogenetic Species
The results of the analysis of niche overlapping for different data sets are presented in Fig. 4. It was demonstrated that the similarity index of ecological niches for different resolutions changes in the range from 86 to 90% depending on the resolution of the layers used. The portion of common habitats (that is, niche overlapping) varies from 97.8 to 98.8%. The similarity test of ecological niches obtained based on different data sets does not deviate for all resolutions (P < 0.02).
Despite the great similarity of ecological niches when using different data sets, however, the centroid is shifted upward in height for reduced data (Figs. 4d-4f). This means that thinning in order to remove the clustered points leads to a significant decrease in the number of points in the habitats of northern populations (Spitak, Gosh, Papanino, and Ijevan-Noyemberyan) and leaves unchanged the number of find points in the southern part of Lake Sevan (Tsovak). Eventually, this leads to an increase in the centroid height and an increase in the habitat suitability at an altitude of more than 2000 m.   TA  CG  TA  CG  TA  CG  TA  CG  TA  CG  TA  CG  TA  CG   TA  CG   T  C   TTC  AAG  TTC  AAG  TTC  AAG  TTC  AAG  TTC  AAG  TTC  AAG  TTC  AAG   T

Comparative Analysis of Centroids of Niches of Parthenogenetic Species Populations
To check the statistical significance of the difference of niche centroids by abiotic environmental factors of different populations, we carried out their comparative analysis. The results of comparison of the mean values of predictor variables (niche centroids) of united northern populations and the isolated Tsovak population are presented in Fig. 5. It follows from Fig. 5 that the ecological conditions of the Tsovak population are statistically significantly different from northern populations in the average temperature and total precipitation of the dry quarter of the year, the average solar radiation, precipitation coefficient of variation, isothermality, and elevation. The southern population differs from the northern populations not only by the distance to the roads, but also by standard deviation of temperature.

DISCUSSION
A comparative analysis of single nucleotide variants outside the microsatellite cluster (SNV) of D. rostombekowi parthenogenetic species with parental spe-  (Osipov et al., 2016;Ryskov et al., 2017) suggests that the parthenogenetic species could be formed as a result of a single hybridization event between close related bisexual individuals of the parental species, since the studied individuals of the parthenogenetic species did not differ by SNV in all four microsatellite loci (Fig. 1).
The analysis of distributions of the most widespread genotypes 1 and 2 and population-specific genotypes 3 and 4 revealed a dependence between the fre-quencies of clones and geographical separation between three northern populations and Sevan population Tsovak. The question remains which of the seven D. rostombekowi clones established was formed first. Assuming that the parthenogenetic species formed as a result of a single hybridization event has a spread and most common "major" clone with several "minor" clones (Parker et al., 1989), the genotype G1 could probably be the initial one. The remaining genotypes are G1 mutations. Mutations de novo in all four loci studied are caused by deletion or insertion of a 1-reduced set of occurrence points 2-full set of occurrence points 1 2 microsatellite repeat (Fig. 1, Table 2). Noticeable differences of the genotypes G3 and G4 from the main genotype G1 and genotypes G2, G5, G6, G7 (found in the populations Gosh, Papanino, and Spitak) can be associated with significant differences in the ecological conditions of their habitats (Fig. 5).
Most likely, other (transitional clonal lines) between genotypes 1-3 and 4 could have been lost as a result of reduction and/or rupture of the range of D. rostombekowi.
In this work, it was demonstrated that SDM is an important instrument, which can be used to detect ecological requirements of D. rostombekowi populations to the environment and to select the most important variables determining their spatial distribution. The additional analysis using an ordination method allows us to conduct a comparative analysis of ecological niches of different populations and to determine the measure of their overlap and similarity and to detect significant shifts of niches in identifying  potential ranges on the basis of different sets of occurrence records and resolutions of the layers of predictor variables. These two approaches, as well as the genetic structure of the populations, allow us to reveal their peculiarities and to estimate the role of ecological factors in the formation of the northern and southern populations of this species.

Most Important Variables Determining the Distribution Model of Parthenogenetic Species
As demonstrated by SDM, the average temperature and total precipitation in the dry and warm quarters of the year, the solar radiation, precipitation coefficient of variation, isothermality, and elevation are important environmental factors that affect the formation of the spatial distribution of lizards on the territory of Armenia, Nagorno-Karabakh, and Azerbaijan (Fig. 3, Table 5). It follows from this that the set of variables determining suitable habitats is largely associated with the thermal conditions and humidity in the dry and warm seasons, respectively, as well as with the annual coefficients of variation of these variables. Moreover, isothermality should be within 28-34%, while the precipitation humidity coefficient is 41-52% (Fig. 4). It is important to note that the isothermality, which reflects sufficiently high temperature stability and its low fluctuations in habitats, has not previously been noted for other lizard species as the most important predictor variable (Kaliontzopoulou et al., 2008;De Pous et al., 2011). The narrow range of variation of these variables means a sufficiently high sensitivity of the species to climate change. The established vari-ables explain the ecological peculiarities associated with the distribution of D. rostombekowi and give an idea of the factors limiting its spread. Solar radiation is another important variable, which provides warming up of the soil in habitats and changes in the narrow range of 4.1-4.2 kW*h/m 2 . The elevation is the important variable only for one (the best) model, which was obtained using layers with the resolution of 90 m (3 s) (Table 5).

Potential Range of Parthenogenetic Species
The occurrence records and potential range of D. rostombekowi demonstrate that the species is distributed quite widely in the northeastern part of Armenia, in the northwestern part of Azerbaijan, and on the territory of Nagorno-Karabakh (Fig. 3a). Most places where D. rostombekowi was previously observed are easily predictable using SDM (Darevsky, 1967;Uzzell, 1975;Martirosyan et al., 2003;Arakelyan et al., 2011). Actually, the model of the potential range of D. rostombekowi predicted a wider range of distribution outside of the already known suitable habitats of the populations Spitak, Gosh, Papanino, Ijevan, Noyemberyan, and Tavush, as well as the southern population on the coast of Lake Sevan, where the species was already registered. In addition, the model highlights new small suitable habitats of the western part of Lake Sevan and the eastern part of Nagorno-Karabakh, and the northern part of Armenia on the border with Georgia. In general, the low demand of the species for vegetation according to the results of MaxEnt modelling is reflected in the potential distribution model, Table 5. Sets of predictor variables determining the models of D. rostombekowi distribution with indication of the indices of suitability (AUC and Boyce) for maps of potential ranges obtained using MaxEnt and EcoSpat. Only variables that contributed more than 5% either by PC or by PI are presented * MeanTempDrQ, average temperature of the dry quarter of the year (°С); HighWay, distance to roads (m); Sradmean, average solar radiation (Kj m -2 day -1 ); PrecipCoefVar, precipitation coefficient of variation (%); IsoTerm, isothermality (%);PrecWarmQ, amount of precipitation in the warmest quarter of the year (mm); StdDivTemp, standard deviation of temperatures (°С); Elevation, height above sea level (m). which includes regions that are located in the forest zone (64%), mountain meadows (12%), mountain steppes (14%), and anthropogenically transformed habitats (9%) at the heights typical for mountain forests (1394 ± 37 m) and mountain meadows and steppes (1957 ± 146 m). Northern populations are mainly located in the forest zones; the isolated Sevan population is found in mountain meadows and steppes.
CONCLUSION This study systematized knowledge about genetic variability and the distribution of parthenogenetic D. rostombekowi species in the native part of the range using four microsatellite-containing loci Du215, Du281, Du323, and Du47G and 118 records with the geographical coordinates of the species presence from the territory of Armenia, Azerbaijan, and Nagorno-Karabakh covering most of its range. During this study, we supplemented the available literature data with our own field observations and constructed the potential range of the species (PRS) based on the analysis of climatic, topographic, and landscape environment variables. The results of the genetic structure (GS) of the populations demonstrated that, despite the presence of seven clonal lines, the parthenogenetic species is characterized by a single monophyletic form; that is, it is formed as a result of a single hybridization event between close bisexual individuals of parental species. The analysis of PRS and GS indicated the difference of an isolated population of the southeastern part of Lake Sevan and northern populations in Armenia. These differences are expressed both in the shift of ecological niches of these populations relative to each other and in the difference of the established genotypes. The currently predicted distribution of D. rostombekowi using PRS identified new small habitats in the western part of Lake Sevan and the eastern part of Nagorno-Karabakh, and the northern part of Armenia on the border with Georgia, where the presence of the species has not previously been reported in the literature. The results of this study and the absence of polyphyletic forms can indicate its regression. The results are valuable to support future studies and can be used to guide environmental organizations in order to preserve the most important habitats of the species studied.