Habitat Partitioning and Overlap by Large Lacertid Lizards in Southern Europe

South-western Europe has a rich diversity of lacertid lizards. In this study, we evaluated the occupancy patterns and niche segregation of five species of lacertids, focusing on large-bodied species (i.e., adults having >75 mm snout-vent length) that occur in south-western Europe (Italian to the Iberian Peninsula). We characterized the niches occupied by these species based on climate and vegetation cover properties. We expected some commonality among phylogenetically related species, but also patterns of habitat segregation mitigating competition between ecologically equivalent species. We used multivariate ordination and probabilistic methods to describe the occupancy patterns and evaluated niche evolution through phylogenetic analyses. Our results showed climate niche partitioning, but with a wide overlap in transitional zones, where segregation is maintained by species-specific responses to the vegetation cover. The analyses also showed that phylogenetically related species tend to share large parts of their habitat niches. The occurrence of independent evolutionary lineages contributed to the regional species richness favored by a long history of niche divergence.


Introduction
Climate is a powerful environmental factor driving the process of niche diversification in reptiles [1,2]. Tolerance to maximum temperatures in reptiles is evolutionarily constrained, possibly because of the importance of external heat sources in maintaining activity and for bodily water balance [3,4]. For these reasons, there is a significant association between the composition of reptile assemblages and thermal latitudinal gradients [5]. However, the thermoregulation efficiency of reptiles is not only mediated by the overall climate conditions, but also by the temperature conditions in microhabitats [6]. Vegetation cover and structure regulate the patterns of reptile occurrence, but at a finer scale than the climate [7,8].
South-western Europe encompasses a relatively rich reptile fauna, favored by its topographic heterogeneity, insularity and mild climate conditions [9]. In this region, several groups of phylogenetically related species display complex patterns of overlap structured by environmental gradients or by interspecific interactions [10][11][12]. In this study we evaluated the niche occupancy patterns of five species of large lacertid lizard (i.e., adults having >75 mm snout-vent length) that occur in south-western Europe and comprise a monophyletic group [13]. We focused on phylogenetically related species because stronger competitive interactions among them can be expected [14]. The target species in this study included one species having a broad circumboreal distribution, Lacerta agilis (Linnaeus 1758), two western Mediterranean endemics Lacerta bilineata (Daudin, 1802) and Timon lepidus (Daudin, 1802), and two Iberian endemics Lacerta schreiberi (Bedriaga, 1878) and Timon nevadensis (Buchholz, 1963). Given these substantial chorological differences, it was also expected that these species would diverge in their environmental associations. For example, T. nevadensis mainly occupies semi-arid steppe-like habitats, L. agilis appears to be confined to sub-alpine meadows, and the other species are possibly more habitat generalists [15,16].
The purpose of this study was to investigate the patterns of habitat occupancy of these five species of large lacertids and evaluate them from an evolutionary perspective. We hypothesized that the diversity of large lacertid lizards in southern Europe would be favored by niche partitioning. However, this partitioning will be phylogenetically constrained because related species tend to occupy similar or equivalent niches [17]. To test these hypotheses, we used multivariate ordination and novel probabilistic methods that enabled quantification of the niche overlap between species, and the niche breadth.

Study Region and Surveys
The study region encompassed most of south-western Europe including the Iberian Peninsula, southern France and the Italian Peninsula ( Figure 1). The region is dominated by Mediterranean climate types, ranging from subtropical warm desert to humid sub-Mediterranean and oceanic, and temperate to tundra-like types in the mountain ranges (Pyrenees, Alps) [18]. This environmental heterogeneity favors the concentration of high biotic diversity in this region, including species having xeric Mediterranean and mesotemperate affinities [19].
ity 2021, 13, x FOR PEER REVIEW 2 of 11 chorological differences, it was also expected that these species would diverge in their environmental associations. For example, T. nevadensis mainly occupies semi-arid steppe-like habitats, L. agilis appears to be confined to sub-alpine meadows, and the other species are possibly more habitat generalists [15,16]. The purpose of this study was to investigate the patterns of habitat occupancy of these five species of large lacertids and evaluate them from an evolutionary perspective. We hypothesized that the diversity of large lacertid lizards in southern Europe would be favored by niche partitioning. However, this partitioning will be phylogenetically constrained because related species tend to occupy similar or equivalent niches [17]. To test these hypotheses, we used multivariate ordination and novel probabilistic methods that enabled quantification of the niche overlap between species, and the niche breadth.

