Connectivity Predicts Presence but Not Population Density in the Habitat-Specific Mountain Lizard Iberolacerta martinezricai

Abstract

The Batuecan lizard Iberolacerta martinezricai is a critically endangered species due to its significantly reduced distribution, which is restricted to the scree slopes (SS) of a few mountain peaks within the Batuecas-Sierra de Francia Natural Park (western Spain). Given its high specialisation in this type of discontinuous habitat, the long-term conservation of the species requires maintaining the connectivity between populations. This study analyses the contribution of connectivity, as well as other patch-related factors, in the distribution and density patterns of the species. With this aim, 67 SS were sampled by line transects from May to October 2018. Each SS was characterised using variables indicative of the microhabitat conditions for the lizard. Inter-SS connectivity was quantified using graph theory for seven distances. Generalised linear models (GLMs) were performed for both presence and density. Model results showed that while connectivity was a relevant factor in the presence of lizards, density only involved patch-related variables. Discrepancies probably occurred because the factors influencing presence operate on a wider scale than those of abundance. In view of the results, the best-connected SS, but also those where the lizard is most abundant and from which more dispersed individuals are likely to depart, seem to be the essential patches in any conservation strategy. The results may also be relevant to other species with habitat-specific requirements.

1. Introduction

Although climate warming poses a major threat to the terrestrial vertebrate species confined to high mountain areas, not all species are equally vulnerable [1,2,3]. Ectotherms are particularly sensitive because they depend on external heat sources to maintain body temperature [4]. In these species, an increase in temperature is likely to induce fast growth and thermal stress, which can accelerate the rate of ageing [5]. To cope with the variation in environmental conditions, species can either try to adapt their physiology or behaviour to extend their thermal niche to these new conditions, or move to new areas where their thermal niche is now located [6,7]. The first option, acclimatisation, seems limited as mountain species tend to be cold-specialists [8], experiencing a greater decline in fitness if body temperatures exceed the optimum [9,10,11,12]. The second one, elevational shifts, are not an option for mountaintop species, as there are no colder areas available to migrate to [13,14]. In addition, mountain species may compete with and be displaced by generalist species that are present in the surrounding lowlands and that may move into higher altitudes with increasing temperature [1,15]. All of these factors make mountain ectotherms, such as the Batuecan lizard Iberolacerta martinezricai, among the most vulnerable animals in the world [2].

Iberolacerta martinezricai is a medium-sized, globally threatened rock lizard of very reduced spatial distribution that lives in less than 20 km² of discontinuous habitat within the limits of the Batuecas-Sierra de Francia Natural Park (Castile and Leon, western Spain) [16]. I. martinezricai is the most restricted, and probably one of the rarest, threatened reptile species in continental Europe.

Due to the projected scenarios of climate warming and drought for the Mediterranean mountain region [17,18], the predicted changes in the altitudinal and latitudinal distribution for the Spanish herpetofauna indicate a significant decline in suitable habitat for the species [19,20,21]. In addition, two other Mediterranean lizard species are increasingly found near the summit of the Peña de Francia mountain: Podarcis guadarramae and Timon lepidus. The presence of these lizards could indicate the advance of Mediterranean forms toward the summit of the mountain, corresponding to a progressive cornering of the Batuecan lizard in the highest areas by competition and climatic causes [22].

Nevertheless, climate warming and probable interspecific competition are not the only threats that the Batuecan lizard faces. Although the perceived remoteness and inaccessibility of these mountain areas give them some natural protection from anthropogenic impacts, human pressure on these ecosystems has increased markedly in recent decades [23]. Thus, recreational activities in mountains are associated with reduced vertebrate richness and abundance [24,25], particularly in lizards [26,27]. In the case of I. martinezricai, within its area of distribution there is a mountain monastery frequented by tourists, which makes human pressure in the area relevant. This impact is particularly important during summer months, which coincides with the active period of this species [22]. The dirt track to the monastery was constructed in 1920 and asphalted in 1961, facilitating massive access for motor vehicles. The highly visible, characteristic antenna was installed in the 1970s. These landmarks show the increasing human pressure upon the area over the past hundred years [28].

