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

Natia Barateli; David Tarkhnishvili; Giorgi Iankoshvili; Luka Kokiashvili; Nikoloz Dvali; Zurab Janiashvili

doi:10.3897/herpetozoa.34.e63072

Abstract

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.

Key Words

competition, ecological niche, parthenogenesis, rock lizards, territoriality

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 inter- and 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; Pfennig and Pfennig 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 (2011, 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 microgeographic 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.

Materials and methods

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 long-term 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.

Figure 1.

Study area. a: satellite image downloaded from Google Earth Pro (2020); b: Orthophoto image prepared with the drone (see the methods section). The arrow points in the northern direction. The brook borders the rocky habitat from the north.

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 PhotoScan 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.

Figure 2.

a: the numbers of studied plots; b: DEM image (0.003 m/ pixel); c: surface types (see the text for the explanations).

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 QGIS Development Team). 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² 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 levels 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.

Estimation of niche overlap

MacArthur and Levins (1967) suggested several baselines for the analysis of competition and niche displacement. They offered to use an asymmetric coefficient of niche overlap,

$O_{21} = \sum_{i = 1}^{n} p_{2 i} p_{11} / \sum_{i = 1}^{n} {(p_{2 i})}^{2}$ .

In this equation, O₂₁ 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.

Results

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.

Figure 3.

Spatial distribution of D. dahli and D. portschinskii in the studied rocky outcrop (aggregated data on the lizard sightings throughout the entire study period).

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.

Figure 4.

Dependence of the average and 95% confidence intervals for the number of recorded individuals (NL) of D. dahli (solid lines) and D. portschinskii (dashed lines) at a plot during a given session on: a: distance from the brook (m), b: year, c: air temperature (°C), d: days passed after 1^st April, e: substrate type and f: average slope index on the 3 cm pixel scale (averaged for a plot).

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

Table 1.

Download as

CSV

XLSX

The outcome of the GLMM analysis*.

Source of variance	Type	F (DD)	F (DP)	df1	df2	Sig (DD)	Sig (DP)
Corrected model		22.34	2.18	20	1.17	<0.001	0.002
Distance from the brook	Fixed in space	9.75	8.79	1	1.17	0.002	0.003
Surface type	Fixed in space	3.66	1.59	11	1.17	<0.001	0.097
Convergence Index	Fixed in space	0.87	0.21	1	1.17	0.352	0.645
Average Slope	Fixed in space	6.03	0.39	1	1.17	0.014	0.535
Year	Fixed in time	30.90	0.67	2	1.17	<0.001	0.513
Date	Fixed in time	28.99	1.32	1	1.17	<0.001	0.250
Air Temperature	Fixed in time	10.17	0.03	1	1.17	<0.001	0.864
Surface Humidity	Varying	16.00	1.50	1	1.17	<0.001	0.221
Surface Temperature	Varying	2.89	0.01	1	1.17	0.089	0.960
	Coefficients (for significant predictors only)
	b (DD)	b (DP)		Sig(DD)		Sig(DP)
Surface = forest	-2.764	N/A		0.004		N/A
Surface = grassy rock	-1.269	N/A		0.050		N/A
Surface = rock	0.765	0.070		0.012		0.016
Surface = Gravel	N/A	0.126		N/A		0.002
Distance from the brook	-0.006	0.001		0.002		0.003
Date (days from 1 April)	-0.013	N/A		<0.001		N/A
Surface Humidity	-0.025	N/A		<0.001		N/A
Air Temperature	-0.089	N/A		0.001		N/A
Average Slope	0.022	N/A		0.014		N/A

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. 2010; Petrosyan 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, 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.

Acknowledgements

