Analysing How Pre-Fire Habitat Legacy and Post-Fire Management Inﬂuence the Resilience of Reptiles to Fire

: At the landscape scale, the Mediterranean region is a mosaic of habitats occupied by plants and animals with different resilience to ﬁre. One of these habitats, the pine plantation, is characterized by its structural simpliﬁcation and susceptibility to ﬁre. Despite its high ﬂammability, few studies have compared the response of animal communities between pine plantations and other autochthonous woodlands. For ﬁve years after a large ﬁre in southwestern Europe, we surveyed reptiles in two natural habitats (oak forest, scrubland) and a pine plantation managed with salvage logging, a post-ﬁre practice which consists of the complete harvesting and removal of death burnt trees. Reptile abundance and species composition were examined to assess differences in the reptile community between these habitats. Differences between burnt and unburnt transects were limited to the ﬁrst year after the ﬁre, while, over the entire ﬁve-year period, differences in species composition and abundance were due to vegetation type instead of ﬁre. The pine logged area showed a delay in the recovery of vegetation and also in the appearance of many reptile species after the ﬁre. At the reptile species level, we found evidence of both positive responses to ﬁre (for lizards with high heliothermic activity) and negative ones (for specialist snake species). Overall, our results conﬁrm the resilience of the reptile community to ﬁre. The mosaic of habitats in the Mediterranean region and the openness caused by ﬁre can increase the reptile biodiversity (landscape- plus pyro-diversity effects), but some practices such as salvage logging coupled with ﬁre regime shifts (larger and more frequent ﬁres) can compromise the conservation of the biodiversity in ﬁre-prone regions.


Introduction
Fire modifies forest ecosystems [1] and has been a source of biodiversity for millions of years [2]. In the temporal sense, wildfire plays a key role in forest dynamics shaping biodiversity composition from microbial to vertebrate communities along the post-fire succession [3,4]. Spatially, fire regimes may lead to heterogeneous landscapes with a mosaic of burnt and unburnt patches at different times after a fire [5]. In recent decades, wildfires have changed in frequency, extent, and severity in many regions worldwide [6], giving rise to a fire regime shift that is threatening biodiversity [2]. This shift is caused by socioeconomic factors such as rural abandonment, with a subsequent reduction of agricultural patchworks or open wooded spaces that are sometimes replaced by extensive pine plantations, resulting in vast accumulations of combustible material, making the landscape more prone to burning [2,7,8]. Some regions, such as the Mediterranean basin, are more susceptible to blazes as a consequence of longer and hotter drought seasons, induced by anthropogenic climate change [2]. In addition, the expansion of the wildland-urban interface increases the risk of fire ignition [7], and firefighting efforts have diminished, as protecting human private property against fire is given priority over wildfire suppression [9]. Even though Mediterranean vegetation species have strategies to survive wildfires [1,10,11], greater fire frequency and intensity may exceed plant resilience [2,12]. Fire regime shifts can cause ecosystem degradation, a simplification of the plant and animal communities, and increase the extinction risk of threatened species [2].
Some species respond positively to recently burnt areas (open-landscape specialists) whereas, others occur only in long-unburnt areas [13][14][15]. A landscape composed of a mosaic of patches at different time periods since burning is more pyro-biodiverse than a single large space caused by a megafire. Similarly, in a landscape mosaic composed of forest patches with different resilience to fire [16,17], the occurrence of a large fire can enhance biodiversity by the combination of species adapted to unburnt and burnt habitats [18].
In burnt landscapes, human management can influence pyrodiversity [19], and shape the response of animal communities to fire and, consequently, their resilience [20]. While an appropriate post-fire intervention can minimize negative impacts (e.g., pest risk reduction; [21]), the harvesting of all burnt trees can induce unpredictable consequences far worse than simply altering the vegetation by logging [22]. Salvage logging (hereafter SL), i.e., the harvesting of trees after natural disturbances [23], is a common practice after blazes in many fire-prone regions including the Mediterranean basin [24,25]. This practice is conducted to recover economic value from the burnt wood and to reduce pest xylophagous insects. As a consequence of the complete biomass suppression, SL impedes several processes in post-fire ecosystems, such as facilitation for seed recruitment [25,26]. Moreover, SL can facilitate the colonization of invasive species and delay the natural regeneration of the plant community [23,26]. The complete removal of dead branches and trunks also eliminates favorable microenvironments that attract large numbers of arthropods, which are key in maintaining the diversity of other animals and plants. Furthermore, SL involves heavy machinery moving over the terrain, damaging soil quality and reducing shelter for ground-dwelling animals [24]. Consequently, SL alters both animal and plant species [26][27][28] and can remove habitats left by fire [29,30]. However, the effects of SL vary between ecosystems and, consequently, more studies are necessary for numerous taxonomic groups, such as reptile species [23].
Our study focuses on the response of the reptile community after a major wildfire that was partially managed with SL. In general, population recovery after a fire may be owed to in situ survival or to immigration from surrounding unburnt areas [31,32]. Reptiles are normally small organisms with restricted home range, limited colonization capacity, and high territorial fidelity [33,34]. Therefore, their survival could depend more on the existence of shelters against fire, such as trunks, roots, and rocks [35] than on the recolonization from unburnt areas. Accordingly, SL can compromise the recovery of the reptile community by eliminating such shelters, delaying plant regeneration, and thereby hindering the recovery of arthropod populations.
Here, we examine the response of reptiles over a five-year post-fire chronosequence in a burnt area composed of natural vegetation (oak forest and scrubland) and a humandesigned area (pine plantation). In this scenario, we propose the following hypothesis: The resilience of the reptile community is habitat-dependent, some habitats being more resilient than others [7,36], and thus the reptile community in oak forest and scrubland will show more resistance to fire than in pine plantation [37]. Accordingly, the reptile community will rapidly return to the pre-fire situation in the most resilient burnt habitats [35,38]. Resilience in our study was assessed by comparing burnt and unburnt patches within the same vegetation type. Unburnt patches were selected in the surrounding unburnt areas near the fire perimeter. The pine plantation was completely logged after the fire. Despite no other habitats being logged, we examined the occurrence of the reptile community in this burnt habitat in order to understand the effect of fire plus the post-fire management on this group of vertebrate organisms.