Access roads to the summit and unpaved roads in the vicinity may also affect the species. Even though small lizards seem not to be the group most affected by roadkill [29,30], other impacts of the roads, such as the barrier effect or fragmentation, can influence their populations [31]. Moreover, the asphalting of the road at the summit caused, at least, necrosis of the fingers of adult specimens, which could be due to sunbathing in contact with the pasty asphalt. These specimens had asphalt-coated toes and had lost parts of those toes, with naked phalanges poking out [32]. Human presence in the area can also favour the occurrence of wildfires [33], which can irremediably change the microclimate of the areas concerned.

All of these pressures and threats to the species are compounded by the fact that I. martinezricai is a strictly saxicolous lizard specialising in a particular type of discontinuous and naturally fragmented habitat, such as scree slopes (SS), which are interspersed with grass and broom mountain vegetation [16,22]. Although the Batuecan lizard has evolved in this naturally fragmented habitat, current human pressures may compromise the dynamics necessary to maintain its viability in the medium and long term. The negative effects of fragmentation on lizards have been widely documented. Fragmentation can alter home ranges and movements [34], demographic structure [35], or community disassembly [36]. Not all species are equally susceptible, and those with the strongest habitat preferences show the greatest response to landscape fragmentation [35,37,38].

Thus, as it is a lizard highly specialised in a specific habitat and that most probably responds to a classical model of metapopulation dynamics, it is crucial to study the importance of the connectivity between the favourable patches for the species. Previous studies have shown that some patch-related variables shape the distribution and the abundance of lizard populations in the SS [16], but the spatial structure of the SS network and the importance of each SS within this network have not been evaluated. With this objective, this study aims to assess whether connectivity between SS, together with other intra-SS factors, contributes to defining the patterns of presence and density of the Batuecan lizard in the Batuecas-Peña de Francia Natural Park. It is hypothesised that connectivity will be more related to presence than to lizard density, as the latter is more likely to be closely linked to intra-patch availability of resources. If such discrepancies exist, it is likely that certain regionalisation within the species’ own range will emerge. Understanding these dynamics will help guide management and conservation plans for this threatened species.

2. Materials and Methods

2.1. Study Area

The study area is in the Las Batuecas-Sierra de Francia Natural Park (BSFNP), Province of Salamanca (Central Spain). BSFNP has been recognised as a biosphere reserve by UNESCO as well as a Special Protection Area and Special Area of Conservation within the Natura 2000 Network. BSFNP is part of the western area of the Sistema Central, with the Pico Hastiala (1735 m) and the crest of the Peña de Francia (1723 m) as its highest peaks (Figure 1). The area is characterised by a Mediterranean mountain climate, with hot dry summers and cold winters in which, along with in spring, the maximum rainfall occurs [39].

Regarding the geology, the area is dominated by Palaeozoic metamorphic materials: quartzites, schists, and slates. The differential erosion of these materials has formed a rugged relief with peaks comprising the hardest Armorican quartzite (Hastiala, Peña de Francia, Mesa del Francés, Rongiero), and deep river valleys in the areas with less resistant rocks. The peaks (hanging synclines) form the watershed divide between the basins of the Tagus and the Duero, with a marked contrast between the northern and southern parts due to the difference in altitude between the Castilian and Extremaduran plateaus. The slopes of the quartzite crests are covered by colluvium and screes that have been partly phytostabilised because of processes that are still active and linked to a periglacial morphogenetic context. These SS are the characteristic habitat of the I. martinezricai lizard.

Due to the differences in altitude, up to three bioclimatic levels are represented, each with a characteristic climate and vegetation, The Oro-Mediterranean forms the highest areas (altitude above 1600 m.a.s.l.), where mountainous shrub vegetation dominates, mainly Cytisus oromediterraneus. The Supra-Mediterranean level (1000–1600 m.a.s.l.) is where the potential vegetation comprises Pyrenean oak Quercus pyrenaica forests and a shrub stratum of heather Erica arborea and Erica australis; some scattered plantations of Pinus nigra and Pinus sylvestris are also present. The lower parts belong to the Meso-Mediterranean level (400–1000 m.a.s.l.), which is characterised by wood pastures of holm oak Quercus ilex and cork oaks Quercus suber.

