Fine-scale analysis of habitat occupancy by Kura lizard (Darevskia portschinskii) and its daughter parthenogenetic form (Darevskia dahli)

Two species of rock lizards, the parthenogenetic D. dahli and the sexually reproducing D. portschinskii, coexist in a rocky outcrop in an area of ca. 1 ha, in the vicinity of Tbilisi, Georgia; the location has been well-known since the middle 1960s. The population density of the parthenogenetic lizard is five times higher than that of the sexual breeder. We studied the distribution of active lizards in space and time over three consecutive years, during the spring and autumn activity periods, to explore spatial and temporal differences between the species on a fine spatial scale. We studied the influence of temperature, humidity, and quantitative characteristics of the surface and the distance from permanent water source on the spatial distribution of D. dahli and D. portschinskii. Darevskia portschinskii was less dependent on the distance from the water source and more evenly distributed in space and time than D. dahli. Despite potential competitive interactions, the species did not avoid each other on the microhabitat scale, suggesting that the observed ecological differences are not caused by a niche shift. More individuals of the sexual breeder than individuals of the parthenogen were found in suboptimal habitats. This feature may increase the evolutionary success of D. portschinskii in a long-term perspective.


Introduction
Comparing ecological niches of closely related species is a challenging task. A species' spatial distribution and activity dynamics depend on its ecological preferences. Simultaneously, the realized niche in Hutchinson's (1957) sense is shaped by competitive interactions, both interand intraspecific. In the case of the presence of a competitor, the space occupied by a population can be shifted, as well as the activity dynamics (Schoener 1983;Abrams 1986;Pfenning and Pfenning 1992;Pianka 2011;Carvalho and Cardoso 2020). This process affects empirically estimated tolerance limits and response curves (Huey and Stevenson 1979;Kassahn et al. 2009) can be different, dependent on the species composition in a biological community (Petren et al. 1993;Pianka 2011;Frishkoff et al. 2015). In general, competition is a significant force causing character displacement and, hence, organismal evolution (Grant 1972;Pfenning and Pfenning 2010).
Parthenogenetic lizards have hybrid origins, and they commonly coexist with their ancestral bisexual forms (Darevsky 1967;Bakradze 1977;Trofimov 1981;Wright and Vitt 1994;Kearney et al. 2009;Tarkhnishvili et al. 2010;Tarkhnishvili 2012;Galoyan 2013). Coexisting parthenogens and their ancestral species are phenotypically similar, they consume similar food, and most of them occupy rocky habitats (Darevsky 1967;Wright and Vitt 1994;Tarkhnishvili 2012). The offspring of sexually breeding males and parthenogenetic females are usually sterile (Freitas et al. 2019;Tarkhnishvili et al. 2020). Hence, the parthenogens and their ancestral species can potentially compete for space or food; niche shift and character displacement may be an important force in their phenotypic evolution.
One of such hybrid parthenogenetic lizards is Darevskia dahli from the Lesser Caucasus Mountains, the patrilineal ancestor of which is the bisexual Darevskia portschinskii and the matrilineal ancestor is Darevskia mixta (Murphy et al. 2000). The extent of occurrence of D. dahli and D. portschinskii almost coincides, and their habitats overlap: D. portschinskii is found at altitudes 400-1,370 m. a.s.l., and D. dahli at higher altitudes between 830-1,950 m. a.s.l. (Tarkhnishvili et al. 2010). A higher altitudinal range is typical for unisexual species compared to their bisexual relatives (Haag and Ebert 2004) and may help avoid competition (Glesener and Tilman 1978;Kearney et al. 2005;Vrijenhoek and Parker 2009). However, the species coexist throughout a broad mountain forest belt and are commonly found in the same location (Bakradze 1977).
Remote sensing analysis showed that the shift of the spatial distribution of the two species is related to different thermal preferences and to different preferred humidity: D. portschinskii prefers higher temperatures than D. dahli, and the distribution of D. dahli is limited to the areas where annual rainfall level is above a certain threshold (Tarkhnishvili et al. 2010;Petrosyan et al. 2020). Simultaneously, Tarkhnishvili et al. (2010) showed that, within equally suitable habitats, the presence of D. dahli might negatively affect D. portschinskii: the latter species was less abundant or absent if D. dahli was also present in equally suitable locations (Tarkhnishvili et al. 2010), a fact suggesting the importance of interspecific competition. It is unclear whether competition may explain niche shift and character displacement in D. dahli and D. portschinskii.
For deeper insight into this question, it is essential to understand how the species share space and other resources on the microgeographic scale and whether competition for the resources or behavioral interactions shape the realized niche of either species. Here, several interacting processes may play a role. Galoyan (2011Galoyan ( , 2013 showed that both males and females of some sexually reproducing Darevskia species may be territorial, whereas the parthenogens show decreased territoriality. One should expect that territorial species do not concentrate so much in areas with optimal living conditions as the species with less expressed territoriality (Brown 1969;Fretwell 1969;Partridge 1978; but see Stamps 1991); this may potentially reduce the mutual negative impact of D. portschinskii and D. dahli. Simultaneously, populations of parthenogenetic species grow faster than those of sexual breeders (Maynard-Smith 1998) and hence may rapidly expand into sink habitats (sensu Pulliam 1988); as a result, the differences in the observed niche dimensions may depend on this spatial dynamics and not on the different ecological preferences.
Here, we examined whether ecological niches of D. dahli and D. portschinskii are displaced on a microgeo-graphic scale (within a semi-isolated habitat covering ca. 1 hectare) and what can be the reason for the differences (if any): type of territorial behavior, population density, or competition. We tested two complementary hypotheses: (1) the differences in the habitat use are related to mutual avoidance of the coexisting species, and/or (2) the differences are due to differential territorial behavior; alternatively, habitat requirements might not differ or differ for reasons independent of the intra-or interspecific competition.

Study area
A rock outcrop near Tbilisi (Georgia) was selected as a study site. This area is described in previous papers (Darevsky 1967;Bakradze 1977;Trofimov 1981;Tarkhnishvili et al. 2010) and presents a good model for longterm stationary observations of the lizards. The area is located at the vicinity of Kojori (1.648939°N, 44.683445°E) (Fig. 1). It is a gentle rocky slope covered with a rare shrub, surrounded by a forest dominated by oriental hornbeam (Carpinus orientalis), rocky oak (Quercus petraea), and maples (Acer campestre, A. letum) and edged with a small brook from the north. The slope has northwestern exposure, and the elevations vary between 1003-1042 m a.s.l. (Fig. 1). The observations were conducted between April and October of 2017, 2018, and 2019; the full data are presented in Suppl. material 1: Table S1.

Description of environmental variables
We divided the study area into 108 squared plots with a surface of 10 × 10 m each ( Fig. 2a; Suppl. material 1: Table S1). For each square, longitude, latitude, and altitude were measured using Garmin GPS (Geographic Positioning System) MAP 64 and then adjusted manually. Geographic Positioning System data were processed in the QGIS 3.12 software and turned into a polygon shapefile (Fig. 2b). DJI Mavic 2 Pro camera drone was used for taking aerial photos of the surface (distance from the ground 60 m). Georeferencing of aerial images was performed using Stonex S9III+ GPS (https://www.stonex.it/project/ s9iii/). Aerial photos were then processed in the Agisoft PhotoScan software to create Orthoimage and DEM (Digital Elevation Model; Agisoft, L. L. C. "Agisoft Pho-toScan User Manual: Professional Edition, Version 1.4", 2018). We used orthoimage for mapping surface types (Boulders, Forest, Grass, Gravel, Leafy Ground, Rock) manually (Fig. 2c). Digital Elevation Model was resampled in two more spatial resolutions (original resolution 0.03 m/pixel and resampled resolution 0.5 m/pixel; 1 m/ pixel). We used SAGA's surface morphometry tool in QGIS 3.12 for estimating slope and convergence index from the three mentioned resolutions of DEM.
Standard deviation, average and maximal values of slope and convergence index in each sample plot were calculated using the Zonal Statistics tool in QGIS 3.12. The slope represents the inclination angle to the horizontal surface in degrees; the maximum value of slope (MSlope), characterizing each plot, was used in final calculations. The convergence index calculates an index of convergence/divergence regarding overland flow (2004-2016. Our study used the standard deviation of convergence index (SDCI), which measures the unevenness of the rocky surface or density of clefts that the lizards are using as hiding shelters.
In each field session (27 field sessions altogether -Suppl. material 1: Table S1), for each of the 108 plots, the land surface temperature (LST, °C) and land surface humidity (LSH, %) was recorded twice -between 1100AM and 1300 PM and between 1700 and 1900 PM. The measurements were taken in the geometric center of each square plot. The temperature was measured with General Tool MMD7NP; the humidity was measured with Etekcity Lasergrip 774. Ambient temperature from https:// www.accuweather.com/ was recorded during each field session. Altogether 2916 observations (field survey/ plot) were conducted during the study period. During each observation, 0-10 individual lizards were recorded.
For each sighted individual, we recorded: age (adult or juvenile); species (D. dahli vs D. portschinskii); sex (for D. portschinskii); time of the day and the number of the plot. All individuals with the distance from the back surface of the hip to the tip of snout exceeding 5 cm were treated as adults.

Statistical procedures: analysis of habitat use in space and time
The numbers of active D. dahli and D. portschinskii were recorded for each of 108 plots 100 m 2 each, during each of 27 field sessions. Hence, altogether we had 2,916 observations for each species. Surface temperature and humidity were scored for most of these observations. The entire observation table is shown in the appendix (Suppl. material 1: Table S1). In 18 plots, the lizards were never recorded, and in 35 plots, fewer than three lizards of either species were recorded during the entire period of observations. In the rest of the plots (55 altogether), the number of observed D. dahli during one session varied between 0-11, and the number of observed D. portschinskii between 0-3.
Because males of D. portschinskii might have spatial distribution different from females and juveniles of the same species, due to potentially more expressed territorial behavior, we first tested significance in differences in distribution between D. dahli, males of D. portschinskii, and females + juveniles of D. portschinskii across the 108 plots, using Kolmogorov-Smirnov Lambda test. Since the differences were significant for different species but insignificant for different ages and sexes of D. portschinskii, we clumped together all individuals of the latter species for further analyses.
In order to evaluate the significance of each of the factors that potentially influence the presence and activity of the lizards (the number of recorded individuals at a plot during a given session; NL; the analysis separately ran for D. dahli and a unified sample of D. portschinskii) we applied generalized linear mixed models (GLMM). Because of the excess of NL equal to zero, before applying this method, we removed from the dataset all plots where fewer than three lizards were recorded during the entire period of the study, reducing the number of observations to 1,483. This was done to achieve better correspondence of the NL distribution throughout the plots and sessions to Poisson distribution, recommended for this kind of data. The plots were set as "subjects" and the sessions as "repeated measures." Target (NL) was set to Poisson distribution with the logarithmic link function. Nine predictors (year, date, time of the day, surface type, humidity and temperature of the ground surface, averaged air temperature, slope, and convergence index) were set as fixed factors; the session combined with the plot was set as a random factor. Covariance type was set to "Variance component" for the subject and to "first-order autoregressive" for the random component (assuming autocorrelation between the sessions close in time). Classification tree analysis (CHAID algorithm) was an additional method applied for inferring predictive models for both species. In this analysis, we included the occurrence of a potential competitor (PC) in the set of predictors to infer whether and how the presence of PC may influence the occurrence of the target species, given similar environmental conditions. In this analysis, we used a complete dataset (2,916 observations, including those on the plots where the lizards were never recorded). Parent and Child nodes were set to 100 and 50. Significance lev-els for splitting nodes and merging categories were set to 0.05, growth limits to 2 levels, the maximum number of iterations to 10,000, and Chi-square statistics calculated using likelihood ratio. For validation, the entire dataset was randomly separated into the training (75%) and test (25%) datasets; squared correlation (R2) between the observed mean in the training and test datasets was used as an index of the model quality. MacArthur and Levins (1967) suggested several baselines for the analysis of competition and niche displacement. They offered to use an asymmetric coefficient of

Estimation of niche overlap
In this equation, O 21 is the number of individuals of species 2 whose resources are consumed by a single individual of species 1; p 1i is the index of the utilization of resource category i by species 1. We considered a single observation as a "resource category" and estimated p xi as NL for species x in this observation. Before that, we conducted this procedure separately for individual plots and survey sessions and estimated the "overall" overlap indices as the product of the two latter ones (Schoener 1974;Colwell and Futuyma 1971;Hurlbert 1978). Furthermore, we used the symmetric Pianka's index of niche overlap (Hurlbert 1978).
All calculations were done using IBM SPSS v. 23.0 and Excel 2010 for Windows.

Response on the individual environmental variables
Darevskia dahli were recorded at 78 out of 108 studied plots, and D. portschinskii at 40 plots (Fig. 3). The average number of the total records per one plot was: 5.62 (SD = 9.07) for D. dahli, 0.44 (SD = 1.07) for males of D. portschinskii, and 0.43 (SD = 0.90) for females and juveniles of D. portschinskii. Kolmogorov-Smirnov statistic D n,m reached 0.382 and 0.345 (P < 0.001) when comparing spatial distribution of D. dahli with males and with females + juveniles of D. portschinskii, respectively. Simultaneously, D n,m reached 0.089 (P > 0.2) when comparing spatial distribution of males with that of females + juveniles of D. portschinskii. In conclusion, males, females and juveniles of D. portschinskii do not differ significantly in the use of space, and simultaneously the differences between the studied species are highly significant.
Both Pearson correlation and Spearman rank correlation coefficients showed a significant negative correlation of the number of records of D. dahli with the distance from the brook, surface humidity, and convergence index, and positive correlation with an average slope of the plot (Suppl. material 1: Table S1). For D. portschinskii, correlation with the distance from the brook was positive (Fig. 4a), whereas correlation with surface humidity was negative. The abundance of D. dahli (but not that of D. portschinskii) decreased with the distance from the brook. In the part of the study area most distant from the brook, the abundance of both species increased, however (Fig. 4a). Both species were more abundant in 2017 than in the two following years, although in D. portschinskii the differences were insignificant (Fig. 4b). NL of D. dahli decreased when the air temperature measured before the afternoon was below 13 and above 21 C; in D. portschinskii the dependence was less clear (Fig. 4c). The peak of activity of both species was in the last days of May (Fig. 4d). Darevskia dahli preferred the plots where the average slope index was around 6 and avoided flat areas (Fig. 4f). Finally, both species depended on the presence of rocky surface (Fig. 4e). The convergence index did not show any trend.
The analysis based on the assumption that the individual plots are "subjects", individual sessions are "repeated measures", and NL has Poisson distribution with Log function, suggested a significant effect of all predictors, except surface temperature and the convergence index on the occurrence of D. dahli, and presence of temporal autocorrelation explaining the occurrence dynamics (Table 1). This species avoided forest and preferred rock as the surface; the occurrence decreased with surface humidity and increased with the slope steepness. The occurrence of D. portschinskii increased on the gravel and rock surface and with distance from the brook. None of the other predictors significantly affected the abundance or activity of this species (Table 1).

Classification tree
R2 value estimated from the node means of the training and test CHAID trees for D. dahli was 0.859 (17 nodes), and for D. portschinskii 0.575 (11 nodes).
The classification tree showed 3-4 variables that significantly influenced NL of D. dahli (Suppl. material 2: Fig.  S1a). The primary variable separating the observations into five categories was the distance from the brook. Given the same distance from the brook, the NL (D. dahli) increased on the plots where rocky surface dominated, and the average slope was higher (>=65). The occurrence of D. portschinskii was not among significant factors, inferred by the analysis, which influences the occurrence of D. dahli.
The distance from the brook was also a primary variable affecting NL for D. portschinskii. Simultaneously, given the same distance from the brook, NL (D. portschinskii) positively correlated with NL (D. dahli), suggesting that similarity of the requirements to the environmental conditions, except distance from the permanent water source, outweigh potential negative interactions between the two species.

Niche overlap indices
The asymmetric index of niche overlap calculated for the utilization of space and time was: for dp/dd 0.543, and for dd/dp 0.053. If we ignore time and only count the overlap in space use, the respective coefficients would be 2.69 (dp/ dd) and 0.10 (dd/dp). Symmetric Pianka's niche overlap index reached 0.224 for individual observations and 0.487 for individual plots irrespective of the observation time.

Discussion
Competition is a significant force, able to lead to one of the competitor's extinction if the differences between the utilized resources are insufficient (Gauze 1934;MacArthur 1970;Pianka 1976;Şahin et al. 2020), or to niche shift and character displacement, if the adaptive strategy is mutual avoidance of the competing species (Grant 1972;Pianka 1976). The distance between the resource utilization curve peaks should exceed the average standard deviation of the curves to prevent competitive exclusion if resources are limited (MacArthur and Levins 1967;Begon et al. 1996). There is also a simpler definition: competition between the individuals of different species should not exceed competition among the conspecifics (Angilletta et al. 2013). Are the observed ecological differences between coexisting D. portschinskii and its daughter parthenogen, D. dahli, driven by interspecific competition? Our results reject this hypothesis.
Remote sensing analysis (Tarkhnishvili et al. 2010;Petrosyan et al. 2020) suggests that summer temperature and rainfall determine different altitudinal ranges for D. dahli and D. portschinskii. However, at a fine-scale, these species did not react differentially to temperature and humidity of the surface, probably because of the high variability of these factors. Our analysis suggests that the most critical environmental predictors influencing the fine-scale spatial distribution of adult lizards include the distance from the permanent water source, steepness of the slope, and presence of the rocky surface. Darevskia dahli relies more on the physical closeness to the water source than D. portschinskii (also reported in Bakradze 1977). The former species is also more selective to the slope steepness, air temperature, surface humidity, and surface type. The presence of a permanent water source determines the availability of water for lizards, similar to the rainfall level on the macrogeographic scale (Tarkhnishvili et al. 2010;Petrosyan et al. 2020), whereas measured surface humidity quickly changes in time and does not determine the suitability of the location, albeit affecting the activity of the individuals (GLMM results).
Why do environmental conditions, including surface type and characteristics, proximity of water source, and air temperature, affect the occurrence of D. dahli more markedly than that of D. portschinskii? Tsellarius (Tsellarius and Tsellarius 2003;Tsellarius et al. 2016) and Galoyan (2011) showed that in bisexual Darevskia, territorial behavior is an important factor influencing spatial distribution. Adults of D. brauneri from the northwestern Caucasus, both males, and females, defend the individual territory from the conspecific outsiders (Galoyan 2013).
In contrast, the parthenogens do not show territorial behavior, and their spatial distribution depends solely on environmental conditions (Trofimov 1981;Galoyan 2013). We did not conduct behavioral observations on D. portschinskii and cannot conclude whether territorial behavior is typical for females of this species. However, there are no significant differences in spatial distribution between the males and females of D. portschinskii. These facts explain the less aggregated distribution of active D. portschinskii, related to D. dahli, in space and time: the latter is concentrated in space with optimal conditions, whereas the bisexual species utilizes space and time more uniformly. The variance of the distribution of D. portschinskii along the environmental gradients was higher than that of D. dahli. The territoriality of D. portschinskii forces the individuals to expand on plots with suboptimal conditions. Hence territoriality, together with the inherited ecological differences, explains the observed differences between D. portschinskii and its daughter parthenogen. Territoriality has another interesting consequence: the territorial species, which is more commonly found in suboptimal habitats and time, may better adapt to changing conditions, if previously optimal parts of the habitat will be lost. Hence, the advantage of sexual reproduction, along with the well-known factor of genetic diversity (White 1970;Case and Bender 1981;Kondrashov 1982;Maynard Smith 1998;Bell 2008), also favors diversity of ecological conditions in a habitat, which is substantially lower in the parthenogenetic form. Case (1990) suggests that parthenogens are weaker competitors than their sexually reproducing relatives; an additional reason for this can be the absence of territoriality and high concentration in optimal areas. This seemingly contradicts an earlier finding suggesting that the negative effect of D. dahli on D. portschinskii is stronger than vice versa (Tarkhnishvili et al. 2010); however, more effective use of suboptimal habitats may have a long-term positive effect challenging to record in a short-term perspective, similar to the advantage of sexual breeders associated with the ability of recombination.
In spite of these differences, the overlap of spatio-temporal niches between D. dahli and D. portschinskii is substantial. This may potentially cause interspecific competition (MacArthur and Levins 1967). If competition exists, it can be strongly asymmetric; because of the exceeding abundance of D. dahli, the latter species can have a more substantial impact on D. portschinskii than vice versa. The niche overlap by itself cannot either support or reject the presence of competition (Holt 1987). However, one can infer active mutual avoidance of potential competitors by applying residual analysis (Tarkhnishvili et al. 2010) or by analyzing the decision-making tree (this study). This analysis helped us to exclude mutual avoidance of the two studied species. Within the studied area, preferred locations for the lizards partly overlapped. On the plots close to the brook, the density of D. dahli was very high, and D. portschinskii were almost absent. In plots with optimal terrain and cover conditions, the density of both species increased even away from the brook. Within an equally suitable area, D. portschinskii tended to occur simultaneously and in the same plots where D. dahli also presents. Hence, given the same parts of the studied outcrop, microhabitat similarities of the two species are more important than the potential competition, and therefore, ecological differences hardly can be the result of character displacement. Why do the differences occur?
These differences, which are related to the temperature and humidity optima, reported here and in earlier observations and modeling results (Darevsky 1967;Bakradze 1977, Tarkhnishvili et al. 2010Petrosyan et al. 2020) are instead related to the different ranges of D. portschinskii and the matrilineal ancestor of D. dahli, Darevksia mixta. The range of D. mixta is limited to mesophylic forests and uplands of central and western Georgia (Darevsky 1967;Tarkhnishvili 2012;Gabelaia et al. 2015). This species is adapted to more humid habitats than D. portschinskii (Petrosyan et al. 2020). Darevskia dahli has phenotypic features intermediate between D. portschinskii and D. mixta, being more dependent on the water source than the former one and more dependent on the rocky surface than the latter.
Although the species do not avoid each other on the microhabitat scale, on the macrogeographic scale the mutual negative impacts are obvious (Tarkhnishvili et al. 2010).
In certain locations with an equally suitable climate, the abundance of each species negatively correlates with the abundance of the other. The reason for this is not clear. It may be competition for food: rock lizards eat arthropods of appropriate size (Darevsky 1967;Tsellarius and Tsellarius 2001;Galoyan et al. 2019). Like other ectothermic organisms feeding on small invertebrates (e.g., Kuzmin and Tarkhnishvili 1992), they hardly can seriously affect the abundance of insects; rather, a high concentration of animals may decrease the hunting success. Competition for space is probably a more inclusive explanation of the negative interactions: it influences the accessibility on food and shelters. An additional (and perhaps even more important) reason for the negative impact of D. dahli on D. portschinskii includes mating attempts of the males of the former species with the parthenogen. Commonly, reinforcement favors pre-zygotic isolation and, hence, decrease risk of unsuccessful mating with non-conspecific individuals (Dobzhansky 1937;Endler 1986). However, parthenogenetic Darevskia commonly mate with their patrilineal ancestors (Bakradze 1977;Trofimov 1981;our observations). Sterile triploid hybrids are relatively common, and sometimes even have developed gonads (Danielyan et al. 2008). In spite of the latter observation, no ongoing gene flow was ever shown (Darevsky and Kulikova 1962;Spangenberg et al. 2017;Freitas et al. 2019), and even if it could happen, this would be an extremely rare event (Tarkhnishvili et al. 2017(Tarkhnishvili et al. , 2020. Many adult D. dahli have tracks of mating bites. We could hypothesize that common mating of D. portschinskii with parthenogenetic females may decrease their reproductive success. This study shows that negative interspecific interactions do not necessarily trigger niche shift and character displacement and that the simultaneous presence of negative interactions and ecological differences between closely related species cannot be taken as an argument for character displacement. The studied species do not actively avoid the competitor, which suggests that their long-lasting coexistence is a result of metapopulation dynamics and sources-sinks type of spatial structure (Pulliam 1988), where the rates of re-population from neighboring source locations compensate the negative trends in population dynamics. The populations of D. dahli and D. portschinskii have coexisted in the studied location for at least 53 years (Darevsky 1967), and during all this time, the abundance of D. dahli in the studied habitat exceeded several times the abundance of D. portschinskii (Trofimov 1981;Tarkhnishvili et al. 2010;this paper). In conclusion, we suggest that a complex spatial structure not only helps species to survive a competitive pressure but also impedes character displacement and reinforcement in multispecies assemblages.  Fig. 2a). dd -the number of observed active of D. dahli, dp -the number of observed D. portschinskii, Humsurface humidity index, Temp -surface temperature, Air.Temp -air temperature, Con-stdev -convergence index, Slope_max -average slope, Time of Day: Mearly afternoon, N -before sunset, Surface -surface type as shown in Fig. 2c, N -northing, E -easting of the center of a plot.

Supplementary material 2 Figure S1
Authors: Natia Barateli, David Tarkhnishvili, Giorgi Iankoshvili, Giorgi Iankoshvili, Nikoloz Dvali, Zurab Janiashvili Data type: PNG figure, classification tree describing habitat preferences of the studied species Explanation note: Classification tree (=decision tree), describing relation of occurrence of D. dahli (dd) and D. portschinskii (dp) to the most important environmental variables. Method -CHAID; Maximum depth -2 steps; minimum size of parental node -100; minimum size of child node -50; Validation: 75% training observations, 25% test observations. The total number of observation -2917. Test tree is shown. Variables/ predictors: dd -number of observed D. dahli; dp-number of observed D. portschinskii; Dis. From River in Meters -the distance of a plot from the brook; slmax -average slope at a plot; cst -convergence index; Surface -surface type (see Fig. 2c).