Study Area and Post-Fire Management
The study area is located in the western central Iberian Peninsula (Figure 1). In August 2015, a large forest fire burned 7833 ha. The fire spread along a large valley (Rivera de Gata; hereafter Gata valley) open to the south and surrounded by steep slopes of the foothills of the Central System to the north. The burnt area covered an elevation gradient from 277 to 1419 m in a transition zone between the mesomediterranean and supramediterranean bioclimatic belts [39]. The most representative plant communities in the burnt area were (1) oak forest of Pyrenean oak Quercus pyrenaica (18.89%); (2) subserial scrub (32.76% of surface), including Arbutus unedo, Cistus ladanifer, Cytisus striatus, Cytisus oromediterraneus, Echinospartum ibericum, Erica arborea, Erica australis, Genista tridentata, and Pteridium aquilinum; and (3) pine forest of maritime pine Pinus pinaster (24.02%); and [40]. Approximately 30% of the burnt area in 2015, corresponding to the scrubland area, had been previously burnt in 2003.

Study Area and Post-Fire Management
The study area is located in the western central Iberian Peninsula (Figure 1). In August 2015, a large forest fire burned 7833 ha. The fire spread along a large valley (Rivera de Gata; hereafter Gata valley) open to the south and surrounded by steep slopes of the foothills of the Central System to the north. The burnt area covered an elevation gradient from 277 to 1419 m in a transition zone between the mesomediterranean and supramediterranean bioclimatic belts [39]. The most representative plant communities in the burnt area were (1) oak forest of Pyrenean oak Quercus pyrenaica (18.89%); (2) subserial scrub (32.76% of surface), including Arbutus unedo, Cistus ladanifer, Cytisus striatus, Cytisus oromediterraneus, Echinospartum ibericum, Erica arborea, Erica australis, Genista tridentata, and Pteridium aquilinum; and (3) pine forest of maritime pine Pinus pinaster (24.02%); and [40]. Approximately 30% of the burnt area in 2015, corresponding to the scrubland area, had been previously burnt in 2003.  The first two years after the fire, the Forest Planning and Management Service of the Regional Government (Extremadura, Spain) implemented a series of management practices in the burnt area [40]: (1) the construction of dikes, girdles, and ridges to reduce erosion; and (2) the construction and maintenance of forest roads to improve access. Moreover, several management practices were selectively applied to the three main habitat types ( Figure 1). In the oak forest, only the completely burnt trees were cut down, although these actions were suspended until 2016 in order to examine the natural regeneration and regrowth capacity of the burnt trees. The scrubland underwent no intervention and was left to natural regeneration. In the burnt pine forest area, salvage logging was carried out, including hauling trees, trimming, stripping, bucking, and stacking for later hauling ( Figure 2). The removal of burnt pine wood was performed to avoid pests and diseases, especially the spread of the pine nematode Bursaphelenchus xylophilus. All trees with more than 50% of the canopy burnt were felled, leaving only a small stand with less than 20% of the canopy burnt and located on a slope greater than 20%. The remediation in the pine forest area began immediately after the fire in September 2015 and spanned the autumn and winter months. The first two years after the fire, the Forest Planning and Management Service of the Regional Government (Extremadura, Spain) implemented a series of management practices in the burnt area [40]: (1) the construction of dikes, girdles, and ridges to reduce erosion; and (2) the construction and maintenance of forest roads to improve access. Moreover, several management practices were selectively applied to the three main habitat types ( Figure 1). In the oak forest, only the completely burnt trees were cut down, although these actions were suspended until 2016 in order to examine the natural regeneration and regrowth capacity of the burnt trees. The scrubland underwent no intervention and was left to natural regeneration. In the burnt pine forest area, salvage logging was carried out, including hauling trees, trimming, stripping, bucking, and stacking for later hauling (Figure 2). The removal of burnt pine wood was performed to avoid pests and diseases, especially the spread of the pine nematode Bursaphelenchus xylophilus. All trees with more than 50% of the canopy burnt were felled, leaving only a small stand with less than 20% of the canopy burnt and located on a slope greater than 20%. The remediation in the pine forest area began immediately after the fire in September 2015 and spanned the autumn and winter months.