2.2. Focal Species

The Batuecan lizard I. martinezricai is one of the 8 species of the genus distributed in different mountain ranges of the Iberian Peninsula and the Alps. The speciation of the genus Iberolacerta Arribas, 1997 is likely to have occurred by allopathic divergence during the late Tertiary because the complex relief of southern Europe provided isolated habitats with limited possibilities of dispersion between them [40]. This may have been helped by competition from wall lizards (genus Podarcis) present in the lower areas [41]. Mountains acted as refuges during the glacial oscillations of the Quaternary [40]. Genetic analyses show that the Batuecan lizard has been present in the mountainous area it currently inhabits for at least 2 ± 0.3 Ma [41].

The area of distribution for I. martinezricai is restricted to an area of 20 km² of discontinuous habitat within the boundaries of the Batuecas-Sierra de Francia Natural Park. Specifically, the lizard is contained in a polygon formed by the peaks of Peña de Francia, Hastiala, Rongiero, and Peña Orconera, on whose steep SS the habitat available for the species is concentrated [16,28]. Although it can be found from 800 m.a.s.l., it mainly occupies medium and high altitudes, above 1300 m.a.s.l. These mountains are surrounded by low-altitude areas with unsuitable habitats for the lizard. In terms of abundance, maximum densities of 50 individuals per hectare are reached [16], although this is lower on most SS. It prefers habitats with a certain amount of humidity and with the presence of mosses and lichens [16,28]. The total number of specimens for the species has been estimated to be between 1200 and 1500 individuals [16], distributed in several fragmented populations with a different degree of connectivity between them. Due to this reduced distribution, it has considered critically endangered since 2006 (CR B2ab(v), C2a(ii) in the IUCN Red List) [32,42,43]. Since 2019, it is also listed as Endangered in the Spanish National Catalogue of Threatened Species. Thus, the species is one of the European reptiles in most urgent need of conservation [40].

The species is a medium-sized lizard with a dark dorsal colouration, usually brown or grey with a black, sometimes greenish, reticulate, and bright blue axillary ocelli (between one and three) in most individuals. Sexual dimorphism is less pronounced than in other rock lizards. Sexual maturation is probably reached when the lizards are three years old, and the longest living recorded specimens were 8 years old [44]. The temperature range within which the species remains active is relatively wide (22.8–39.2 °C). It is a good thermoregulator in relation to substrate temperature, but not so good with respect to air temperature [28]. Regarding reproduction, copulation begins at the end of April and oviposition takes place during the summer months (especially at the end of July). The size of the clutches ranges from two to six eggs, with an average of 4.23. The incubation period lasts between 33 and 42 days [45]. There are still no data on species movement, its capacity to disperse, or its willingness to cross different habitats of the landscape matrix in which the SS are located.

2.3. Field Surveys and Density Estimation

Sampling was carried out from May to October 2018, covering the activity period of the lizard. During the spring, the Batuecan lizard shows a unimodal pattern of activity. In summer, the pattern is more bimodal, with a main peak during the first hours of sunlight and a marked reduction afterwards as the temperature increases. Evening activity is much less [28]. Patches of the characteristic habitat of these lizards (mountain SS) were positively selected as the preferred study sites. Less favourable areas for the species were a priori excluded from the surveys according to our current knowledge of its habitat preferences, namely low-elevation areas with warmer temperatures. At each SS, surveys were performed between 08:00 h and 19:00 h GMT along linear transects of about 100 m [46] (average of 96.89 m; differences are due to accessibility difficulties in some SS, given the abrupt topography of the terrain). At least one transect per 10 hectares of scree surface was carried out in each SS (average of 1 transect per 3.58 ha). For SS with more than one transect, sampling was conducted on the same day and transects were separated by at least 50 m to avoid double counting. All individuals seen in a 3 metre wide band on either side of the transect were recorded. The distance at which individuals were observed was measured perpendicularly to the transect line with a Bosch DLE50 Professional laser distance meter. Whenever possible, the sex and age of the individuals reported were identified.

