Ecological niche models applied to post-megare vegetation restoration in the context of climate change

Climate change and land-use changes are the main drivers altering re regime and leading to the occurrence of megares. Current management policies focus on short-term restoration without considering how climate change affects regeneration dynamics. We aimed to test the usefulness of ecological niche models (ENMs) in post-re restoration plans, assessing the short and long-term effects (climate change) of tree species distribution. We, also, examined different important conceptual and methodological aspects. We executed ENM for the four main tree species in an area affected by a megare in Central Spain following a hierarchical approach and precise resolution (25 meters), at two scales (local and regional), two sampling strategies (regular 2 years after the re, and opportunistic 14 years after the re) at local scale, and under two future climate change scenarios. m resolution: aspect, northerliness (orientation along the north-south axis), monthly solar radiation, distance from streams, compound topographic index, heat load index, linear aspect, elevation, curvature and slope. We used the same approach as for the regional predictors based on Pearson's pairwise correlations to prevent multicollinearity, nally retaining the six most uncorrelated and ecologically meaningful variables for the plant species in the study area: (1) Distance to streams. For each pixel, we calculated the shortest surface distance (in m) to any stream (Fleming and Doan 2013); (2) Heat load index. We used the method proposed by McCune and Keon (2002), which assumes warmer temperatures on southwest-facing slopes than on southeast-facing slopes and takes into account the steepness of the slope; (3) Elevation, derived from a 25 m resolution digital elevation model (DEM, Spanish National Geographic Information Centre); (4) Slope, derived from the DEM; (5) Solar radiation in August, derived from the DEM using ArcGIS 10.3; and (6) Curvature, derived from the DEM using ArcGIS 10.3. All independent variables were projected at ETRS89