Reptile Sampling Method
Reptile sampling started in spring 2016, i.e., less than one year after the fire, and continued yearly until spring 2020. Surveys were made by the line-distance sampling technique following a stratified design (see a similar procedure in [18]). In the three habitat types, transects were laid out ( Figure 3) to coincide with the post-fire management units: pine plantation (salvage logging), oak forest (thinning out of burnt trees), and scrubland (no intervention). For each of these habitat types, four linear transects were established at least 1000 m apart within the perimeter of the fire. The same number of transects were

Reptile Sampling Method
Reptile sampling started in spring 2016, i.e., less than one year after the fire, and continued yearly until spring 2020. Surveys were made by the line-distance sampling technique following a stratified design (see a similar procedure in [18]). In the three habitat types, transects were laid out ( Figure 3) to coincide with the post-fire management units: pine plantation (salvage logging), oak forest (thinning out of burnt trees), and scrubland (no intervention). For each of these habitat types, four linear transects were established at least 1000 m apart within the perimeter of the fire. The same number of transects were drawn outside the perimeter of the burnt area with equivalent characteristics: type of pre-fire vegetation, topographic orientation, and elevation range ( Figure 1). The comparison of the reptile community between burnt and unburnt transects over the five sampling years was used to assess reptile resilience within each vegetation type. Thus, a total of 24 transects (average length 1465.6 ± 336.8 m) were surveyed from the year after the wildfire in 2016 until 2020. drawn outside the perimeter of the burnt area with equivalent characteristics: type of prefire vegetation, topographic orientation, and elevation range ( Figure 1). The comparison of the reptile community between burnt and unburnt transects over the five sampling years was used to assess reptile resilience within each vegetation type. Thus, a total of 24 transects (average length 1465.6 ± 336.8 m) were surveyed from the year after the wildfire in 2016 until 2020. During these five years of sampling, surveys were made in late spring, between the months of May and June, when environmental temperatures were above 20 °C in order to ensure maximum biological activity. Each transect was sampled for 45 min each spring to ensure a reliable comparative analysis. Within the transects, searches concentrated on microhabitats where reptile individuals were most likely to be found. For each observation, the species and its perpendicular distance from the transect line at the time of sighting were recorded.
All reptiles were identified to the species level. However, the wall lizard Podarcis hispanica has emerged as a complex case with several species inhabiting the Iberian Peninsula, two of these, P. guadarramae and P. virescens, showing a contact zone between both species ranges in our study area [41]. Given their morphological similarity and the lack of data on the limits of their respective distributions, both species are included in this manuscript under the common genus Podarcis spp. When it was not possible to identify a reptile to species level, records were noted as "not identified" and computed only for abundance purposes.
To infer aboveground biomass, we used the Normalised Difference Vegetation Index (NDVI) taken from Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data. Environmental variables such as those captured by remote sensing techniques are correlated to above-ground primary productivity [42,43] or landscape change caused by fire [44]. NDVI is a proxy of photosynthetic activity and green biomass [45] and it has been used to monitor post-fire dynamics of the habitat occupied by reptile species and other animal organisms [46]. In this study, we used the imagery from days 129 to 177 (spring), so that it included the season of maximum reptile activity and also the dates of the surveys, to calculate the mean of pixel values spatially coinciding with transects for the sampling years (2016-2020) and for the spring before the fire (2015). Thus, we had one NDVI During these five years of sampling, surveys were made in late spring, between the months of May and June, when environmental temperatures were above 20 • C in order to ensure maximum biological activity. Each transect was sampled for 45 min each spring to ensure a reliable comparative analysis. Within the transects, searches concentrated on microhabitats where reptile individuals were most likely to be found. For each observation, the species and its perpendicular distance from the transect line at the time of sighting were recorded.
All reptiles were identified to the species level. However, the wall lizard Podarcis hispanica has emerged as a complex case with several species inhabiting the Iberian Peninsula, two of these, P. guadarramae and P. virescens, showing a contact zone between both species ranges in our study area [41]. Given their morphological similarity and the lack of data on the limits of their respective distributions, both species are included in this manuscript under the common genus Podarcis spp. When it was not possible to identify a reptile to species level, records were noted as "not identified" and computed only for abundance purposes.
To infer aboveground biomass, we used the Normalised Difference Vegetation Index (NDVI) taken from Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data. Environmental variables such as those captured by remote sensing techniques are correlated to above-ground primary productivity [42,43] or landscape change caused by fire [44]. NDVI is a proxy of photosynthetic activity and green biomass [45] and it has been used to monitor post-fire dynamics of the habitat occupied by reptile species and other animal organisms [46]. In this study, we used the imagery from days 129 to 177 (spring), so that it included the season of maximum reptile activity and also the dates of the surveys, to calculate the mean of pixel values spatially coinciding with transects for the sampling years (2016-2020) and for the spring before the fire (2015). Thus, we had one NDVI value per transect and year. The NDVI values were transformed as a percentage in relation to the pre-fire value of 2015 to identify the anomaly caused by the fire and the dynamics prevailing over the five-year post-fire chronosequence.