Distance 6.0 software [47,48,49] was used to estimate the density and abundance of lizards from the observations collected during the fieldwork. The number of data were sufficient to achieve an adequate adjustment of the detection function [50]. Using the lowest value of the Akaike Information Criterion and the smallest proportion of the χ² goodness-of-fit divided by its degree of freedom, the detection function model that best fit the study data was selected [48,49]. In this case, a half-normal key function with cosine adjustment was adopted.

2.4. Connectivity Analysis for I. martinezricai

Spatial connectivity is defined as the ability of the landscape to facilitate or impede the movement of individuals of a species between fragments of habitat. It is usually divided into two components: structural and functional [51]. The structural component is determined by the spatial connection of different types of habitat in the landscape, and the functional component refers to the response of individuals and species to the physical structure of the landscape. The latter is influenced by the species’ habitat requirements, tolerance of altered habitats, and life stage. In this sense, species, although living in the same habitat, have different behavioural responses and therefore experience different levels of connectivity [52].

The study of ecological networks and connectivity is often carried out using graph theory [53,54,55]. In terms of landscape ecology, a network consists of a set of favourable habitat patches, the nodes (in this case the SS, the habitat of I. martinericai), embedded in a non-favourable habitat matrix [56]. If any ecological flow occurs between any two nodes of the network, both nodes are said to be “connected” [57]. These analyses are useful for the design of conservation strategies [58,59].

For the connectivity analysis of the metapopulations of lizards within the SS network, the probability of connectivity index (PC), based on this graph theory, was used [60]. As the Batuecan lizard is strongly saxicolous and SS are discrete geographical units, the definition of the nodes was relatively simple. The analyses were carried out using the CONEFOR software [61]. This index describes the probability that two organisms situated at random in the landscape are in patches (in this case, SS) that are interconnected. Therefore, PC indicates the importance of each SS within the network. The mathematical expression for PC is:

$$PC = \frac{{\sum }_{i = 1}^n{\sum }_{j = 1}^n{a}_i{a}_j \times p_{ij}^{*}}{A_L^2} = \frac{P{C}_{num}}{A_L^2}$$

in which n, in this case, is the number of SS in the study zone (not only those sampled in the study, but all of the mapped SS, a total of 156); a_i, a_j are the areas of the SS i and j; A_L is the area of the study zone, including the SS and the landscape matrix located between them; and p*_ij is the maximum probability of dispersion between i and j, which is determined by the following negative exponential function [60]:

$$p_{ij} = e^{-k{d}_{ij}}$$

where the p_ij value is usually calculated from d_ij distances obtained from mark-recapture data or genetic analysis. The software adjusts the previous function for these values, obtaining a constant k, from which the dispersion probabilities are calculated and, therefore, the connectivity probability for each SS. In this case, as there is no information on the capacity of dispersion and mobility between SS, connectivity was simulated for different distances. Thus, connectivity probabilities of 50% for seven distances were tested: 25, 50, 100, 250, 500, 1000 and 10,000 (practically unlimited connectivity) metres. Considering the characteristics of the species, we believe that this distance interval includes both the most frequent events of short distance dispersal and events of colonisation of distant areas.

Saura and Rubio [62] showed that PC can be divided into three fractions quantifying the different ways in which a habitat patch (in this case, a SS) can contribute to connectivity: dPC_intra shows the contribution of a given SS in terms of intrapatch connectivity, and it is related to the habitat availability; dPC_flux corresponds to the attribute-weighted dispersion flow (usually area-weighted) through the maximum-likelihood paths of a given SS with all other SS in the network when it is the point of origin or destination of that flow; and dPC_connector for a given SS quantifies the contribution of that SS to the connectivity between the rest of the network. A SS can only contribute to this fraction if it is part of the maximum-likelihood paths between SS other than itself. dPC_flux and dPC_connector are measured in the same units and can be directly compared to each other.