Introduction
Wild res play a major role in landscape dynamics in Mediterranean ecosystems (Bond et al. 2005;Keeley 2012). However, re regimes are affected by current climate change (Jolly et al. 2015;Dupuy et al. 2020) and human activities (Knorr et al. 2014;Aquilué et al. 2019). The extreme drought events caused by climate change lead to alterations in the behaviour and physiological conditions of plant species, as well as in their distribution and recruitment, thus affecting forest resilience (Lindner et al. 2010;Stevens-Rumann et al. 2018;Resco de Dios 2020). In addition, abandonment of the rural environment has led to an increase in forest mass and fuel connectivity (Ortega et al. 2012;Madrigal et al. 2017). All of these factors exacerbate the wild re regime (Ortega et al. 2012), increasing the number and frequency (Pellegrini et al. 2021) of extreme events such as mega res (Stephens et al. 2014;Tedim et al. 2018) and resulting in loss of human lives, along with severe economic and environmental damage (San-Miguel- Ayanz et al. 2013) that threatens the conservation of different ecological communities in the Mediterranean region (Moriondo et al. 2006;Nogués-Bravo and Rahbek 2011;Alexander et al. 2016).
Decision-making regarding the restoration of vegetation in areas affected by mega res is currently carried out according the known autoecological characteristics of the species selected for restoration and to static criteria based on maps of potential vegetation and forest maps (Pemán García et al. 2006;Vallejo et al. 2012). However, the ways in which plant community dynamics will adapt to climate change in the future are not taken into account (Thom et al. 2017). Post-re management strategies mainly focus on short-term recovery by curbing the effects of erosion and re-establishing the community existing before the re (Lucas-Borja et al. 2021). This approach can give rise to great controversy because ecological and management preferences often oppose each other (e.g. salvage logging, DeRose and . In addition, long-term management is complicated by the uncertainty of potential changes in the species distribution and ecological conditions at different scales (microsite stand and landscape) during regeneration processes . Hence, in post-re management it is important to consider how climate change affects the potential distribution of tree species, with the aim of enhancing the natural evolution of these species (Vallejo and Alloza 2019).
Ecological modelling is a powerful tool in numerous elds , and ecological niche models (ENMs) have increasingly been used in recent years. One of the main applications of these models is to study how climate change affects vegetation (e.g. Thuiller et al. 2008;Engler et al. 2011;Felicísimo et al. 2011). Although ENMs are potentially useful in ecological restoration, they have only occasionally been applied (Gastón et al. 2014;Mateo et al. 2019a) and have not included future climate scenarios. Many conceptual and methodological issues therefore remain unresolved in this context. In this study, we raise the main hypothesis that ENMs are useful and applicable tools in post-re dynamic restoration programmes, considering the consequences of climate change on the spatial distribution of the plant species. The resulting species suitability maps (i.e. potential distribution), in combination with other sources of information (i.e. remote sensing, eld work and drone data, etc.), increase the amount of information available to managers and could lead to a signi cant improvement in post-re restoration plans in the current context of climate change (San-Miguel- Ayanz et al. 2012). To corroborate this hypothesis, we addressed three critical conceptual and methodological issues. Page 4/23 First, processes that are important for plant post-re regeneration occur in two sequential phases: (1) the immediate post-disturbance phase; and (2) later succession, when species interactions and resource competition determine abundance (Noble and Slatyer 1980). However, ENMs assume that there is an equilibrium (or at least a pseudo-equilibrium) between the vegetation and the environment dynamics (Guisan and Zimmermann 2000;Hampe 2004) but are not able to represent the progress of these sequential phases. In the present study, we investigated this aspect by collecting data at two different times: 1) shortly after the disturbance (2 years after the mega re) and 2) at a later stage of vegetation regeneration (14 years after the mega re). Our hypothesis is that data collected at an early stage of regeneration, when the vegetation community dynamics are not stable, will only provide a partial view of the potential distribution of plant species. Therefore, the ENMs generated will not accurately re ect the potential distribution of the species considered (Guisan and Zimmermann 2000), and the information derived will not be useful for post-re restoration plans.
Second, spatial scale plays an important role in identifying the ecological niche of species during the calibration of the ENMs (Thuiller et al. 2004). Thus, ENMs tted across small geographic areas will not include the complete range of climatic conditions suitable for the species, and the climate change projections thus generated will be inaccurate (Petitpierre et al. 2016). A larger scale, such as a regional or continental scale, will cover the complete climatic niche of the species (Pearson et al. 2004;Petitpierre et al. 2016). Nevertheless, in restoration programmes, the accuracy of predictions must be able to support decision-making at stand level by including information on the composition of the species and considering local factors (e.g. microclimate, lithology) affecting species distribution. Hence, forest management requires relatively ne grain resolution at local scale (Randin et al. 2009;Gastón and García-Viñas 2010;Gastón et al. 2014). To address this issue, we developed ENMs at two scales (regional and local). We initially assumed that the models obtained at these scales will be different. Nevertheless, ENMs developed at both scales can produce important information for post-re restoration. Regional models are expected to provide information about species distribution according to the whole climatic niche and more accurate future climate projections (Petitpierre et al. 2016). However, local models are expected to represent how species are distributed under the current climate scenario, and they will capture important local conditions for vegetation (Mateo et al. 2019b).
Third, eld data are usually obtained through regular systematic sampling strategies (Hirzel and Guisan 2002) for which a large amount of resources are necessary. However, budgets in conservation and restoration programmes are always limited, making it di cult to apply this type of data in combination with modelling tools. In a previous study, an opportunistic sampling strategy was proposed in which existing roads and tracks are used for data collection, thus enabling development of ENMs with fewer resources ). In the present study, we assumed that reliable ENMs can be produced by collecting data from roads and tracks in areas affected by a mega re.
The main aim of the present study was, therefore to develop and validate high-resolution ENMs for use in a post-re restoration programme in the current context of global change and considering different future climate scenarios. With this aim, the following steps were carried out: (i) evaluation of short and long term potential distribution of plant species after a mega re by the use of ENMs and data collected at different times (2 and 14 years after the mega re); (ii) assessment of the importance of the scale in the calibration of ENMs by using two different scales (regional and local); (iii) testing the reliability of the ENMs derived from an opportunistic sampling strategy involving the use of roads and tracks.

Study area
The study was conducted in a site burned by a mega re that began in "Riba de Saelices", located in Central Spain (40º 46'-41º06 'N 1º 52'-2º 39'W), in the autonomous region of Castilla-La Mancha. The re occurred on July 2005 and burned an area of 12,874 ha. The elevation of the area affected by the re varied between 1,200 and 1,370 m. The mean annual temperature in the area is 10.5 ºC and the main annual precipitation is 468 mm. The substrate is characterized by red Buntsandstein sandstones (rich in iron), secondary limestones and siliceous sediments. Before the re, the forest stands were generally composed by conifers, dominated by Pinus pinaster Ait. together with a variety of Fagaceae fam. (Quercus pyrenaica Willd., Q. faginea L., Q. ilex L. subsp. ballota (Desf.) Samp.). These oak species all occurred in the understory and forming more or less extensive shrublands. The stands are considered to be of natural origin (Rodenales de Molina natural provenance, according to Alía and Moro 1996). Pinus pinaster stands were intensively managed from the eaerly 1900s until the beginning of the 1980s for resin production, increasing the area of natural P. pinaster forest in a landscape probably dominated by oaks in the past (Osorio Vélez et al. 2003).