Statistical Analyses
We checked whether NDVI values varied spatiotemporally in the study area, i.e., between burnt and unburnt transects, elevation and vegetation type, and over the five sampling years. For this, we modeled NDVI values with a General Linear Model (GLM) using Year, Vegetation type, and Fire condition (i.e., burnt and unburnt) as predictor variables.
On the basis of reptile species occurrence and abundance, we assessed whether differences in the composition of the reptile community were caused by vegetation type (i.e., oak forest, pine plantation, and scrubland classes), fire condition (i.e., burnt and unburnt), and the interaction of the two factors. Thus, we performed a non-parametric permutational method of multivariate analysis of variance (PERMANOVA; [47]) to determine the similarity in reptile composition between pairs of transects. The similarity was calculated with the Bray Curtis index, which takes into account species abundance values, and the Jaccard index, which uses only the presence of species. Two temporal points were considered to perform PERMANOVAs: only the first year after the fire, i.e., 2016 data, to analyze the impact of the fire over the short term; and the entire period, i.e., 2016 to 2020, to explore whether the fire impact persists over the medium term.
General Linear Models (GLMs) were used to assess the response of the reptile community to the wildfire. For dependent variables, we selected the following reptile metrics: (1) Abundance as the total sum of reptile observations per transect; (2) Lizard abundance as the sum of all lizard observations per transect, to reduce the possibility of false negatives when considering more cryptic species; (3) Species richness as the number of species per transect); and (4) the Shannon index, representing both species diversity and abundance in each transect. The species taken into account to calculate Lizard abundance were the spiny-footed lizard Acanthodactylus erythrurus, wall lizards of the genus Podarcis sp., the Algerian sand lizard Psammodromus algirus, and the western sand lizard Psammodromus occidentalis. Although GLM results for Abundance and Lizard abundance were very close, they revealed that the latter returned the best-fit models, i.e., the lowest Akaike Information Criterion AIC and the p-values obtained for the best model (see below). Therefore, in the results section, we described only the p-values generated with Lizard abundance. After checking for overdispersion, a Poisson distribution was assumed for the response variables, except for the Shannon index with a Gaussian distribution (after calculating the inverse of the natural logarithm of the sum of its values plus 1 2 ). The predictors selected to build GLM models were NDVI scores, Elevation, Year, Vegetation type, and Fire condition. Post-fire management was not included as an independent variable in the GLMs since it was represented by the categories of vegetation type in burnt transects: pine plantation (SL), oak forest (Thinning burned out trees), and scrubland (no intervention). Likewise, in the GLM, the following interactions between pairs of variables were included: Fire × Vegetation type, Fire × NDVI evolution, Fire × Year, NDVI evolution × Vegetation type, NDVI evolution × Year, and Vegetation type × Year.
The dredge function enabled us to perform a model-averaging approach, testing all possible models and then ordering them by AIC values. Based on this information, the best models were considered those with AIC delta below 2, and the best model the one with AIC delta below 2 and higher R 2 . Furthermore, each of the fixed factors and their interactions were ordered according to the sum of the weight of all the models where each factor occurred [48]. Finally, the variance inflation factor (VIF) was examined in advance to evaluate biases due to multicollinearity. All explanatory variables registered a VIF < 2, far from exceeding the value of 10, indicating that explanatory variables have independent information and therefore do not lead to biased estimation [49]. We used the vegan package [50] to run the GLM, and MuMln package version 1.43.17 [51] to introduce the dredge function, both within the R statistical environment (version 4.0.2) [52].
In a preliminary analysis, we checked the suitability of using a reptile-detection function (see a similar approach in [36]) to improve the model fit. We used Generalized Linear Models (GLZs) to test differences in reptile detection (perpendicular distance values) between burnt and unburnt transects and among vegetation types. This model was run only for Psammodromus algirus, because it was the commonest species (76.3% of the total observations) and showed a balanced number of sightings among vegetation types and burnt vs. unburnt transects. GLZ showed no differences in the perpendicular distance for the factors considered (Fire p-value = 0.440; Vegetation type p-value = 0.177; Fire × Vegetation type p-value = 0.340), and therefore the use of distance functions was not justified to fit the models according to the species detectability. All models in this preliminary analysis were run by the SPSS Statistics software [53].