All of the potential SS habitats (not only those sampled) for the species in its small natural range were included in the connectivity analysis. SS cartography was obtained from the Spanish Land Occupancy Information System (SIOSE) developed by the National Geographic Institute (IGN), which is the most detailed database of land occupation for all of Spain, at a reference scale 1:25,000. Although SIOSE shows the main SS in the area with remarkable precision, some tesserae are mapped as mixed units that include smaller SS, which are probably also interesting for the species. The well-defined SS in the orthophotos within these mixed units were mapped and joined with the rest of the SS for analysis.

2.5. Presence and Density Models for I. martinezricai

A species inhabits those areas where abiotic conditions are favourable and where the community of species allows its coexistence, and in accessible places that could be colonised in both evolutionary and ecological terms [63]. All of these factors interact dynamically to produce a geographic range that is unique to each species [64]. These distribution ranges may shift, expand, or contract in response to changes in environmental conditions [65]. On the other hand, population density is also related to environmental characteristics, individual traits, and intra- and interspecific interactions [66,67]. However, the factors involved and/or their intensity do not necessarily have to be the same to determine presence and abundance.

Together with the probability of connectivity index (PC), the presence and density models included as explanatory variables those already identified as relevant for the species in previous studies [16]. Although climatic data from sources such as WorldCLIM2 [68] are often used in presence/abundance models, the resolution of these sources of information is not valid to account for the different microclimates that contribute to explaining the specific and reduced geographical range of I. martinezricai. For the same reason, climatic data from one or more nearby weather stations would not have been useful either. To overcome these limitations, a set of variables related to topography that could well reveal microclimatic aspects relevant to the species presence or density was selected (altitude ALT and orientation ORI). As it seems that the species positively selects those SS under which there is more humidity, the Topographic Humidity Index (TWI) was introduced into the analyses. TWI describes the tendency of a “cell” to accumulate surface water [69]. The abundance of lichens and bryophytes in SS rocks (LIC) (related to humidity levels, the slope (SLP), and rock size (ROS)) were also considered. The description of all of the variables considered is included in

Table 1. Description of the variables used in the presence and density models for I. martinezricai in the Batuecas-Sierra de Francia Natural Park.

Variable	Abbrev.	Description
Altitude	ALT	Obtained from a 5 metre resolution DEM, developed by the Spanish National Geographic Institute through interpolation from LIDAR data. The minimum (ALTmin), average (ALTave), maximum elevation (ALTmax), and altitudinal range (ALTran) were considered.
Slope	SLP	Obtained from the same DEM, expressed in degrees.
Orientation	ORI	From the 5 metre DEM using the ArcGIS 10.6 “aspect” tool. Five orientations were defined: North (N), South (S), East (E), West (W), and flat.
Probability of connectivity	PC	An indicator of landscape connectivity based on habitat availability, dispersal probabilities between habitat patches and graph structures (see Section Section 2.4). Obtained from SS cartography and CONEFOR software.
Topographic wetness index	TWI	Describes the tendency of an area to accumulate water. Obtained from the 5 m DEM using ArcGIS geoprocesses.
Lichen cover	LIC	Average coverage of lichens and bryophytes on the rocks.
Rock size	ROS	Each SS was classified in relation to the average size of the rock fragments. Three categories were considered: 1 (<50 cm), 2 (50–100 cm), 3 (>100 cm).

.