Study Region and Surveys
The study region encompassed most of south-western Europe including the Iberian Peninsula, southern France and the Italian Peninsula ( Figure 1). The region is dominated by Mediterranean climate types, ranging from subtropical warm desert to humid sub-Mediterranean and oceanic, and temperate to tundra-like types in the mountain ranges (Pyrenees, Alps) [18]. This environmental heterogeneity favors the concentration of high biotic diversity in this region, including species having xeric Mediterranean and meso-temperate affinities [19]. Species records were obtained opportunistically, based on random habitat surveys conducted in the region from spring to autumn, following the annual activity periods. The surveys were planned to capture the maximum heterogeneity of habitats within the distribution ranges of these species, but with no a priori selection of the most suitable habitats (i.e., a random sampling of available habitats). In total, 823 sites were surveyed throughout south-western Europe by 1-3 observers; each site was visited only once. The occurrence of each species was assessed based on visual surveys and rock flipping, because both techniques have been used to build inventories of diurnal lizards [20]. The surveys were conducted on sunny days between 10:00 a.m. and 5:00 p.m. local time. The visual surveys were complemented with rock flipping in cases where identification of the species could not be visually ascertained. In total, 188 records of five species of large lac- Species records were obtained opportunistically, based on random habitat surveys conducted in the region from spring to autumn, following the annual activity periods. The surveys were planned to capture the maximum heterogeneity of habitats within the distribution ranges of these species, but with no a priori selection of the most suitable habitats (i.e., a random sampling of available habitats). In total, 823 sites were surveyed throughout south-western Europe by 1-3 observers; each site was visited only once. The occurrence of each species was assessed based on visual surveys and rock flipping, because both techniques have been used to build inventories of diurnal lizards [20]. The surveys were conducted on sunny days between 10:00 a.m. and 5:00 p.m. local time. The visual surveys were complemented with rock flipping in cases where identification of the species could not be visually ascertained. In total, 188 records of five species of large lacertids were obtained ( Figure 1) and were distributed as follows: L. agilis (number of records = 29), L. bilineata (59), L. schreiberi (8), T. lepidus (80), and T. nevadensis (12).

Environmental Data
The niches for the species were described based on climate and vegetation. We used three variables to describe the climate: the maximum temperature of the warmest month, the minimum temperature of the coldest month, and the accumulated annual precipitation, provided by the WorldClim database [21]. These variables represent climatic properties that reportedly influence reptile ranges because they describe thermal extremes and available environmental moisture [22].
The influence of plant cover was assessed using a surrogate for vegetation primary productivity, the enhanced vegetation index (EVI; [23]); EVI values correlate positively with the density of trees [24]. The EVI data were obtained for the 2009-2019 period from the high-resolution (250 m pixel -1 ) MODIS Collection EVI composite images [25]. The MODIS data were first checked to remove atmospheric artefacts, and then used to generate a series of variables describing the seasonal variability among habitats [26] including: the mean value (EVImean), the coefficient of variation (EVIcv), and the range (EVIrange), considering the inter-annual (mean value for 10 years) and spatial variability (mean value for 50 points, generated randomly within a maximum radius of 5 km). This larger area assessed the effect of the environment around the core habitat where specimens were found and took into account that species occurrence is sustained by isolated suitable habitat patch, but also by the interconnection of habitat patches that support the entire population [27].