NDVI Spatiotemporal Variation
A GLM (overall R 2 = 0.842; F = 20.772; p < 0.001) showed that NDVI values were affected by all fixed factors and interactions of the model: Fire condition (p < 0.001), Vegetation type (p < 0.001), Year (p < 0.001), Elevation (p < 0.001), Fire × Vegetation type (p < 0.001), Fire × Year (p < 0.001), and Vegetation type × Year (p = 0.014). Based on these results, NDVI reflects an environmental impact by fire which is nuanced by the rest of the variables, so we included it in the reptile models as an explanatory variable.
NDVI was strongly determined by the interaction of fire occurrence and vegetation type, with a smoothing effect in the oak forest zones (Supplementary Material Figure S1). The third year after the fire, NDVI values did not differ between burnt and unburnt transects in oak and scrubland areas. By contrast, burnt transects in pine plantations had lower NDVI values than did unburnt transects even five years after the fire (Supplementary Material Figure S1).

Reptile Community Composition
In transects laid out between 2016 and 2020, 728 reptiles were observed across 12 species (Supplementary Materials Table S1). Among these, four lizard species (Acanthodactylus erythrurus, Podarcis sp., Psammodromus algirus and Psammodromus occidentalis) represented 90.2% of the records. Natrix astreptophora and Vipera latastei were detected only in unburnt transects, while Psammodromus occidentalis was detected only in burnt ones (Supplementary Materials Table S1). Overall, more species were detected in transects laid out in the burnt area for all vegetation types (Supplementary Materials Table S1).
The lizard species Psammodromus algirus and Timon lepidus were the only species detected in all years. For burnt scrubland and oak forest transects, most species were detected in the first years after the fire, while in the burnt pine transects most species appeared only after the third year post-fire (Supplementary Materials Tables S1 and S2).
For the first year after the fire, the PERMANOVA results revealed differences for the abundance of each species (Bray Curtis index) due to the type of vegetation, fire condition, and the interaction between the two factors ( Table 1). The results for species composition (Jaccard index) were equivalent for all factors analyzed. For the period 2016-2020, the PERMANOVA showed differences due only to the type of vegetation, both in abundance and in species composition (Table 1). Table 1. Permutational analysis of variance (PERMANOVA) of the reptile community between unburnt and burnt plots at Gata valley. The reptile community was compared in terms of species composition (using species presence/absence with the Jaccard index) and species abundance (using the relative abundance of species with the Bray Curtis index). The PERMANOVA was performed one year after the fire (2016) and all the chronosequence (2016-2020). The independent variables were vegetation type (oak forest, pine plantation, and scrubland), fire condition (burnt and unburnt), and their interaction.