To investigate the relative importance of each variable in species presence and density, generalised linear models (GLMs) were performed using R [70]. A preliminary statistical comparison between SS with and without the lizard was carried out using non-parametric Mann–Whitney U tests for the quantitative continuous variables and the Pearson Chi-Square test for the categorical one (ORI). GLM with a logit-link function and binomial distribution were used to model presence/absence data and GLM with a log-link function gamma distribution to model species density. To avoid introducing noise into the GLMs, potential explanatory variables that did not provide information were previously eliminated. Collinearity among explanatory variables was explored using tolerance levels, in addition to Goodman and Kruskal’s gamma, to test for categorical variables against other variables. Based on this, the variables that were highly correlated were removed. Tolerance levels were sufficiently high (>0.1), and Goodman and Kruskal’s gamma coefficients were low enough (<0.4) to allow the inclusion of the variables included here.

For both presence and density, a full model with all predictors was tested first. We then removed variables that were clearly insignificant. To select from the multiple models, the Akaike’s information criterion (AICc) was used as a measure of overall model fit [71]. For the best fit GLM models (∆ICc = 0), the calc.relipm procedure was implemented in R-package “relaimpo” [72] to decompose the variance of the final models among the different predictors and interactions.

3. Results

Sampling Results

A total of 67 SS were sampled, with the species present in 30 of them (44.7%). The total number of individuals observed was 154, including adults (49%, 23% males and 26% females), sub-adults (22%), and juveniles (13%). No proper identification could be made in 16% of the cases. In general, the specimens were detected when the temperature was not very high (around 22–25 °C) and in the early morning hours. The average density of lizards for the SS with confirmed species presence was 41.44 ± 27.81 individuals/ha, with a maximum of 99.39 individuals/ha. Figure 2 shows the SS sampled and those where the presence of the species was detected. Table S1 (see Supplementary Material) shows the results for each SS including area, number of transects, mean length, number of individuals located, and density.

The results of the Mann–Whitney U and Pearson Chi-Square tests are shown in

Table 2. Average values in the SS with and without the presence of the species and the results of the Mann–Whitney U (quantitative variables) and Pearson Chi-Square (categorical variables) tests comparing both groups. The significance level was 0.05. Statistically significant differences are shown in bold.

	SS Positive (Average)	SS Negative (Average)	U Mann Whitney	Z	p-Value
ALTmin	1330.1	1180.5	241.5	−3.612	<0.001
ALTmax	1510.7	1344.8	189	−4.319	<0.001
ALTave	1421.6	1263.9	198	−4.197	<0.001
ALTran	180.7	164.3	463	−0.632	0.527
SLP	24.77	26.50	340	−2.287	0.022
TWI	9.41	9.03	408	−1.372	0.170
PC10000	2.05	1.13	318	−2.583	0.010
PC1000	3.46	1.31	264	−3.310	0.001
PC500	3.61	1.30	268	−3.256	0.001
PC250	3.43	1.26	287	−3.001	0.003
PC100	3.39	1.29	228	−3.311	0.001
PC50	3.53	1.26	253	−3.331	0.001
PC25	3.49	1.29	228	−3.211	0.001
PCintra *	0.13	0.06	318	−2.583	0.010
PCflux *	2.50	0.92	262	−3.336	0.001
PCconnector *	0.98	0.32	327	−2.464	0.014
LIC			275.5	−2.542	0.011
ROS			205.5	−3.334	0.001
			χ2	d.f.	p-value
ORI			3.871	3	0.276

* The PC_intra, PC_flux, PC_connector values shown in the table are those obtained for the PC₅₀, the distance for which a better adjustment of the presence models was obtained (see Table 3).

. The SS with lizards present showed a higher maximum, average and minimum altitude, lower slope, and higher PC (at all distances considered) than SS without the species. The three elements that compound the PC also showed significant differences. Rock size and lichen coverage were greater in the positive SS. No differences were found for the range of altitudes and TWI (Table 1), nor for the orientation according to the Pearson Chi-Square test.

The results of presence and density models are shown in

Table 3. Ranking of generalised linear models for Batuecan lizard presence and density, using Akaike’s Information Criterion corrected for small sample sizes (AICc). Total deviance explained (%∑dev) for the best models, and % explained by the main factor are also shown.

