The organization of biotic communities is an ongoing area of research in ecology. These studies have shown that the diversity of communities is determined by the environment (e.g., by primary productivity or habitat complexity; Koh, Lee & Lin, 2006) and by the interactions established among co-occurring species (Petchey & Gaston, 2002). Phylogenetically related species share functional traits and tend to use similar resources (e.g., basking spots and prey), thereby becoming mutually exclusive (Stuart et al., 2014). Thus, biotic communities frequently exhibit a phylogenetic structure (Cavender-Bares, Keen & Miles, 2006).
Lizards comprise one of the most diverse reptile groups in warm temperate regions where they form rich communities in terms of species (Arnold, 1987). In these communities, congeneric species tend to be segregated, thereby reducing negative interactions via the differential use of microhabitats (Capula, Luiselli & Rugiero, 1993; Vanhooydonck, Van Damme & Aerts, 2000). Mediterranean skinks are a group of lizards characterized by semi-fossorial habits (Caputo, Lanza & Palmieri, 1995), with their communities being structured according to species morphology and the type of substrate (Attum, Eason & Cobbs, 2007).
In this study, we investigated habitat partitioning in two species of skinks that are endemic to south-western Europe i.e., Bedriaga’s Skink Chalcides bedriagai (Boscá, 1880) and the Western Three-toed Skink Chalcides striatus (Cuvier, 1829). These species occur throughout most of the Iberian Peninsula and the Mediterranean regions of southern France (Salvador, 1998; Cheylan et al., 2012). The ranges of both species overlap in a large part of this region (Pollo, 2004a; Pollo, 2004b) and they can occur in syntopy (Salvador, 1998; Cabana, 2010), thereby suggesting that they use the same habitats, at least occasionally. Chalcides striatus has a slender body shape and small limbs with a reduced number of digits, which facilitate its movement through dense layers of grass (Caputo, Lanza & Palmieri, 1995). Chalcides striatus is frequently associated with wet grasslands and it has been suggested that is a habitat specialist (Pollo, 2004a; Cheylan et al., 2012). By contrast, C. bedriagai exhibits generalist morphology, with robust limbs and pentadactyl feet, and occurs in a broad range of habitats (Malkmus, 2004).
We characterized the niches occupied by these skinks based on macroecological (climate, topography, and substrate) and habitat/microhabitat (plant associations, vegetation structure and shelter types) variables. Some level of niche partitioning could be expected (hypothesis i) in a similar manner to other congeneric lizards that overlap geographically and share similar resources (Arnold, 1987). However, the fact that both species occur in syntopy suggests that there may also be habitat overlap, which could be restricted by ecotonal habitats (Leache & Cole, 2007) or a lack of spatial structure (hypothesis ii).
Materials and Methods
Study region and surveys
The study region was the Iberian Peninsula and south-east France, which encompass a large part of the distribution of both species (Fig. 1). Most of this region has a Mediterranean climate (types Csa and Csb; Köppen classification). However, humid/oceanic climates (types Cfa and Cfb) occur in the extreme north, semi-arid climates (types Bsk/Bsh) in the central and south-eastern parts, and desert climates (types Bwk/Bwh) on the south-eastern coast (AEMET, 2011).
We considered all possible habitat types (ranging from alpine grasslands to coastal dunes and peri-urban environments) throughout the known distributions of both species (Fig. 1). The occurrence of species was assessed based on visual surveys and rock flipping because both techniques have been used in previous studies of the presence/absence of diurnal reptiles (Losos, Marks & Schoener, 1993; McDiarmid et al., 2012). All of the sites were visited at least two times by two surveyors during the period of maximum activity for these species (early spring to early summer; Salvador, 1998). Field work was conducted over a period of 12 years as part of a study of the region’s herpetofauna (Escoriza et al., 2016). Permits for field work were provided by the Departamento de Conservação da Natureza e Florestas do Algarve and the Departament de Medi Ambient de Catalunya (ref. SF/574).
At the macroscale level, we used several climate and topography variables that influence the occurrence of reptiles in temperate-subtropical regions (Anderson, 1999). We also included the substrate because it forms part of the habitat of lizards (Zaady & Bouskila, 2002; Escoriza, in press). The climate was characterized based on an aridity index (mean annual precipitation/mean annual potential evapotranspiration; Trabucco & Zomer, 2009) and the mean annual temperature (Hijmans et al., 2004). The aridity index ranges between >0.65 (humid), 0.50–0.65 (dry sub-humid), 0.20–0.50 (semi-arid), 0.03–0.20 (arid) and <0.03 (hyper arid). The climate data were obtained at a resolution of 1,000 m pixel−1 (Hijmans et al., 2004). The topography was characterized based on elevation, which was determined in situ using a global positioning system (Garmin Etrex 10; Garmin Ltd., Olathe, KS, USA). We examined the substrate texture at the first standard depth (0–5 cm) based on the sand (grain size 50–2,000 µm) and clay (grain size <2 µm) content, and the soil depth to bedrock. The soil data were obtained at a resolution of 250 m pixel−1 (Hengl et al., 2017). The data of these variables were extracted by the package QGIS vs 2.18 (QGIS Development Team, 2017).
We also characterized the habitats and microhabitats of the species, where we measured 37 parameters to describe the composition and structure of the habitats at a fine spatial resolution (50–0.12 m; Cerqueira & Freitas, 1999; Guénnette & Villard, 2005). Data regarding plant density and types were collected by sampling along two transects, where each measured 50 m in length and they radiated out from the center (i.e., where we found a specimen of the target species) in two opposite directions (north and south). Spermatophytes were classified according to a standard botanic classification based on growth type, stem lignification, and life cycle (Flora Iberica, 2017). The growth types were trees (species with a maximum height >5 m), bushes (1–5 m), sub-shrubs (<1 m), and lianas (climbing plants). Stem lignifications separated woody plants (trees, bushes, sub-shrubs, and some lianas), suffrutices (similar to sub-shrubs, but only lignified at the base), and grasses (plants without woody parts). The life cycle types were perennial, deciduous, and semi-deciduous for woody plants, and annual, biennial, semi-perennial, and perennial grasses (Flora Iberica, 2017). Monilophytes (ferns and equisetums) were classified as semi-perennial or perennial types (Flora Iberica, 2017).
We also assessed some characteristics of the microhabitat (at 5–0.12 m), where these variables described the density of vegetation cover and the number and type of potential shelters (Cerqueira & Freitas, 1999). We measured the distance to a forest edge (in meters up to a maximum of 1,000 m), rock surface exposure (%), soil surface exposure (%), number of stones, number of stumps, number of trees with a diameter at breast height (DBH) greater than 0.75 m, canopy cover (%, measured at 1.5 m using a spherical crown densitometer; Forestry Suppliers, Inc., Jackson, MS, USA), woody plant and grass species and stem density, and litter depth (cm).
We first visualized the multidimensional niches of species using principal component analysis (PCA). In subsequent analyses, we tested the variables separately at three spatial levels, i.e., at the macroscale (1,000–250 m), mesoscale (50 m), and microscale (5–0.12 m) levels. This hierarchical approach allowed that the number of cases per variable was always greater than two, which is the minimum required for an adequate estimation of the coefficients in linear regression analyses (Austin & Steyerberg, 2015). The associations between the predictor variables and the presence of C. bedriagai and C. striatus were modeled using distance-based linear models (DistLM), where we transformed the dependent variable into a matrix of Sørensen distances (Clarke & Gorley, 2006). This analysis was used to generate the best subset of predictor variables based on a stepwise selection procedure and Akaike’s information criterion corrected for finite sample sizes (AICc; Burnham & Anderson, 2002). To assess whether the predictors had positive or negative effects on the dependent variable, we generated XY scatter plots using trend lines (Clarke & Gorley, 2006). These analyses were performed using PRIMER-E (PRIMER-E Ltd., Plymouth, UK).
We also estimated the site-level suitability. We expected that if the species occupied a well-defined ecological space (i.e., parapatric niches), then the predictions would have low cross-classification errors (Evans & Cushman, 2009). Niche suitability was estimated using random forests, which is a classification method that was designed to avoid overfitting (Broennimann et al., 2012). This property was important for our study because of the high number of variables included in the models. We set the model parameters for classification, where the variables were sampled randomly at each split as the square root of the total number of variables included in the model and we set the number of training trees to 10,000 (Liaw & Wiener, 2002). We assessed the performance of the model using the area under the receiver operating characteristic curve (AUC). The AUC values ranged from 0.5 for models with predictive ability similar to chance to 1.0 for models with perfect predictive ability (Araújo et al., 2005). These analyses were carried out with the package ‘randomForest’ (Liaw & Wiener, 2002) and ‘ROCR’ (Sing et al., 2015) in R (R Development Core Team, 2015).
The target species were detected at 49 sites in Spain, Portugal, and southern France (Fig. 1). At three sites, both species were detected syntopically (Fig. 1). Chalcides striatus was detected using visual surveys and rock flipping whereas C. bedriagai was detected mostly by rock flipping. The descriptive statistics for the environmental variables are shown in Table 1. These data showed that both species overlap in a broad range of environmental parameters, occupying sparsely wooded habitats with a dense cover of grasses and shrubs (Table 1). Figure 2 shows some examples of these habitats.
|Scale (m)||C. bedriagai||C. striatus|
|Elevation (m asl)||616 (1–1,601)||476 (1–1,491)|
|Mean annual temperature (° C)||1,000||14.1 (7.8–18.2)||13.7 (8.6–17.1)|
|Aridity index||1,000||0.55 (0.24–1.34)||0.74 (0.33–1.63)|
|Sand texture%||250||43.0 (33.0–57.0)||41.7 (31.0–56.0)|
|Clay texture%||250||24.6 (16.0–32.0)||23.1 (17.0–28.0)|
|Soil depth (cm)||250||1,297 (431–2,053)||1,290 (337–2,112)|
|Needle-leaved tree%||50||0.7 (0.0–4.6)||0.1 (0.0–1.2)|
|Deciduous broad-leaved tree%||50||0.9 (0.0–15.1)||1.6 (0.0–6.4)|
|Perennial broad-leaved tree%||50||0.9 (0.0–6.6)||1.7 (0.0–5.3)|
|Needle-leaved bush%||50||0.5 (0.0–6.0)||0.3 (0.0–2.8)|
|Deciduous broad-leaved bush%||50||0.6 (0.0–6.0)||2.6 (0.0–13.9)|
|Perennial broad-leaved bush%||50||14.6 (0.0–58.4)||11.0 (0.0–33.0)|
|Subaphylla bush%||50||1.5 (0.0–24.4)||1.0 (0.0–8.8)|
|Sub-shrub%||50||0.7 (0.0–8.1)||0.8 (0.0–12.4)|
|Suffrutex%||50||22.8 (0.0–52.1)||7.8 (0.0–31.2)|
|Deciduous liana%||50||0.3 (0.0–2.4)||0.9 (0.0–6.8)|
|Semi-deciduous liana%||50||2.0 (0.0–17.2)||3.7 (0.0–17.2)|
|Perennial liana%||50||0.4 (0.0–6.6)||0.8 (0.0–6.4)|
|Perennial equisetum%||50||0.0||0.1 (0.0–2.3)|
|Semi-perennial fern%||50||2.0 (0.0–24.8)||2.7 (0.0–25.9)|
|Perennial fern%||50||0.1 (0.0–1.8)||0.4 (0.0–5.5)|
|Annual grass%||50||19.8 (0.0–85.1)||24.8 (0.0–79.1)|
|Biennial grass%||50||0.4 (0.0–3.2)||0.2 (0.0–1.5)|
|Semi-perennial grass%||50||1.8 (0.0–15.3)||4.2 (0.0–33.8)|
|Perennial grass%||50||29.7 (1.2–59.9)||34.8 (9.5–63.5)|
|Distance to forest edge||251 (3–1,000)||169 (7–1,000)|
|Rock surface exposure%||5||30 (1–70)||14 (0–40)|
|Soil surface exposure%||5||23 (0–75)||6 (0–30)|
|Number of stones||2.5||6 (0–15)||5 (0–22)|
|Number of stumps||2.5||0.7 (0.0–12.0)||0.1 (0.0–1.0)|
|Number of trees DBH > 0.75 m||2.5||0.1 (0.0–1.0)||0.2 (0.0–1.0)|
|Canopy cover||0.5||1.6 (0.0–20.0)||1.3 (0.0–20.0)|
|Wood species density||0.5||1.1 (0.0–3.0)||0.6 (0.0–2.0)|
|Wood stem density 0.2 m||0.5||0.7 (0.0–4.0)||0.7 (0.0–5.0)|
|Wood stem density 0.4 m||0.5||0.6 (0.0–9.0)||0.2 (0.0–2.0)|
|Wood stem density 0.6 m||0.5||0.2 (0.0–2.0)||0.1 (0.0–1-0)|
|Wood stem density 0.8 m||0.5||0.4 (0.0–6.0)||0.0|
|Grass species density||0.5||3.3 (1.0–7.0)||3.7 (2.0–6.0)|
|Grass stem density 0.2 m||0.5||15.8 (2.0–60.0)||36.3 (10.0–65.0)|
|Grass stem density 0.4 m||0.5||0.5 (0.0–10.0)||1.0 (0.0–7.0)|
|Grass stem density 0.6 m||0.5||0.2 (0.0–5.0)||1.8 (0.0–20.0)|
|Grass stem density 0.8 m||0.5||0.2 (0.0–2.0)||0.04 (0.0–1.0)|
|Litter depth (cm)||0.1||0.2 (0.0–1.0)||0.1 (0.0–1.0)|
The ordination plot obtained by PCA based on the first three axes (explained variance = 28.05%) showed that a large part of the 0.95 confidence ellipsoid generated for C. striatus was included within that generated for C. bedriagai (Fig. 3). The DistLM analysis showed that at the macroscale, the aridity index negatively influenced the occurrence of C. bedriagai (Table 2), where it was the most important variable for explaining the differences between the two skinks (explained variance, EV = 10.6%; Table 3). At the mesoscale level, the proportion of suffrutex and needle-leaved trees had a positive relationship with the presence of C. bedriagai whereas the proportion of deciduous broad-leaved bushes had a negative relationship (Table 2). The best model at the mesoscale level included the proportion of suffrutices (EV = 26.2%) and deciduous broad-leaved bushes (EV = 4.0%). At the microscale level, rock and soil surface exposure had a positive relationship with the presence of C. bedriagai whereas the grass stem density (height = 0.2 m) had a negative relationship (Table 2). The best explanatory model at the microscale level included the grass stem density (height = 0.2 m) (EV = 27.9%) and rock (EV = 8.3%) and soil surface exposure (EV = 9.6%).
|Mean annual temperature||+||0.20||0.645||0.004|
|50 m||Needle-leaved tree%||+||5.66||0.018||0.107|
|Deciduous broad-leaved tree%||−||0.92||0.396||0.019|
|Perennial broad-leaved tree%||−||2.05||0.163||0.042|
|Deciduous broad-leaved bush%||−||6.48||0.010||0.121|
|Perennial broad-leaved bush%||+||1.02||0.322||0.021|
|5–0.12 m||Distance to forest edge||+||0.72||0.408||0.015|
|Rock surface exposure%||+||9.70||0.004||0.171|
|Soil surface exposure%||+||15.29||0.0002||0.246|
|Number of stones||+||0.64||0.438||0.013|
|Number of stumps||+||1.59||0.231||0.033|
|Number of trees DBH > 0.75 m||−||0.68||0.470||0.014|
|Wood species density||+||3.56||0.071||0.070|
|Wood stem density 0.2 m||+||0.03||0.906||0.001|
|Wood stem density 0.4 m||+||1.24||0.382||0.026|
|Wood stem density 0.6 m||+||0.47||0.742||0.010|
|Wood stem density 0.8 m||+||2.09||0.106||0.042|
|Grass species density||−||1.00||0.350||0.021|
|Grass stem density 0.2 m||−||18.22||0.0001||0.279|
|Grass stem density 0.4 m||−||0.65||0.487||0.014|
|Grass stem density 0.6 m||−||2.78||0.090||0.056|
|Grass stem density 0.8 m||+||1.61||0.364||0.033|
±indicates the sense (positive or negative) of the association between the species presence and the environmental predictor. The models were generated independently for each spatial scale. Significant relationships are marked in bold.
proportion of explained variance
|1,000–250 m||Aridity index||5.56||0.021||0.106|
|Deciduous broad-leaved bush%||2.64||0.111||0.040|
|5–0.12 m||Grass stem density 0.2 m||18.22||0.0001||0.279|
|Rock surface exposure%||5.96||0.019||0.083|
|Soil surface exposure%||7.92||0.007||0.096|
The models were generated independently for each spatial scale.
proportion of explained variance
|Venta del Aire||40.11||−0.71||914||CBE||0.83||0.16|
|Saldaña de Burgos||42.25||−3.67||887||CST||0.57||0.42|
|La Vall d’en Bas||42.16||2.39||886||CST||0.10||0.89|
|Losar de la Vera||40.13||−5.60||643||CBE||0.43||0.56|
|Losar de la Vera||40.12||−5.60||510||CST||0.18||0.81|
|Losar de la Vera||40.11||−5.60||479||CST||0.24||0.75|
|Cabo de Gata||36.72||−2.16||174||CBE||0.95||0.04|
|Sant Carles de la Ràpita||40.63||0.56||147||CBE||0.68||0.31|
|Cap de Creus||42.25||3.22||18||CST||0.36||0.63|
Datum for all coordinates, WGS84. Elevation, m asl.
The performance of the random forests model was good (AUC = 0.875). The predictions for the species and sites are shown in Table 4. The model erroneously classified 26% of the C. bedriagai sites and 17% of the C. striatus sites. The false negatives (predicted probability <0.5) for C. bedriagai included localities with a moderate/high herbaceous density (48.1%–87.9%) but that also included ombrophilous taxa (e.g., ferns). All of the syntopic sites were false negatives for C. bedriagai. The false negatives for C. striatus were sites with low/moderate herbaceous densities (30.4%–49.9%). At seven sites, the model predicted similar probabilities for both species (Table 4).
In this study, we investigated the habitat use by two species of reptiles in a broad geographical framework, where we considered the spatial scales that could hierarchically structure their occurrence. Our results showed that the inclusion of macro-, meso-, and microhabitat descriptors could be useful for understanding the complex patterns of niche occupancy by lizards.
These skinks occupied a wide environmental range, where they occurred in high mountains and coastal habitats. At the macroscale level, aridity was the factor that affected species distributions most clearly. Thus, at one of the extremes of the environmental gradient, C. bedriagai occupied semi-desert habitats in the coastal areas of south-eastern Spain (Pollo, 2004b). At the opposite extreme, C. striatus occupied the margins of the temperate-oceanic climate belt in regions with high precipitation (Pollo, 2004a). However, between these extremes, our data showed that both species overlapped in a wide range of ombroclimates. Our analyses also indicated that the characteristics of the substrate did not differentiate the sites occupied by these species.
At the meso-habitat level, we found that the habitats occupied by these skinks also differed in terms of their relative plant compositions. In the habitats occupied by C. bedriagai, the plant associations included higher proportions of suffrutex and needle-leaved (= conifers) trees than those where C. striatus was found. By contrast, in the localities occupied by C. striatus, the plant associations included a higher proportion of deciduous shrubs. Nevertheless, these differences could only reflect clinal variations in the vegetation as a response to the aridity gradient. Suffrutex (e.g., Helichrysum, Lavandula, Thymus) and conifer (e.g., Pinus halepensis) type plants dominated the arid landscapes of this region (Flora Iberica, 2017). Similarly, most of the deciduous shrubs corresponded to ombrophilous genera (e.g., Coriaria, Ligustrum, Salix, Vitex; Flora Iberica, 2017). At the microhabitat level, the sites occupied by C. bedriagai were characterized by higher substrate exposure and lower herbaceous density.
Overall, these results suggest that the potential overlap between the two species is not limited to a narrow ecotonal strip, which is consistent with the syntopy found in this study in geographically distant areas located at different altitudinal levels (El Torcal de Antequera, Málaga, at 1214 m; Monchique, Algarve, at 777 m; Louro, La Coruña, at 6 m). The syntopy sites were characterized by highly suitable conditions for C. striatus but they could be occupied opportunistically by C. bedriagai. However, our surveys showed that C. striatus could also occupy habitats with relatively low herbaceous cover where the models predicted highly suitable conditions for C. bedriagai. Therefore, syntopy could be frequent but it might not always be detected due to the elusiveness of these species (Pollo, 2004a; Pollo, 2004b).
The syntopy between these species could be attributable to their thermal requirements. Chalcides bedriagai and C. striatus exhibit helio- and thigmothermal regulation, and thus both must occur in open habitats. However, within the same habitat, the more specialized morphology of C. striatus could favor this species occupying areas with greater herbaceous cover whereas C. bedriagai could occupy the habitat randomly. This hypothesis was supported by our results at this spatial level. The differences in microhabitat use could reduce mutual interactions (because they predate mainly upon small arthropods; Salvador, 1998) and favor syntopy between the two species, which is relatively unusual in congeneric Mediterranean lizards (Arnold, 1987).
The spatial overlapping patterns could be also explained by the high genetic divergence between both species (25.1 Ma; Zheng & Wiens, 2016). It is suggested that complementary niches only appear between distantly related species because unrelated species are less likely to use the same type of resources (Losos et al., 2003; Srivastava et al., 2012). In this sense, there are also examples of sibling species within the genus Chalcides (C. chalcides/ C. striatus, estimated divergence 10.1 Ma; Zheng & Wiens, 2016) showing parapatric distribution patterns where they meet (Sindaco et al., 2006).