Lizard-Abundance Models
Four explanatory factors were important in explaining lizard abundance per transect (in descending order): vegetation type, NDVI, fire condition, and the interaction NDVI and vegetation type (Figure 4). The best-fitting models (Supplementary Materials Table S3) explained a high proportion of the variance of lizard abundance (R 2 range from 0.625 to 0.610), and the best model confirmed the effect of these factors plus the interaction NDVI × Fire into lizard abundance ( Table 2). Regarding vegetation type, oak forest was the habitat with the highest lizard abundance (Figure 5a Figure 5b). Lizard abundance increased in transects with higher NDVI values, and this relationship was strong in oak transects and weaker in pine and scrubland transects (Figure 5c).

Reptile-Diversity Models
Vegetation type best-captured variability in the Shannon index, followed by fire condition (Figure 6a). The best-fitting models explained a low proportion of the variance of the Shannon index (maximum R 2 = 0.120; Supplementary Materials Table S4), and the best model showed differences on the Shannon index only for Vegetation type and the interaction between Fire x Vegetation type (Table 3). Thus, the Shannon index was higher in  Table 2. Results of the best model for lizard abundance in the Gata fire based on model averaging procedure and data dredging. Significant codes: "***" <0.001; "**" <0.01; "*" <0.05; "." <0.1.

Reptile-Diversity Models
Vegetation type best-captured variability in the Shannon index, followed by fire condition (Figure 6a). The best-fitting models explained a low proportion of the variance of the Shannon index (maximum R 2 = 0.120; Supplementary Materials Table S4), and the best model showed differences on the Shannon index only for Vegetation type and the interaction between Fire x Vegetation type (Table 3). Thus, the Shannon index was higher in scrubland (average 0.52; range 0-1.39), followed by oak forest (average 0.35; range 0-1.03) and finally pine plantation (average 0.19; range 0-1.01). From burnt to unburnt areas, this index decreased in scrubland and pine plantation, while in oak forest it showed a slight increase (Figure 6b).

Discussion
Our study demonstrates that the habitat context (i.e., vegetation type) is the most important factor explaining differences in the composition of the reptile community. Fire operated only as a major disturbance, prompting variation in the reptile community for the first year after the fire. This result shows that the reptile community in fire-prone areas  Table 3. Results of the best model for Shannon index in the Gata fire based on model averaging procedure and data dredging. Significant codes: "***" <0.001; "**" <0.01; "*" <0.05; "." <0.1. Elevation and vegetation type were the most important explanatory factors in speciesrichness models (Supplementary Materials Table S4), even though no variable had a high cumulative weight (all below 0.6), and the best fitting models explained a low proportion of variance (maximum R 2 = 0.059; Supplementary Materials Table S5).

Discussion
Our study demonstrates that the habitat context (i.e., vegetation type) is the most important factor explaining differences in the composition of the reptile community. Fire operated only as a major disturbance, prompting variation in the reptile community for the first year after the fire. This result shows that the reptile community in fire-prone areas such as the Mediterranean basin is resilient to fire (see similar results in [15,54]). This resilience however depends on the dominant vegetation type: in natural habitats (e.g., oak forest), NDVI and lizard abundance recover swiftly after the fire. By contrast, pine plantations tend to reach their pre-fire state slowly (according to the comparison to unburnt transects), and this gradual recovery affects both NDVI and reptile occurrence.