Candidate Presence Models	∆AICc	%∑dev	Candidate Density Models	∆AICc	%∑dev
ALTave+PC50+ROS+SLP	0	48.29	LIC+ALTave	0	30.63
ALTave+PC50+ROS+SLP+LIC	1.63	ALTave (38.1)	LIC	0.50	LIC (85.7)
ALTave+PC50+ROS	3.64		LIC+ALTave+PC25	2.00
ALTave+PC50	10.35		LIC+ALTave+PC25+SLP	3.50
ALTave	20.46		LIC+ALTave+PC25+SLP+ROS	5.50
PC50	22.91		ALTave	7.00
PC25	24.05		PC25	7.00
PC100	24.29		SLP	7.10
PC1000	25.85		ORI	7.10
PC500	26.12		PC500	7.40
PC250	26.54		PC250	7.40
PC10000	26.87		ROS	7.40
ROS	27.11		PC10000	7.50
SLP	33.14		PC100	7.50
LIC	34.33		PC1000	7.60
ORI	38.90		PC50	7.60

. The most parsimonious models included four variables (ALTave, PC₅₀, ROS, and SLP) for the presence model and 2 (LIC and ALTave) for the density model. The best presence model showed a total deviance of 48.29%. ALTave contributed 38.1%, PC₅₀ 24.4%, ROS 27.8%, and SLP 9.7%. The total deviance explained by the best density model was 30.63%. LIC contributed 85.7% and ALTave 14.2%.

4. Discussion

According to the results, it seems that the species presence involves variables related to connectivity as well as SS characteristics. Species density is only influenced by the microhabitat conditions present in each SS. Thus, the probability of lizard presence in an SS is correlated with the importance of that SS for the connectivity of the SS network, higher altitudes, less pronounced slopes, and large rock sizes. The variable that explained the highest percentage of the deviance was the altitude of the SS. A large rock size favours the presence of cavities that serve as refuges for the lizards, storing fresh air and generating a more humid environment, probably with a greater density of prey [16]. As mass movements are gravitational phenomena, less pronounced slopes favour larger rock sizes in the SS. The probability of occurrence for the species is minimal for SS in low and thermophilic areas, far from and poorly connected to the distribution centre. On the contrary, lizard density depends not on the position of the SS within the network but on the characteristics of each SS. The populations of I. martinezricai are more abundant in the SS with greater coverage of lichens and mosses (largely responsible for the explained deviance), indicators of greater humidity and possibly greater availability of prey. This abundance and presence respond to different factors, or at least there are variations in the weight of their contribution, which has already been described for other species of lizards [73]. The discrepancies are probably due to the fact that the factors influencing presence operate on a wider scale than those of abundance [74]. This seems to be the case for the Batuecan lizard, where the presence is shaped by both inter-SS spatial relations and intra-SS characteristics, and density only by the latter. However, it depends on the species, since, for example, in forest lizards, it has been found that habitat quality is also the most relevant variable in their spatial distribution [75].

The importance of the connectivity indicates that the Batuecan lizard probably acts as a metapopulation in which each SS constitutes a small subpopulation, reaching an equilibrium over time by compensating for local extinctions through recolonisations from other SS [76]. The best fit of the models was obtained for PCs with small dispersal distances, which could indicate a limited capacity of the species to reach distant SS. This limited capacity, together with the importance of intra-SS characteristics in lizard density, is likely to lead to the formation of clusters of interacting individuals, called neighbourhoods, that are regionally organised into continuous networks within the total range of the species [77,78]. In addition, connectivity was shown to be important for all distances tested, which may be indicating population flows operating at different scales. Structural connectivity shapes the spatial distribution of the lizard. The three fractions in which PC can be divided showed differences between SS with and without lizard presence. This means that both habitat availability and dispersal flows are important for species presence in a given SS.