Species distribution data
We selected the four most representative tree species in the study area: maritime pine (Pinus pinaster), Pyrenean oak (Quercus pyrenaica), Portuguese oak (Q. faginea subsp. faginea) and holm oak (Q. ilex subsp. ballota). We chose to study tree species because of their importance in restoration plans (Mucina et al. 2017). Based on our objectives and hypotheses, the ENMs were calibrated with data obtained at two different spatial scales: 1. Regional scale: occurrence data for the four tree species were downloaded at the European level from the Global Biodiversity Information Facility (GBIF, https://www.gbif.org) database. The administrative information was corroborated with information from administrative maps, and mismatching occurrences were removed. In addition, points with coordinates of simple zeros or a precision greater than 1000 m were excluded. Finally, the maps were visually revised and compared with bibliography. The total numbers of presences used for the calibrated models were 4,048 (P. pinaster), 2,849 (Q. faginea), 2,171 (Q. ilex) and 2,400 (Q. pyrenaica). 2. Local scale: species presence/absence data were obtained in two eld studies by using different sampling strategies and time periods. In 2007 (two years after the re), the regional government of Castilla-La Mancha conducted a vegetation inventory following a systematic sampling method, based on a grid layout over the forest area (equivalent to 71% of the study area), with nodes separated by 200 m (a total of 2,458 points were sampled). In 2019 (14 years after the re), additional opportunistic sampling (following Mateo et al. 2018) was carried out, and a total of 224 circular plots of radius 10 m were established at intervals of 800 m along the tracks and roads, with the starting point established 100 m from the re perimeter (Fig. S1). The distance of the interval (800 m) was selected considering the amount of eldwork (minimum number of plots) required to capture the environmental variability in the study area, determined by a straightforward analysis considering different distances and a principal component analysis of the environmental local variables (see Mateo et al. 2018). All of the sampling points were georeferenced using the opensource mobile application Epicollect5 (https:// ve.epicollect.net).

Environmental and climatic variables
Following our hypotheses, we used different data sets of environment variables and spatial resolutions at the two scales considered, with the premise that the spatial distributions of the plant species are expected to be determined by different drivers at different scales (McGill 2010). 1. Regional scale: regional models were calibrated with large scale variables (e.g. macroclimate) with the expectation that the whole climatic niche of the species would be captured and projections to future scenarios would be improved (Petitpierre et al. 2016;Mateo et al. 2019b (2019a), the climatic variables were downscaled from 1 km to 25 m resolution in the study area. For temperature variables, this was carried out using local linear regressions. First, the relationship between temperature and elevation was tted at 1 km resolution by using the localReg function. The intercepts and slopes were then disaggregated at 25 m and spatially smoothed with a conic density kernel. The tted relationship was then applied at 25 m across the study area. Finally, the precipitation variables were downscaled using a simple bilinear approach. All the analyses were carried out in R. 2. Local scale: local models were calibrated with environmental variables (e.g. microtopography and microclimate) with the expectation of re ecting the microsite conditions available for the plant species (Lortie et al. 2004, Mateo et al. 2017). We generated different environmental variables at 25 m resolution: aspect, northerliness (orientation along the north-south axis), monthly solar radiation, distance from streams, compound topographic index, heat load index, linear aspect, elevation, curvature and slope. We used the same approach as for the regional predictors based on Pearson's pairwise correlations to prevent multicollinearity, nally retaining the six most uncorrelated and ecologically meaningful variables for the plant species in the study area: (1) Distance to streams. For each pixel, we calculated the shortest surface distance (in m) to any stream (Fleming and Doan 2013); (2) Heat load index. We used the method proposed by McCune and Keon (2002), which assumes warmer temperatures on southwest-facing slopes than on southeast-facing slopes and takes into account the steepness of the slope; (3) Elevation, derived from a 25 m resolution digital elevation model (DEM, Spanish National Geographic Information Centre); (4) Slope, derived from the DEM; (5) Solar radiation in August, derived from the DEM using ArcGIS 10.3; and (6) Curvature, derived from the DEM using ArcGIS 10.3. All independent variables were projected at ETRS89 UTM 30N.

Ecological niche modelling
For each species, we generated an ENM for the three sources of data available at both scales considered: 1. Regional scale. The regional models (Regional) were calibrated at the European scale by using data from the GBIF database and climatic variables at 1 km resolution. They were projected to the local study area by using downscaled climate variables at 25 m resolution (see Mateo et al. 2019b).
2. Local scale. Two local models were calibrated and projected in the local study area by using environmental variables at a resolution of 25 m: (i) with the data obtained from the regular sampling conducted two years after re (Local1), and (ii) based on the opportunistic sampling carried out 14 years after the re (Local2).
For each database (Regional, Local1, Local2) and tree species, an ensemble model (Araújo and New 2007) was generated by combining three statistical techniques: generalized linear models (McCullagh and Nelder 1989), boosted regression trees (Friedman 2001) and random forests (Breiman 2001). All the models were generated using the biomod2 ) R package with default parameters. For each modelling technique, the corresponding data set was randomly split ten times into an 80% data set, to generate the models, and a 20% data set, to estimate the predictive accuracy of each (i.e. split-sample cross-validation), thus producing a total of 30 models for each species and scale. The AUC statistic was calculated for all of the models. After elimination of all models with an AUC < 0.8, we generated an ensemble model consisting of the weighted mean of the different model predictions, where the weight of each model was proportional to its predictive accuracy. We compared the consensus models generated for each species among approaches, applying the Pearson correlation coe cient. Finally, the three ensemble models generated for each species were then projected onto the two climate future scenarios.
To assess our initial assumption about the differences in the extent of the suitable distribution of the species depending on the scale, we converted original ensemble models (continuous suitability index) into binary models (presence/absence) by applying two different threshold criteria: a threshold derived from the AUC statistics and the threshold maximizing the Kappa statistics (Scherrer et al. 2018). Finally, we calculated the potential area predicted (number of pixels in the binary model) and the percentage relative to the total study area for each species, scale and time period.

Results
The nal ensemble models across all species and scales were statistically reliable, yielding AUC values > 0.9 (Table S1).

Regional and local scale models for the current scenario
Comparison of all ENMs showed large differences in the models generated across the two spatial scales (local and regional) regarding the four tree species ( Table 1). The correlations between scales were lower than 0.36 for all species except Q. pyrenaica (Local1 vs Local2 and Local2 vs Regional), which yielded a value of 0.57 (Table 1). Regional binary models predicted higher percentages of suitable area than local binary models (Local1 and Local2), reaching 100% for Q. pyrenaica, Q. faginea and P. pinaster and 57.73% for Q. ilex (Table 2). In general, the Local2 models showed greater similarity to regional models than Local1 (Tables 2, S2). Thus, the distribution maps of suitable habitat showed different spatial patterns for local and regional models (Figs. 1, S2).

Fieldwork sampling strategies
Evaluating the performance of the different strategies used to calibrate the local ENMs, i.e. Local1 (regular-systematic sampling strategy) and Local2 (opportunistic strategy), showed that both produced reliable predictions (Tables 2, S2) but with large differences between each. Local2 models (14 years after mega re) predicted larger areas than those predicted by the Local1 models (2 years after mega re): twice as large for Q. pyrenaica and P. pinaster, three times as large for Q. faginea and more than six times as large for Q. ilex (Table 2, S2). Thus, e.g. the percentage of suitable area predicted by Local1 for Q. pyrenaica was less than 37%, while the suitable area predicted by Local2 for this species was around 70% (Table 2). Table 1 Pearson's correlations for the pairwise comparison of the ecological niche models (ENMs) generated (Local1, Local2 and Regional) under the current scenario for the species considered in this study Table 2 Percentage of suitable area predicted in the study area by binary ENMs projected under the current scenario and at the different future scenarios (RCP 4.5 and RCP 8.5) according to the area under the curve (AUC) statistic for each species Projections of future scenarios The regional scale ENMs projections for two future scenarios (RCP 4.5 and RCP 8.5) showed a general decrease in the suitable area predicted (Tables 2, S2, Figs. 2, S3). The Regional binary models projected to the RCP 4.5 scenario showed a drastic reduction in the suitable area predicted for holm oak (percentages = 0), a one-third reduction for maritime pine (percentages = 31.22) and a reduction of half for Portuguese oak (percentages = 51.37; Table 2). Maintenance of close to 100% suitable areas was only predicted for Pyrenean oak (Tables 2, S2). The ENMs projected to the RCP 8.5 scenario showed a signi cant reduction in the suitable area predicted across the four species in the regional models (percentages = 0; Tables 2, S2), representing future climatic conditions that would not allow optimal development of the species. On the other hand, local models (Local1 and Local2) predicted a general increase in suitable area, except for Pyrenean oak at RCP 8.5 in both Local1 and Local2 and at RCP 4.5 in Local2 (Tables 2, S2). Figure 1 (a) Regional ENM for Pinus pinaster (b) Local1 ENM for P. pinaster (short-term post-re regeneration) (c) Local2 for P. pinaster ENM (long-term post-re regeneration). Suitable distribution index under the current scenario. All of the models were rescaled from 0 to 100 potential suitability values Figure 2 Suitable distribution index for ENMs developed for maritime pine (P. pinaster) under the different future scenarios (RCP 4.5 and RCP 8.5). All of the models were rescaled from 0 to 100 potential suitability values. Regional, Local1 (short-term post-re regeneration) and Local2 (long-term post-re regeneration) are shown. The regional model under the RCP 8.5 scenario did not predict the presence of the species

Discussion
In this study, we have shown the potential application of ENMs as a support tool in decision-making in a post-mega re restoration programme, by representing the potential distribution of the most representative tree species in an area burned by a mega re.
Short and long term effects of re on tree species distribution Regeneration of the plant community is guided by two sequential phases (Kayes et al. 2010), and the design of forest management plans after mega res should contemplate both phases. In addition, the information should be collected during the most appropriate period for designing the most accurate ecological models. Here, we compared two eld trials that differ temporarily by 12 years (Local1 and Local2) at the same local scale and which showed differences in the spatial patterns (Figs. 1, S2), as we hyphothesized. The cause of the variations in these distributions can be attributed to the stabilization of the plant dynamics based on the regeneration traits of the different resprouting and seeding strategies (San-Miguel-Ayanz et al. 2012;Vallejo et al. 2012;Pausas and Keeley 2014).
Among the species studied, only maritime pine is an obligate seeder. During the rst phase of seed recruitment, immediately after re, the distribution of this species is related to the available seed bank and re severity (Vega et al. 2008, Madrigal et al. 2010, with a remarkable regeneration phase that depends on new individuals with lower growth than resprouters (Pausas and Keeley 2014). In fact, the resprouter species (Quercus spp.) showed a rst consistent regeneration phase with signi cantly greater growth than maritime pine (Madrigal et al. 2011), which explains the high correlation observed for Pyrenean oak on comparing both Local1 and Local2 models ( Table 1). Survival of resprouters (Quercus spp.) (at least the root) depends on the individuals present before the re, while consolidation of the seeded progeny (maritime pine) depends on other factors such as post-re meteorology, severity, post-re mortality and ecophysiological factors (Calvo et al. 2003;Vega et al. 2008;Madrigal et al. 2010). Thus, in the representation of the distribution of the suitable areas for maritime pine (Fig. 1) and Quercus spp. (Fig. S2), the initial dominance of the resprouters is not re ected by the models due to its dispersed distribution (presence of shoots and roots) relative to seeders (seed rain after the re).
The study ndings con rm that eldwork at an advanced phase of vegetation succession provided more useful information for decision-making in post-re restoration plans than when eldwork was only carried out shortly after the re, due to stabilization of the regeneration phase.

Importance of the scale and effect of different climate change scenarios
The results were consistent with our initial hyphothesis regarding the effect of scale under different climate change scenarios. The dissimilarities in future projections according to the different scales (Fig. 2, S3) were generated by different climatic and environmental drivers (Mateo et al. 2019a) and by differences in the extent used for their calibration. Thus, the ENMs calibrated with local scale occurrence data and landscape enviromental variables (Local1 and Local2; Figs. 1, S2) re ected the local ecological conditions for the species (Randin et al. 2009), but failed to represent the whole ecological niche because they were calibrated with a partial climatic niche of the species, generating biased future projections (Thuiller et al. 2004;Petitpierre et al. 2016). This nding implies that ENMs calibrated at local scale are limited by partial representation of the species ecological niche, and therefore future projections may be biased.
The ENMs at a regional scale covered a greater extension of the species climatic niche (Pearson et al. 2004) and were calibrated with macroclimatic variables. Therefore, the effect of climate change on the species can be projected to different future climate scenarios. However, the regional models provided a poorer representation of the territory affected by the mega re, over-predicting species occurrences (Petitpierre et al. 2016;Mateo et al. 2019a). By way of illustration of this, the distribution maps for maritime pine models (Fig. 1) show different spatial patterns, where the highest suitable presences are more widespread in the regional models than in the local models, while in the local models the most suitable areas are mainly concentrated in the southeast of the burned area. The wide range of information provided by the ENMs calibrated across the different scales can improve decision-making in post-re management. Models calibrated at a regional scale represent the whole species climatic niche, providing information to managers about the presence of species based on climatic suitability. These data can help to characterize burned areas on the basis of predictions in future scenarios. Thus, we consider that the regional models yield more accurate ts than local models due to greater ecological niche representation. On the other hand, models calibrated at local scales provide information about the current distribution of the community after the mega re, which is useful for managers in regard to the most suitable restoration measures.
The regional ENM scale projections for all species projected in the future scenario with higher concentrations of greenhouse gases (RCP 8.5) were not represented in suitable distribution maps due to the low percentage of suitable area predicted according to climatic conditions in this scenario (Table 2), which is not favourable for the development of this species.
However, despite the strong capacity of the ENMs projected to future climate change scenarios, some limitations can be recognised. Future model predictions are based on the current distributions of the species and the micro/macroclimatic conditions and also environmental factors, and they are limited by not considering other types of variables that are critical in the dynamics of plant community distributions, e.g. soil type, land-cover type and biotic interactions (Pearson and Dawson 2003;Mateo et al. 2011). In this study, the suitable area distribution maps generated by ENMs projected to future climate change scenarios are proposed as a support tool that can be used in combination with other sources of information and techniques (remote sensing, soil studies by eldwork, etc.) to expand the availability of information. The ENMs re ect how species will behave according to climate change in the near future in order to design resistant and resilient landscapes, making post-re restoration plans more dynamic, which is not currently contemplated (Jacobs et al. 2015;Löf et al. 2019).
Reliability of the ENMs designed using an opportunistic sampling strategy The orographic di culties (elevation range 1200-1370 m) in the burned area led us to propose a type of sampling strategy that differs from traditional methodology (i.e. systematic sampling, Hirzel and Guisan 2002). Here, we carried out opportunistic sampling making use of existing roads and tracks ) during 12 eldwork days, achieving the objective of representing the environmental and climatic variability of the study area with the least possible number of points based on the distance between them (800 m). We tted the Local2 models to the sample data, obtaining statistically reliable (Table S1) and ecologically consistent results, corresponding to the hypothesis.
Opportunistic sampling also meant a reduction in efforts and economic costs due to the following: (1) The increased performance obtained by the reduction in time required to access to the sample plots, located close to existing roads and tracks.
(2) The use of an open-source application on a mobile phone platform (Epicollect5) for georeferencing the sample plots. It should be noted that the georeferencing accuracy will depend on the capacity of the mobile phone used, and its geolocation. In this case, the geospatial precision was 5 to 10 m in all sample plots.
(3) The small number of people required to carry out the eld work. The ndings suggest that opportunistic sampling could be used to develop ENMs as an alternative to traditional systematic grid design, at least when the budget and/or time available to complete the sampling process are limited ). In addition, this sample strategy may be effective for validating post-re regeneration maps based on remote sensing (i.e. Gitas et al. 2012, Fernández-Guisuraga et al. 2019.

Implications for post-re management
We designed a decision support scheme (Table S3) based on ENMs outputs (Figs. 1, 2, S2, S3), the binary models overlapping maps for each species depending on the scenario (current, RCP 4.5 and RCP 8.5) (Fig. S4) and prediction of forest stands by combining the binary ENMs for each species for the current and RCP 4.5 scenario (Fig. S5), in order to show how these models could help to support strategic post mega re restoration-related decisions.
This theoretical exercise shows how the use of regional and local ENMs in different scenarios can help managers to make decisions regarding the restoration strategy (active, passive) and the management strategy in the case of active restoration actions: silviculture (thinning, shrub clearing, wild re prevention), grazing management, wildlife management, monitoring and afforestation. For example (Fig. S5), in the current scenario all species are present and biomass accumulation is ongoing (Madrigal et al. 2017) showing a tendency towards mixed forest of hardwoods and conifers. Therefore the restoration strategy suggests passive (no action) (Castro et al. 2011) or active management (Moreira et al. 2020) in order to dosify competence and to reduce the vulnerability of the area to wild res Alfaro-Sánchez et al. 2014). In the RCP 4.5 scenario, the intermediate probability of the presence of Q. faginea and P. pinaster and the absence of Q. ilex at the regional scale ( Table 2) together suggest that these species should be focused on in developing an active management strategy including actions to promote resilience after wild res. Finally, in the RCP 8.5 scenario, the tree vegetation cover is expected to decrease until disappearing (regional model) and the probability of new mega res will be less important due to the absence of available fuels (Resco de Dios 2020).
Observing the different potential distributions of the overlapping species (Fig. S4), the models predict that in a future scenario the stand will evolve towards a mixed forest with a greater presence of hardwood trees relative to conifers (Fig. S5). Under these conditions, the restoration strategy should focus on active management to promote mixed forests (Bravo et al. 2019), active and passive strategies to promote biodiversity (Leverkus et al. 2019) and the need to monitor regenerated species to decide on potential local changes in species adapted to the new climatic niches in the near future (Del Campo et al. 2020).