Reptile Community Composition in Relation to Vegetation Type
Our results imply that the composition of the reptile community in this study area is determined primarily by vegetation type and not by an ordered succession linked to fire occurrence [14,[35][36][37]55]. That is, the reptile community is closely linked to the environmental attributes of the habitat. For example, vegetation structure (such as canopy cover or vegetation density) is directly related to the potential for thermoregulation, and therefore it is a central feature for understanding the capacity of the environment to accommodate reptiles [56]. In these habitats, reptile species tend to select adequate microhabitats [57] in line with their thermal preference, water balance, and physiological attributes [58]. Habitat characteristics also determine the availability of food resources, and this constitutes a prime factor explaining the occurrence of predator species such as snakes [59]. In conclusion, environmental attributes of different vegetation types in the Gata valley shape reptile assemblages and bolster reptile beta diversity at the landscape scale.
In our study, natural habitats (oak forest and scrubland) appeared to offer better conditions in both parameters of reptile abundance and diversity. In contrast, pine plantation was characterized by sparse understory vegetation, a high canopy, and low food (arthropod) availability to lizards [20,60]. Moreover, the forestry management in monospecific pine plantations creates a homogeneous landscape with a lack of refuges for fauna. All these features discourage the appearance of reptile species and define the pine plantation as a low-quality habitat for these and many other organisms [61][62][63].

Reptile-Community Resilience to Fire
Habitat shifts after a wildfire result in greater impacts on reptiles than does direct mortality by fire [64]. Thus, we point to the high regeneration ability of the Mediterranean vegetation as a key factor supporting the recovery of animal communities throughout the food web (e.g., the quick regrowth of ground vegetation facilitates recolonization by arthropods, which in turn, as prey, promote the resilience of the reptile community; [54]). Most Mediterranean ecosystems are well adapted to fire, and plant communities can recover by seeder or resprouter strategies [7,8,65]. The reptile community in the Gata valley responded negatively one year after the fire in accordance with the sweeping habitat transformation but recovered later over the five-year surveys. This resilience is, however, habitat dependent as the resilience changed among the three vegetation types in the study area: the scrubland was composed of several resprout and facultative plants (i.e., species with thermos-dehiscent seeds, plus the capacity of replacing the aerial part) while the oak forest was dominated by Quercus pyrenaica, which has an ample resprout capacity as well as a high tree survival rate [66]. In both cases, the pre-fire plant community was quickly re-established. The post-fire regeneration of Pinus pinaster, a fire-sensitive habitat in the Mediterranean basin [7], relies on the serotinous cones and seed recruitment. The reduction in lizard abundance (over the short term) and differences in reptile assemblage (over the medium term) were far more marked in pine plantation than in the oak forest and scrubland. In short, differences in the resilience of the three vegetation types in turn influence the resilience of the reptile community living in these habitats.

SL Implications of Salvage Logging in the Pine Plantation
The NDVI has been used to monitor post-fire habitat dynamics [46,67], and it is considered as a reliable predictor of habitat recovery after fire for many ecosystem components [68]. Our results support this conclusion as we detected a general reduction of the NDVI index in all burnt areas one year after the fire. Moreover, NDVI reduced drastically in the pine plantation the first year post-fire, and its recovery was far below the levels of scrubland and oak forest. The NDVI provides a proxy of photosynthetic activity and green biomass [45]. The increase in primary plant production and, associated with it, greater food availability for lizards (i.e., most Mediterranean lizard species are insectivorous) can be the main cause to explain the positive relationship between NDVI values and lizard abundance. SL was applied in the whole burnt pine plantation. This vegetation type exhibited a slow NDVI recovery compared to the unburnt control plantation. In contrast, NDVI recovered fast on the other two vegetation types. We acknowledge that we have not sampled reptiles in unlogged burnt pine plantation due to the absence of unlogged burnt patches. However, the comparison with the rest of the vegetation types affected by the same disturbance and environmental conditions suggests that SL (i.e., the complete removal of burnt wood and additional mechanical actions on the soil) caused the slow NDVI recovery in pine plantation burnt patches, and indirectly, the slow reptile recovery.
No universal consensus is available concerning how SL affects the regeneration of woody plants after the fire, as its possible impacts depend on local characteristics and species [69]. For example, SL may benefit open-habitat species [70,71]. In the Gata valley, some logged lowlands were quickly occupied by scrublands (pers. observation), although SL globally delayed the woody regeneration, and most reptile species were found in burnt pine transects only three years after the fire. The scientific literature reports multiple direct and indirect effects on habitat characteristics that appear to support our results [26]. For instance, the drastic reduction of refuges caused by the use of machinery has a negative impact on the survival of fire-sensitive species within large burnt areas [72]. In addition, the absence of biological legacies (such as wood debris, burnt trunks, or survival vegetation) that provide services to the ecosystem recovery (such as soil retention, seed traps, energy, and nutrients; [23,26]) may be detrimental to the main taxonomic groups (i.e., arthropods and gastropods) preyed on by reptiles [73]. Reptiles are low-dispersal organisms, and the impact of SL on many biotic and abiotic ecosystem components can expose reptile species to local extinction and compromise their conservation whether or not the frequency and extent of fire increase.