Data Analysis
Species associations with two niche dimensions (climate and habitat) were visualized using outlying mean index (OMI) analysis [28]. The OMI ordination describes the species responses by quantifying their ecological marginality (the distance between the species centroid and the mean environmental conditions [28]). An OMI value close to zero indicates a higher similarity between the species position and the background environmental conditions. The OMI analysis also provided an estimate of the niche breadth of the species (tolerance index), and the proportion of the environmental variance explained by the OMI axes (residual tolerance; [28]). These analyses were carried out using the software package ade-4 [29] for R [30].
The niche overlap was estimated between pairs of species for a probabilistic niche region [31]. This overlap index was generated after 10,000 Monte Carlo draws for a niche region (alpha = 0.95), using the predictor variables. This analysis evaluates the probability that an individual of species X is found in the habitat of species Y, and vice versa, and produces two index values for a single pair of species (i.e., X→Y and Y→X) [32]. This method has the advantage of being weakly sensitive to the sample size [31], which is useful when evaluating the ecological overlap between species having dissimilar distributions, as was the case in this study. These analyses were carried out using the nicheROVER software package [32] for R.
We also compared the habitat characteristics between the pairs of species. Before modelling these associations, we tested the predictor variables for spatial autocorrelation using Moran's I correlograms [33]. Moran's I values were statistically significant, varying from 0.16 (EVIcv) to 0.55 (maximum temperature). To remove spatial autocorrelation, we built Binomial Generalized Linear Auto-Covariate Models (BGLAMs; [34]) selecting the best candidate model using the Akaike information criterion corrected for small sample sizes (AICc; [35]). In general, the best candidate models show lower AICc values, and an AICc weight ≥0.1 [35]. These analyses were carried out using the software packages spdep [36] and AICcmodavg [37] in the R environment.

Molecular Phylogenetics and Biogeography
Evolutionary relationships among the species of Lacerta and Timon were assessed by building a phylogenetic tree generated using Bayesian analysis of mitochondrial cytochrome b and 12s genes, obtained from GenBank (Supplementary Materials; Appendix I). The sequences were assembled and aligned using Bioedit 7.09 [38]. Our dataset comprised 32 sequences of variable length combining the two genes, and represented nine Lacerta and six Timon species, and four outgroups. The analysis of DNA evolution was conducted using jModelTest [39] and showed that the GTR+I+G model was the best for both the 12s and cytochrome b genes. We used three points of calibration to establish the times of divergence among species. The emergence of the Canary Island of El Hierro representing the divergence of Gallotia caesaris caesaris and Gallotia caesaris gomerae (1.0 Mya [40]) was used selecting a normal prior distribution (sigma = 0.02). The split between Lacerta and Timon was dated based using as a prior a gamma distribution based on a minimum age of 17.5 Mya [41] with shape and scale set to 1.0. The same prior distribution was used to calibrate the divergence between L. viridis and L. bilineata 8.7 Mya [42]. Bayesian analyses were performed using BEAST v 2.6.3. [43] running two chains of 5 × 10 8 iterations, sampling every 10,000 iterations. Chains were checked for convergence and ESS using Tracer 1.5 [44] and were combined after a burn-in of 99%. To reconstruct the biogeographic history of large lacertid lizards in south-western Europe we used ancestral range estimation using Bio-GeoBEARS implemented in RASP 4.2 [45]. The regions included were south-western Asia (Caucasus, Anatolia, Iran and the Middle East), Europe (excluding the Iberian Peninsula), north-western Africa and the Iberian Peninsula. The models of vicariance, dispersal and extinction allowing all the combinations of ancestral areas (except the Iberian Peninsula plus Asia) were evaluated using the Akaike criterion [46].