Conclusions
To our knowledge, this is the rst study proposing the use of ENMs as a tool for post-mega re management of forest land. Our research approach produced accurate ENMs (spatial resolution of metres) through an opportunistic sampling strategy approach, with optimization of resources, representing a more affordable and dynamic working option for research groups. Furthermore, we demonstrated that, in addition to carrying out eldwork shortly after the mega re, it is also important to carry out eldwork at a later date. The latter considers a more advanced phase in vegetation regeneration dynamics, implying an improvement in representation of the suitable area of distribution of species by the models in the current scenario.
Our study ndings show that calibrating ENMs at different scales (local and regional) provides useful information for management decision-making, and that the models can be used as an objective tool in restoration plans or in other post-mega re management actions. We propose that post-re restoration plans should contemplate the effects of climate change on regeneration dynamics. Thus, we suggest the use of ENMs as a proactive support tool for managers, as these models provide greater dynamism to restoration plans (despite some recognised limitations), in addition to other types of management actions implemented after mega res.

Declarations
Funding This study was supported by the Spanish Ministry of Science and Innovation through the projects VIS4FIRE (Integrated Vulnerability of Forest Systems to Wild re: Implications on Forest Management Tools, RTA2017-00042-C05-01) and FORESTCHANGE (In uence of natural disturbance regimes and management on forests dynamics, structure and carbon balance, AGL2016-76769-C2-1-R). VIS4FIRE is co-funded by the EU through the FEDER program. C.CG was nanced by the Ministry of Science and Innovation, the European Social Fund and the State Research Agency through the predoctoral aid programme (PRE2018).
Con icts of interest/Competing interests The authors declare no con ict of interest.
Ethics approval Not applicable Consent to participate Not applicable

Consent for publication Not applicable
Availability of data and material The databases required to calibrate the models were provided by the regional government of Castilla-La Mancha solely for the purpose of this work Tables   Table 1   Local1 vs Local2 Local1 vs Regional Local2   Figure 1 (a) Regional ENM for Pinus pinaster (b) Local1 ENM for P. pinaster (short-term post-re regeneration) (c) Local2 for P. pinaster ENM (long-term post-re regeneration). Suitable distribution index under the current scenario. All of the models were rescaled from 0 to 100 potential suitability values Suitable distribution index for ENMs developed for maritime pine (P. pinaster) under the different future scenarios (RCP 4.5 and RCP 8.5). All of the models were rescaled from 0 to 100 potential suitability values. Regional, Local1 (short-term post-re regeneration) and Local2 (long-term post-re regeneration) are shown. The regional model under the RCP 8.5 scenario did not predict the presence of the species

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. renamed769fb.pdf