The research was supported by Shota Rustaveli National Science Foundation (Award# 217478).The authors appreciate Akaki Nadaraia for assisting in drone operation in field conditions and Eduard Galoyan for useful suggestions during the field work. Günter Gollmann and two anonymous reviewers provided useful comments on the first draft of the paper.

References

Abrams PA (1986) Character displacement and niche shift analyzed using consumer-resource models of competition. Theoretical Population Biology 29: 107–160. https://doi.org/10.1016/0040-5809(86)90007-9

Angilletta MJ, Zelic MH, Adrian GJ, Hurliman AM, Smith CD (2013) Heat tolerance during embryonic development has not diverged among populations of a widespread species (Sceloporus undulatus). Conservation Physiology 1: cot018. https://doi.org/10.1093/conphys/cot018

Bakradze MA (1977) Sympatric populations of parthenogenetic and bisexual species of rock lizards of genus Lacerta in Georgia. PhD Dissertation. Iv. Javakhishvili Tbilisi State University. [In Russian]

Begon M, Harper JL, Townsend CR (1996) Ecology. Individuals, populations and communities. NY: Blackwell Scientific Publications.

Bell G (2008) Selection: the mechanism of evolution. Oxford University Press, Oxford.

Brown JL (1969) Territorial behavior and population regulation in birds. Wilson Bulletin 81: 293–329.

Carvalho JC, Cardoso P (2020) Decomposing the causes for niche differentiation between species using hypervolumes. Frontiers in Ecology and Evolution 8: e243. https://doi.org/10.3389/fevo.2020.00243

Case TJ (1990) Invasion resistance arises in strongly interacting species-rich model competition communities. Proceedings of the National Academy of Sciences of the USA 87: 9610–9614. https://doi.org/10.1073/pnas.87.24.9610

Case TJ, Bender EA (1981) Testing for higher order interactions. American Naturalist 118: 920–929. https://doi.org/10.1086/283885

Colwell RK, Futuyma DJ (1971) On the measurement of niche breadth and overlap. Ecology 52: 567–576. https://doi.org/10.2307/1934144

Danielyan F, Arakelyan M, Stepanyan I (2008) Hybrids of Darevskia valentini, D. armeniaca and D. unisexualis from a sympatric population in Armenia. Amphibia-Reptilia 29: 487–504. https://doi.org/10.1163/156853808786230424

Darevsky IS, Kulikova VN (1962) Taxonomic characters and certain peculiarities of the oogenesis of hybrids between bisexual and parthenogenetic forms of Lacerta saxicola Eversmann. Cytologia 5: 160–170.

Darevsky IS (1967) Skal’nye yashcheritsy Kavkaza. Leningrad: Nauka (Translated as: Rock lizards of the Caucasus. New Delhi: Indian National Scientific Documentation Centre).

Dobzhansky T (1937) Genetic nature of species differences. American Naturalist 71: 404–420. https://doi.org/10.1086/280726

Endler JA (1986) . Natural selection in the wild. Princeton University Press, Princeton.

Freitas SN, Harris DJ, Sillero N, Arakelyan M, Butlin RK, Carretero MA (2019) The role of hybridisation in the origin and evolutionary persistence of vertebrate parthenogens: a case study of Darevskia lizards. Heredity 123: 795–808. https://doi.org/10.1038/s41437-019-0256-5

Fretwell SD (1969) On territorial behavior and other factors influencing habitat distribution in birds. Acta Biotheoretica 109: 45–52. https://doi.org/10.1007/BF01601955

Frishkoff LO, Hadly EA, Daily GC (2015) Thermal niche predicts tolerance to habitat conversion in tropical amphibians and reptiles. Global Change Biology 21: 3901–3916. https://doi.org/10.1111/gcb.13016

Gabelaia M, Tarkhnishvili D, Murtskhvaladze M (2015) Phylogeography and morphological variation in a narrowly distributed Caucasian rock lizard, Darevskia mixta. Amphibia-Reptilia 36: 45–54. https://doi.org/10.1163/15685381-00002975

Galoyan EA (2011) The role of social interactions in the density regulation in populations of parthenogenetic and bisexual species of rock lizards. PhD Dissertation. M. Lomonosov Moscow State University. [In Russian]

Galoyan E (2013) Joint space use in a parthenogenetic Armenian rock lizard (Darevskia armeniaca) suggests weak competition among monoclonal females. Journal of Herpetology 47: 97–104. https://doi.org/10.1670/11-242

Galoyan E, Bolshakova A, Abrahamyan M, Petrosyan R, Komarova V, Viсtor S, Arakelyan M (2019) Natural history of Valentin’s rock lizard (Darevskia valentini) in Armenia. Zoological Research 40: 277–292. https://doi.org/10.24272/j.issn.2095-8137.2019.036

Gauze GF (1934) The mathematical theory of the struggle for existence and its application to populations of yeast cells. Bulleten Moskovskogo Obshchestva Ispytatelei Prirody, Otdel Biologicheskii 43: 69–87.

Glesener RR, Tilman D (1978) Sexuality and the components of environmental uncertainty: clues from geographic parthenogenesis in terrestrial animals. American Naturalist 112: 659–673. https://doi.org/10.1086/283308

Google earth Pro (2020) Google earth V 7.3.3.7786. Kojori, Georgia. 41°38'56.65"N, 44°41'00.13"E, Eye alt 753 feet. http://www.earth.google.com [September 25, 2020]

Grant PR (1972) Convergent and divergent character displacement. Biological Journal of the Linnean Society 4: 39–68. https://doi.org/10.1111/j.1095-8312.1972.tb00690.x

Haag CR, Ebert D (2004) A new hypothesis to explain geographic parthenogenesis. Annali Zoologii Fennici 42: 539–544.

Holt RD (1987) On the relation between niche overlap and competition: the effect of incommensurable niche dimensions. Oikos 1987: 110–114. https://doi.org/10.2307/3565696

Huey RB, Stevenson RD (1979) Integrating thermal physiology and ecology of ectotherms: a discussion of approaches. American Zoologist 19: 357–366. https://doi.org/10.1093/icb/19.1.357

Hurlbert SH (1978) The measurement of niche overlap and some relatives. Ecology 59: 67–77. https://doi.org/10.2307/1936632

Hutchinson GE (1957) Concluding remarks. Cold Spring Harbor Symposium in Quantitative Biology 22: 415–27. https://doi.org/10.1101/SQB.1957.022.01.039

Kassahn KS, Crozier RH, Pörtner HO, Caley MJ (2009) Animal performance and stress: responses and tolerance limits at different levels of biological organisation. Biological Reviews 84: 277–292. https://doi.org/10.1111/j.1469-185X.2008.00073.x

Kearney M, Wahl R, Autumn K (2005) Increased capacity for sustained locomotion at low temperature in parthenogenetic geckos of hybrid origin. Physiological and Biochemical Zoology 78: 316–324. https://doi.org/10.1086/430033

Kearney M, Fujita MK, Ridenour J (2009) Lost sex in the reptiles: constraints and correlations. In: Schön I, Martens K, Dijk P (Eds) Lost sex. Dordrecht: Springer; 2009, 447–474. https://doi.org/10.1007/978-90-481-2770-2_21

Kondrashov AS (1982) Selection against harmful mutations in large sexual and asexual populations. Genetics Research 40: 325–332. https://doi.org/10.1017/S0016672300019194

Kuzmin SL, Tarkhnishvili DN (1992) Food competition in Anura after metamorphosis. Bulleten Moskovskogo Obshchestva Ispytatelei Prirody, Otdel Biologicheskii 97: 53–63.

MacArthur R, Levins R (1967) The limiting similarity, convergence, and divergence of coexisting species. American Naturalist 101: 377–385. https://doi.org/10.1086/282505

MacArthur R (1970) Species packing and competitive equilibrium for many species. Theoretical Population Biology 1: 1–11. https://doi.org/10.1016/0040-5809(70)90039-0

Maynard-Smith J (1998) Evolutionary Genetics, 2^nd edn. Oxford University Press, Oxford.

Murphy RW, Fu J, Macculloch RD, Darevsky IS, Kupriyanova LA (2000) A fine line between sex and unisexuality: the phylogenetic constraints on parthenogenesis in lacertid lizards. Zoological Journal of the Linnean Society 130: 527–549. https://doi.org/10.1111/j.1096-3642.2000.tb02200.x

Partridge L (1978) Habitat selection. In: Krebs JR, Davies NB (Eds) Behavioral ecology: and evolutionary approach. Blackwell, Oxford, 351–376.

Petren K, Bolger DT, Case TJ (1993) Mechanisms in the competitive success of an invading sexual gecko over an asexual native. Science 259: 354–358. https://doi.org/10.1126/science.259.5093.354

Petrosyan V, Osipov F, Bobrov V, Dergunova N, Omelchenko A, Varshavskiy A, Danielyan F, Arakelyan M (2020) Species distribution models and niche partitioning among unisexual Darevskia dahli and its parental bisexual (D. portschinskii, D. mixta) Rock Lizards in the Caucasus. Mathematics 8: e1329. https://doi.org/10.3390/math8081329

Pfennig DW, Pfennig KS (2010) Character displacement and the origins of diversity. American Naturalist 176: S26–S44. https://doi.org/10.1086/657056

Pianka ER (1976) Natural selection of optimal reproductive tactics. American Zoologist 16: 775–784. https://doi.org/10.1093/icb/16.4.775

Pianka ER (2011) Evolutionary ecology. 6^th edn. Pearson.

Pulliam HR (1988) Sources, sinks, and population regulation. American Naturalist 132: 652–661. https://doi.org/10.1086/284880

Şahin MK, Candan K, Caynak EY, Kumlutaş Y, Ilgaz Ç (2020) Ecological niche divergence contributes species differentiation in worm lizards (Blanus sp.) (Squamata: Amphisbaenia: Blanidae) in Mediterranean part of Anatolian peninsula and the Levantine region. Biologia 2020: 1–8. https://doi.org/10.2478/s11756-020-00548-1

Schoener TW (1974) Competition and the form of habitat shift. Theoretical Population Biology 6: 265–307. https://doi.org/10.1016/0040-5809(74)90013-6

Schoener TW (1983) Field experiments on interspecific competition. American Naturalist 122: 240–285. https://doi.org/10.1086/284133

Spangenberg V, Arakelyan M, Galoyan E, Matveevsky S, Petrosyan R, Bogdanov Y, Danielyan F, Kolomiets O (2017) Reticulate evolution of the rock lizards: meiotic chromosome dynamics and spermatogenesis in diploid and triploid males of the genus Darevskia. Genes 8: e149. https://doi.org/10.3390/genes8060149

Stamps JA (1991) The effect of conspecifics on habitat selection in territorial species. Behavioral and Ecological Sociobiology 28: 29–36. https://doi.org/10.1007/BF00172136

Tarkhnishvili D, Gavashelishvili A, Avaliani A, Murtskhvaladze M, Mumladze L (2010) Unisexual rock lizard might be outcompeting its bisexual progenitors in the Caucasus. Biological Journal of the Linnean Society 101: 447–460. https://doi.org/10.1111/j.1095-8312.2010.01498.x

Tarkhnishvili D (2012) Evolutionary history, habitats, diversification, and speciation in Caucasian rock lizards. In: Jenkins R (Ed. ) Advances in Zoology Research 2: 79–120.

Tarkhnishvili D, Murtskhvaladze M, Anderson CL (2017) Coincidence of genotypes at two loci in two parthenogenetic rock lizards: how backcrosses might trigger adaptive speciation. Biological Journal of the Linnean Society 121: 365–378. https://doi.org/10.1093/biolinnean/blw046

Tarkhnishvili D, Yanchukov A, Şahin MK, Gabelaia M, Murtskhvaladze M, Candan K, Galoyan E, Arakelyan M, Iankoshvili G, Kumultas Y, Ilgaz C, Matur F, Colak F, Erdolu M, Kurdadze S, Barateli N, Anderson CL (2020) Genotypic similarities among the parthenogenetic Darevskia rock lizards with different hybrid origins. BMC Evolutionary Biology 20: 1–25. https://doi.org/10.1186/s12862-020-01690-9

Trofimov AG (1981) Population spatial structure of some parthenogenetic and gonochoristic rock lizards from the Caucasus. Ph.D. diss., Taras Shevchenko Kyiv University, USSR. [In Russian]

Tsellarius AY, Tsellarius EY (2001) Alteration of spacing pattern of population of Lacerta saxicola in broad-leaved forests of Navagir mountain ridge. Zoologicheskii Zhurnal 80: 1–8. [in Russian]

Tsellarius AY, Tsellarius EY (2003) An access to the females as a resource of male’s territory in Lacerta saxicola. In Herpetologia Petropolitana. Proceedings of the 12^th Ordinary General Meeting of the Societa Europaea Herpetologica, 222–225.

Tsellarius AY, Tsellarius EY, Galoyan EA (2016) Social relationships between males and females in the rock lizard (Darevskia brauneri, Lacertidae) 1. Friendly monogyny in males and polyandry in females. Biology Bulletin 43: 1077–1086. https://doi.org/10.1134/S106235901609017X

Vrijenhoek RC, Parker ED (2009) Geographical parthenogenesis: general purpose genotypes and frozen niche variation. In: Schön et al. (Eds) Lost Sex. Springer, Dordrecht, 99–131. https://doi.org/10.1007/978-90-481-2770-2_6

White MJD (1970) Heterozygosity and genetic polymorphism in parthenogenetic animals. In Essays in evolution and genetics in honor of Theodosius Dobzhansky. Springer, Boston, MA, 237–262. https://doi.org/10.1007/978-1-4615-9585-4_8

Wright JW, Vitt LJ (1994) Biology of whiptail lizards (genus Cnemidophorus). Systematic Biology 43: 148–151. https://doi.org/10.1093/sysbio/43.1.148

Supplementary materials

Supplementary material 1

Table S1

Natia Barateli, David Tarkhnishvili, Giorgi Iankoshvili, Giorgi Iankoshvili, Nikoloz Dvali, Zurab Janiashvili

Data type: Microsoft Excel table, full data of the records and physical conditions

Explanation note: Full record for the three years of the observation. Pl – the number of a plot (see Fig. 2a). dd – the number of observed active of D. dahli, dp – the number of observed D. portschinskii, Hum – surface humidity index, Temp – surface temperature, Air.Temp – air temperature, Con-stdev – convergence index, Slope_max – average slope, Time of Day: M – early afternoon, N – before sunset, Surface – surface type as shown in Fig. 2c, N – northing, E – easting of the center of a plot.

This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.

Download file (279.04 kb)

Supplementary material 2

Figure S1

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

This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.

Download file (375.23 kb)