Results
The first two axes of the OMI explained 96.98% of total inertia (axis 1: 70.47%, axis 2: 26.51%). The first axis described a gradient from higher to lower temperature and precipitation and differentiated those species that occur under humid-cold conditions from those that occur under hot-dry conditions (Figure 2). The second axis described a transition between habitats having different vegetation cover, typically distinguishing habitats having a relatively high EVI and a low seasonal coefficient of variation (CV) (e.g., forests) from those having a relatively low EVI and a high seasonal CV (i.e., grasslands/cultivated lands; Figure 2). The genus Timon was separated from the genus Lacerta mainly along the first axis, the former showing a positive association with dry-warm climates (Figure 2). The niche indices indicated that a large part of the variation in the occurrence of species was explained by the environmental variables, with residual tolerance values between 19.2% (L. agilis) to 66.9% (L. schreiberi) ( Table 1). In general, the species showed moderate to high distances from the environmental centroid, ranging from 20.2% (L. schreiberi) to 74.4% (L. agilis) (Table 1), which indicated that they occupied confined subspaces within the available environment ( Figure 2). The tolerance indices consistently showed moderate to low values, ranging from 22.6% (T. lepidus) to 6.4% (L. agilis) ( Table 1), indicating that these species differed in their niche sizes (Figure 2).   Figures 3 and 4 and Table 2). Among phylogenetically more distant species the patterns were complex, including high (e.g., L.  Table 2). Between some pairs of species, the overlap was highly asymmetric (e.g., L. agilis-L. bilineata, T. lepidus-L. schreiberi, and T. lepidus-T. nevadensis) indicating that the niche of species Y (smaller niche) was partially nested within that of species X (larger niche) (Figure 2).  The niche overlap indices showed high values (>60) between pairs of sister species; for example, T. nevadensis → T. lepidus: 83.03 (Figures 3 and 4 and Table 2). Among phylogenetically more distant species the patterns were complex, including high (e.g., L. schreiberi → L. bilineata: 61.59; L. bilineata → T. lepidus: 86.03), moderate (L. agilis → L. bilineata: 53.54), and low (L. bilineata → T. nevadensis: 4.28) levels of overlap, and in some cases no overlap (L. agilis → T. nevadensis: 0.0) (Figure 3 and 4 and Table 2). Between some pairs of species, the overlap was highly asymmetric (e.g., L. agilis-L. bilineata, T. lepidus-L. schreiberi, and T. lepidus-T. nevadensis) indicating that the niche of species Y (smaller niche) was partially nested within that of species X (larger niche) (Figure 2). ity 2021, 13, x FOR PEER REVIEW 6 of 11       The BGLAMs showed than the niche separation between L. agilis-L. bilineata was mainly related to the by maximum-minimum temperatures (R 2 = 0.257) ( Table 3). The separation between L. agilis-T. lepidus was also attributed to temperature, and to a lesser degree to plant cover (R 2 = 0.469). In contrast, the separation between L. bilineata-T. lepidus was mainly related to plant cover and to a lesser extent to the maximum temperature and annual precipitation (R 2 = 0.181). Comparisons that included L. schreiberi and T. nevadensis produced statistically poorly supported models and are not shown. Table 3. Binomial Generalized Linear Auto-Covariate Models (BGLAMs) evaluating the environmental separation between pairs of species which geographically contact. Lacerta schreiberi and T. nevadensis were not included in the analyses. AIC, Akaike information criterion; AICWt, AIC weights. LAG = L. agilis; LBI = L. bilineata; TLE = T. lepidus. T, temperature; Prec, precipitations; EVI, Enhanced Vegetation Index; m, mean; cv, coefficient of variation. The Bayesian phylogenetic analysis strongly supported the monophyly of the genera Timon and Lacerta, and most of the relationships across species within these genera ( Figure 4). The most supported biogeographic model was DIVALIKE +J (AICc = 79.3, AICc weight = 0.630), indicating a founder event (j = 0.047) and low rates of dispersal and extinction (both <0.0001). Based on ancestral range reconstruction, large European lacertids arose in western Asia and independently colonized south-western Europe on four occasions. The oldest colonization event occurred approximately 8.9-15.5 Mya, during the invasion of the western Iberian Peninsula from Europe by the ancestor of L. schreiberi. The ancestor of Timon emigrated from western Asia to north-western Africa, and subsequently invaded the Iberian Peninsula approximately 6.5-13.0 Mya. The split of the ancestral Iberian species into T. lepidus and T. tangitanus probably occurred 4.6-11.5 Mya. The split of L. bilineata from the shared ancestor with L. viridis occurred approximately 8.7-9.4 Mya, possibly after their isolation in the Balkan and Italian peninsulas. Lacerta bilineata recently colonized the Iberian Peninsula from the Italian Peninsula (300,000-550,000 years ago) and the widespread palaearctic species L. agilis, colonized to the eastern Pyrenees from central Europe approximately 0.4-1.8 Mya.