The results for lizard presence are maintained over time, since the current distribution is very similar to that obtained 10 years ago [16], which shows a remarkable stability. The scarce changes between the two periods are mainly located in SS with low density for the species in the low altitude peripheral areas of its distribution area. These differences could be explained by the difficulty of detection, conditioned by the meteorology of the sampling days. In contrast, the comparison for density does not show such a sustained pattern over time. This may be due to several causes. First, the sampling in this study was carried out using transects rather than plots in order to cover a greater number of SS. Secondly, this is a species present in very low densities, so changes in the detection of one or two individuals, strongly influenced by weather conditions, can lead to significant changes in the densities of lizards obtained. These uncertainties mean that the results for the density models should be taken with caution.

Connectivity is especially relevant in the context of climate change. With increasing temperatures in summer and a reduction in humidity in the SS, it is possible that the species may experience an elevational shift, abandoning the patches located at lower altitudes where habitat conditions seem better for thermophilic generalists [15]. Individuals need to be given the opportunity to reach other, more suitable SS, under new environmental conditions. In addition, mountain species may be preyed on, compete with, and be displaced by generalist species present in the surrounding lowlands that may go to higher altitudes. For example, during the fieldwork, several species that can compete or prey on I. martinezricai were found within the same altitudinal range, such as the smooth snake Coronella austriaca, the Lataste’s Viper Vipera latastei, the Lusitanian lizard Podarcis guadarramae, or the Algerian psammodromus Psammodromus algirus. P. guadarramae is likely to be the main candidate to outcompete or substitute I. martinezricai, as it is also saxicolous and occupies the same habitat. However, at present, there are not enough data available to describe the potential impact of these species on the Batuecan lizard populations.

Furthermore, considering connectivity in the specific case of the Batuecan lizard is important because, unlike other mountain species that inhabit remote areas with low human pressure, the presence of the monastery makes the mountains that constitute the small area of distribution of the lizard very popular for tourism. Although tourists do not usually have access to SS, tourism pressure on the Peña de Francia can affect the quality of habitat for I. martinezricai (i.e., with food remains) and attract opportunistic species that can predate upon lizards (foxes, stone martens, and rats, among others). Perhaps more importantly, this human presence can favour the occurrence of forest fires, which are frequent in the area. The population consequences of fragmentation due to human pressures on different lizard species have been widely documented [36,79,80].

The results, therefore, highlight the need to include, with particular emphasis, connectivity criteria in the management and conservation plans of threatened species that are highly specialised in discontinuous habitats, to maintain their populations in the long term. Active or passive strategies need to be developed to promote accessibility among different SS, ensuring movement of individuals and gene flow between patches [81,82,83]. Connectivity analyses also allow the identification of key patches [84]. However, not only the best-connected SS, but also those where the lizard is most abundant and from which more dispersed individuals are likely to depart, seem essential in any conservation strategy for the species.

This connectivity analysis should be complemented and extended by other detailed studies based on genetic analysis [85] and movement resistance in relation to the different land uses between the SS, quantifying the effect of the configuration and quality of the matrix in gene flow and genetic variation [55,86]. Furthermore, future studies on the effects of climate change, thermal ecology, predation, and competition with other lizards are essential for the design of appropriate management tools for the species.

5. Conclusions

Connectivity, together with patch-related variables, are important factors in explaining the presence of the critically endangered Batuecan lizard I. martinezricai in the SS of the Batuecas-Sierra de Francia Natural Park. In contrast, lizard density in each SS is not related to the spatial configuration of the SS network, but only patch-related variables have an influence, especially those indicative of a higher degree of humidity. Considering the metapopulation dynamics that emerge in these discontinuous habitats, the discrepancies in the variables related to lizard presence and density mean that, in terms of conservation, it is necessary to give special importance to two types of SS. On the one hand, those SS that act as the main connectors within the network favour the recolonisation processes (which counteracts local extinctions). On the other hand, those SS whose internal characteristics favour abundant populations that feed the stochastic processes are involved in the long-term maintenance of all these metapopulation dynamics. Based on these two groups of SS, we would expect to see the formation of “neighbourhoods” of interacting individuals, organised regionally within the total range of the lizard. These spatial dynamics should be taken into account when designing conservation plans for the species. Moreover, the results may also be relevant to other species with habitat-specific requirements.