Species-Specific Responses
Reptiles display contrasting species-specific responses to fire, with some species disappearing or becoming less frequent [35,38,74] and others benefiting from post-fire habitat changes [13,75]. For example, two species of snakes did not occur in burnt transects: Natrix astreptophora occurs in areas of low thermal exposure and high moisture [76] and Vipera latastei is an ambush predator that inhabits undisturbed areas [77]. Despite their anecdotal occurrence in the transects, they are specialist species and consequently can be more affected by fire than are generalist snakes [78]. On the contrary, Psammodromus occidentalis, a lizard detected only in burnt transects, takes advantage of degraded habitats consisting of open scrubland instead of the original tree cover [79]. One of the potential benefits for ectothermic animals such as reptiles is a higher incidence of solar radiation at levels closer to the ground due to the reduction of the tree canopy [80], which benefits heliothermic species such as the lizard Timon lepidus, which has been detected only in the burnt forest transects.

Conclusions
Vegetation type was the main factor explaining differences in the reptile community composition. The structure of each habitat studied created microhabitats with particular environmental features which are adequate for reptile species according to their physiological characteristics. In contrast, fire impacted the reptile community only in the short term, and this result indicated that the reptile community in the study area is resilient to fire. We have also verified that natural forests had more resilient communities of reptiles to fire compared to monospecific plantations for commercial purposes of Pinus sp. Moreover, salvage logging delayed vegetation recovery and, consequently, the recovery of associated fauna (i.e., reptile community). This result is important as salvage logging is a generalized practice in many burnt pine plantations in the Mediterranean biome. Both pre-fire habitat legacy and post-fire management are key factors to understanding reptile resilience to fire. Thus, post-fire management (e.g., the amount of burnt wood removal) is a fundamental issue for conservation purposes under the current fire regimes shift in the Mediterranean region.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/f12111487/s1, Figure S1. Temporal variation of NDVI values in burnt (red) and unburnt (green) transects along the chronosequence (2016-2020) plus the reference year prior to the fire (2015), at each category of vegetation type. For each year, the average value ± standard deviation (whiskers) is given; Figure S2. Photographs of reptile species found only in unburnt transects: Natrix astreptophora [a] and Vipera latastei [b]; and not found in unburnt forest transects: Psammadromus occidentalis [c] and Timon lepidus [d]; Table S1. Total amount of reptile species found by vegetation type in burnt and unburnt tran-sects in Gata valley from 2016 to 2020. In bold, species only detected either in burnt (BU) or un-burnt (UN) transects within each vegetation type; Table S2. Temporal variation in the occurrence of reptile species according to the three vegeta-tion types in Gata valley. The presence of a species in each type of transect is marked with an X and those detected for the first time are highlighted in bold. The columns on the left indicate the burnt transects, and on the right the unburnt transects. The last row indicates the total number of species found. The acronyms of the species are indicated in Table S1; Table S3. Model ranking of the variables that best explain lizard abundance per transect and year in Gata valley. Lizard sampling was conducted in burnt and unburnt plots for five years af-ter the fire. The ranking was based on AICc values and only models with ∆ AICc < 2 are consid-ered. Each row represents a model, and includes a description of the parameters and the varia-bles used. Best model (∆ AIC < 2 and higher R2) is shaded; Table S4. Model ranking of the variables that best explain reptile diversity (Shannon index) per transect and year at Gata valley. Reptile sampling was conducted in burnt and unburnt plots for five years after the fire. The ranking was based on AICc values and only models with ∆ AICc < 2 are considered. Each row represents a model, and includes a description of the parameters and the variables used. Best model (∆ AIC < 2 and higher R2) is shaded; Table S5. Model ranking of the variables that best explain reptile species richness per transect and year in Gata valley. Reptile sampling was conducted in burnt and unburnt plots for five years after the fire. The ranking was based on AICc values and only models with ∆ AICc < 2 are considered. Each row represents a model, and includes a description of the parameters and the variables used. Best model (∆ AIC < 2 and higher R2) is shaded.