Discussion
The OMI analysis showed that in south-western Europe the large lacertid lizards differed in niche marginality and tolerance. In general, environmental variables explained a major part of the variability in the occurrence of these species (residual tolerance 19.2% to 53.3%); an exception was L. schreiberi (residual tolerance 66.9%), possibly as a result of the low number of records. The ecological diversification of these species has possibly been driven by physiological tolerance, particularly variability in preferred temperature and the level of environmental moisture [47].
The analysis also indicated that T. lepidus and L. bilineata had the highest ecological tolerances, potentially enabling them to occupy parts of the niches of the other species. In contrast, T. nevadensis and L. agilis, occupied relatively narrow niches, at opposite extremes of the environmental range. Ecological partitioning was evident, but incomplete, between L. agilis and L. bilineata (maximum overlap probability 53.54%) and between L. agilis-T. lepidus (maximum overlap probability 14.76%). These species pairs consistently showed either geographical overlap (L. agilis-L. bilineata 170,990 km 2 ) or were almost completely parapatric (L. agilis-T. lepidus).
Regression models showed that the segregation between L. agilis and these species of lacertids (L. bilineata and T. lepidus) was mainly influenced by the variables describing temperature conditions, with the occurrence of L. agilis being negatively associated with temperature. This is consistent with the relict status of this species in the region, where it appears mainly isolated to mountain ranges, at altitudes above 1300 m [48]. The BGLAMs revealed that vegetation cover was unrelated to the niche partitioning between L. agilis and L. bilineata, possibly because of the use of open habitats in harsh subalpine environments by the later species [15,49].
There was wide but asymmetric overlap in the niches of T. lepidus and L. bilineata, and the niche of the latter was partially nested within that of T. lepidus. Both species also coexist over a wide geographic range (110,203 km 2 ), but T. lepidus is more widespread in the Iberian Peninsula. The later arrival of L. bilineata to the Iberian Peninsula may have been a disadvantage for this species, with it being excluded from potentially suitable habitats by the other lacertid species that evolved in this region (T. lepidus, L. schreiberi). Our results indicated that in the zone of coexistence between T. lepidus and L. bilineata, the species were largely segregated according to the vegetation cover. In general, forest habitats constitute unfavorable habitat for temperate lacertid lizards, which regulate their body temperature by sun basking [50]. The canopy of temperate forests greatly reduces light transmittance to the lower understory layers [51]. For this reason, closed forests are usually avoided by lacertid lizards, although they can colonize discontinuities in these habitats including path edges, interspersed meadows, forest margins and rocky outcrops [52][53][54]. Lacerta bilineata exploits microhabitats with a very dense vegetation cover, because this lizard thermoregulates efficiently using bushes and tree logs as basking platforms [47,55].
OMI and overlap analyses showed similarities in the patterns of habitat occupancy, between the species pairs T. lepidus-T. nevadensis and L. schreiberi-L. bilineata, and T. nevadensis niche was almost completely nested within that of T. lepidus. Timon nevadensis typically occupies relatively sparsely vegetated semiarid habitats, but T. lepidus is also able to exploit open, de-vegetated habitats in regions where T. nevadensis does not occur (e.g., in the shrub-steppes of the Ebro valley; [56]). These findings suggest that the range limits of both species may be sustained by interspecific interactions in the contact zones rather than by ecotonal transitions [57]. Our results showed a substantial habitat overlap between L. schreiberi-L. bilineata, which may trigger competitive interactions where both species contact geographically [58,59]. The results did not indicate that the species of Iberian origin (L. schreiberi, T. lepidus, and T. nevadensis) have greater habitat overlap than those of recent arrivals (e.g., L. bilineata and L. agilis). It is possible that during the prolonged isolation in the Iberian Peninsula these species diverged in their environmental niches, to reduce interspecific competition. Together these results indicate that the interaction of several mechanisms (i.e., interspecific competition, species evolutionary history, ecophysiological tolerance) determine the occurrence of large lizard species in the region, in a similar way to that observed for other lizard assemblages [60,61]. However, one limitation of our study was that it was not able to discern the relative importance of autoecological and synecological aspects of the distribution patterns.

Conclusions
The Iberian Peninsula has the richest lacertid fauna in south-western Europe, supporting five species of large lacertids. Our results revealed that this species richness is favored by climate niche partitioning, but with transitional areas of overlap where the segregation is maintained by species-specific responses to the level of vegetation cover. Interspecific competition may also play a key role in the patterns of occurrence, with the species that have arrived more recently on the Iberian Peninsula having been excluded from potentially favorable habitats because of prior habitat occupancy by species that evolved in this region. The occurrence of several independent evolutionary lineages has partly contributed to the species richness, which has been favored by a long history of ecological divergence between subclades having distinct geographic origins. The analyses applied in this study have the advantage of being weakly sensitive to sample size and are robust to spatially aggregated records, so they may be useful in disentangling the patterns of niche partitioning in assemblages including ecologically (or phylogenetically) related species structured by the effect of various interacting factors (i.e., competition, evolutionary history, environmental tolerance), which is characteristic of biotic communities in the Mediterranean region [62,63].