Climate- and fire-smart landscape scenarios call for redesigning protection regimes to achieve multiple management goals

Journal of Environmental Management 322 (2022) 116045 Available online 5 September 2022 0301-4797/© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Research article Climateand fire-smart landscape scenarios call for redesigning protection regimes to achieve multiple management goals Miguel C´ anibe Iglesias a , b , c , * , Virgilio Hermoso d , e , Jo˜ ao C. Campos f , g , Cl´ audia Carvalho-Santos h , Paulo M. Fernandes i , Teresa R. Freitas i , Jo˜ ao P. Honrado f , j , k , Jo˜ ao A. Santos i , ˆ Angelo Sil a , b , f , i , k , Adri´ an Regos c , d , f , k , Jo˜ ao C. Azevedo a , b a Centro de Investigaç˜ ao de Montanha (CIMO), Instituto Polit´ ecnico de Bragança, Campus de Santa Apol´ onia, 5300-253, Bragança, Portugal b Laborat´ orio Associado para a Sustentabilidade e Tecnologia em Regi˜ oes de Montanha (SusTEC), Instituto Polit´ ecnico de Bragança, Campus de Santa Apol´ onia, 5300253, Bragança, Portugal c Departamento de Zooloxía, Xen´ etica e Antropoloxía Física, Universidade de Santiago de Compostela, 15782, Santiago de Compostela, Spain d Centre de Ci` encia i Tecnologia Forestal de Catalunya (CTFC), Ctra. Sant Llorenç de Morunys, km2. 25280, Solsona, Lleida, Spain e Departamento de Biología Vegetal y Ecología, Universidad de Sevilla, 41012, Sevilla, Spain f InBIO/CIBIO - Centro de Investigac˜ ao em Biodiversidade e Recursos Gen´ eticos, Campus Agr´ ario de Vair˜ ao, Rua Padre Armando Quintas, nº 7, 4485-661, Vair˜ ao, Portugal g CICGE – Centro de Investigaç˜ ao em Ciˆ encias Geo-Espaciais, Faculty of Sciences, University of Porto, Alameda do Monte da Virgem, 4430-146, Vila Nova de Gaia, Portugal h Centre of Molecular and Environmental Biology (CBMA) & Institute for Bio-Sustainability (IB-S), University of Minho, 4710–057, Braga, Portugal i CITAB - Centro de Investigaç˜ ao e de Tecnologias Agro-Ambientais e Biol´ ogicas, Universidade de Tr´ as-os-Montes e Alto Douro, 5001-801, Vila Real, Portugal j Departamento de Biologia, Faculdade de Ciˆ encias da Universidade do Porto, Rua do Campo Alegre, S/N, Edifício FC4, 4169-007, Porto, Portugal k BIOPOLIS Program in Genomics, Biodiversity and Land Planning, CIBIO, Campus de Vair˜ ao, 4485-661, Vair˜ ao, Portugal ARTICLE INFO Keywords: Conservation planning Fire management Biodiversity conservation Ecosystem services Fire-smart Climate-smart ABSTRACT Integrated management of biodiversity and ecosystem services (ES) in heterogeneous landscapes requires considering the potential trade-offs between conflicting objectives. The UNESCO’s Biosphere Reserve zoning scheme is a suitable context to address these trade-offs by considering multiple management zones that aim to minimise conflicts between management objectives. Moreover, in Mediterranean ecosystems, management and planning also needs to consider drivers of landscape dynamics such as wildfires and traditional farming and forestry practices that have historically shaped landscapes and the biodiversity they host. In this study, we applied a conservation planning approach to prioritise the allocation of management zones under future landscape and climate scenarios. We tested different landscape management scenarios reflecting the outcomes of climate-smart and fire-smart policies. We projected the expected landscape dynamics and associated changes on the distribution of 207 vertebrate species, 4 ES and fire hazard under each scenario. We used Marxan with Zones to allocate three management zones, replicating the Biosphere Reserves zoning scheme (“Core area”, “Buffer zone” and “Transition area”) to address the various management objectives within the Biosphere Reserve. Our results show that to promote ES supply and biodiversity conservation, while also minimising fire hazard, the reserve will need to: i) Redefine its zoning, especially regarding Core Areas, which need a considerable expansion to help mitigate changes in biodiversity and accommodate ES supply under expected changes in climate and species distribution. ii) Revisit current management policies that will result in encroached landscapes prone to high intensity, uncontrollable wildfires with the potential to heavily damage ecosystems and compromise the supply of ES. Our results support that both climateand fire-smart policies in the Meseta Ib´ erica can help develop multifunctional landscapes that help mitigate and adapt to climate change and ensure the best possible maintenance of biodiversity and ES supply under uncertain future climate conditions. * Corresponding author. Centro de Investigaç˜ ao de Montanha, Instituto Polit´ ecnico de Bragança, Campus de Santa Apol´ onia, 5300-253, Bragança, Portugal. E-mail addresses: [email protected] (M.C. Iglesias), [email protected] (V. Hermoso), [email protected] (J.C. Campos), c.carvalho.santos@bio. uminho.pt (C. Carvalho-Santos), [email protected] (P.M. Fernandes), [email protected] (T.R. Freitas), [email protected] (J.P. Honrado), [email protected] (J.A. Santos), [email protected] (ˆ A. Sil), [email protected] (A. Regos), [email protected] (J.C. Azevedo). Contents lists available at ScienceDirect Journal of Environmental Management journal homepage: www.elsevier.com/locate/jenvman https://doi.org/10.1016/j.jenvman.2022.116045 Received 30 May 2022; Received in revised form 22 July 2022; Accepted 17 August 2022 Journal of Environmental Management 322 (2022) 116045 2 1. Introduction The loss and degradation of ecosystems is leading to a global decline of biodiversity at rates 100 times higher than the background extinction rates for many taxa (Ceballos et al., 2015; IPBES, 2019). Furthermore, the increasing demand for food, water, raw materials and energy, mostly driven by an increasing growth of human population, is leading to the unsustainable use of renewable and non-renewable resources, further degrading ecosystems. Despite the remarkable conservation efforts in the last decades, biodiversity and ecosystem condition are still declining due mainly to the impacts of regional land-use and global climate changes that lead to shifts in disturbances regimes (IPBES, 2019). All these pressures compromise the persistence of biodiversity (Habel et al., 2019) and consequently the supply of ecosystem services (ES) (Cabral et al., 2021). Climate-smart management originated in agricultural systems (Scherr et al., 2012) and was later extended to forest management (Bowditch et al., 2020). Climate-smart landscape management aims to help mitigate the effect of climate change by increasing carbon stocks and sequestration rates via measures such as direct afforestation and rewilding initiatives (Perino et al., 2019). However, rewilding and afforestation can reduce landscape heterogeneity and increase fuel load and connectivity in the landscape, which in turn reduce ecosystem resilience and adaptability (Holl and Brancalion, 2020) and consequently increase wildfire hazard (Hermoso et al., 2021). This therefore requires the implementation of fire-smart management solutions with the goal of building fire resilient landscapes alongside climate-smart policy and practice while maintaining high levels of biodiversity and ecosystem services delivery (Hirsch et al., 2001). Fire-smart management is especially relevant in Mediterranean regions, where wildfires are a key driver of landscape dynamics (Lloret et al., 2002). These regions are currently experiencing an increased risk of large and high-intensity fires due to the combined effect of climate warming, longer drought periods, and long-standing land abandonment or afforestation processes (Moreira et al., 2020) that increase fuel load and connectivity at site and landscape levels (e.g., Fernandes et al., 2014). Additionally, traditional fire management policies based on fire exclusion have increased the risk of catastrophic fires under extreme weather conditions (Pausas and Fern´ andez-Mu˜ noz, 2012). Fire-smart management can therefore play a key role in reducing fire hazard by promoting fire-smart landscapes in Mediterranean regions (Pais et al., 2020). Under favourable conditions, fire-smart management can even create opportunities for fire to provide benefits to some species (Regos et al., 2018) or increase fire resilience by reducing fuel load and connectivity (e.g., prescribed burning; Fernandes et al., 2013). Another benefit of fire-smart management is its contribution to maintaining landscape heterogeneity, minimising the negative impacts of wildfires on biodiversity and ES. The integration of ES in landscape management is vital for the development of more integrative conservation frameworks and financing mechanisms such as the “Reducing Emissions from Deforestation and forest Degradation” (REDD+) or the EU Green Infrastructure Strategy (Benedict and McMahon, 2002). Despite efforts, the amount of land required to fulfil human needs continues to grow (Foley et al., 2011) whereas the ecosystem’s ecological integrity continues to decline (Plumptre et al., 2021). This highlights the need for landscape planning and management approaches that aim to protect land while integrating conservation and development, avoiding, or at least minimising, the conflicts that are likely to arise in the presence of limited resources for different, and often conflicting, objectives. The successful implementation of holistic landscape management approaches relies on our capacity to address trade-offs and synergies among different ES and between ES and biodiversity conservation (Carvalho-Santos et al., 2016; Mor´ an-Ord´ o˜ nez et al., 2017; Sil et al., 2016). Studies on the spatial relation between biodiversity and ES supply have shown that synergies between conservation and socioeconomic development objectives exist and are a great opportunity for integrated landscape-level planning and management (Chan et al., 2006; Egoh et al., 2009; Nelson et al., 2009; Ramel et al., 2020). However, the relationship between biodiversity and ES is complex. Biodiversity provides and regulates ecosystem processes responsible for the supply of ES or it can be considered as an asset (Mace et al., 2012). This is also highly dependent on the socioeconomic and biophysical context of the regions of interest, as well as on the scale of the analysis (Duncan et al., 2015). In recent studies, spatial trade-offs and synergies between ES and biodiversity in systematic conservation planning have been addressed using management zones (Barbosa et al., 2019; Hermoso et al., 2018). These zones aim to achieve ES supply and biodiversity conservation goals simultaneously, enhancing co-benefits between objectives, minimising potential trade-offs. Some of these management zones can be designed to simultaneously address compatible objectives, such as biodiversity conservation and carbon storage, and others to allow uses that are not compatible, such as conservation and timber production (Lanzas et al., 2019) or fisheries and conservation (Beger et al., 2015). This multi-zoning approach allows for more flexibility in planning, securing larger amounts of resources while minimising conflicts between objectives when compared to a single management zone (Hermoso et al., 2018). UNESCO’s Biosphere Reserves (BR) offer an innovative and integrated management model to preserve biodiversity along with the sustainable use of natural resources and research (UNESCO, 2021). To efficiently achieve these multiple goals, BRs have popularised a flexible multi-zoning scheme based on three management zones (UNESCO, 2017a): i) Core Areas, with stricter biodiversity conservation goals usually encompassing areas already included in protected areas with higher conservation value; ii) Buffer Zones, aiming to buffer and connect Core Areas while allowing for traditional activities compatible with conservation; and iii) Transition Areas, where sustainable resource management is promoted. Although zoning should be case-specific, this scheme presents a good starting point for holistic landscape and conservation management where climateand fire-smart policies can be accommodated while avoiding or minimising critical trade-offs between objectives. In this study, we aimed to prioritise the allocation of different management zones within the Meseta Ib´ erica Transboundary Biosphere Reserve (NW Iberian Peninsula) on the Portugal-Spain border (UNESCO, 2017b) to achieve multiple objectives for biodiversity conservation and ES supply while minimising wildfire hazard under different landscape management scenarios. The area features a Mediterranean landscape, affected by land use change and wildfires. We assessed contrasting landscape management scenarios (from climate-to fire-smart strategies based on likely outcomes of forestry and agricultural policies) and their potential effects on biodiversity, ES supply and fire regime. 2. Material and methods 2.1. Study area Our study was conducted in the Meseta Ib´ erica Transboundary Biosphere Reserve (Fig. 1) located in the north-west Iberian Peninsula. This BR, designated in 2015, has a total extent of 11,326 km 2 (UNESCO, 2017b) including territories from both Portugal (58% of the BR) and Spain (42%; Trillo Santamaría and Paül Carril, 2018). The BR comprises 12 Portuguese municipalities in the district of Bragança plus the municipality of Figueira de Castelo Rodrigo in the district of Guarda. On the Spanish side the BR includes 75 municipalities, 48 in the province of Zamora and 27 in the province of Salamanca, all of them in the Castilla y Le´ on autonomous community (ZASNET, 2021). The BR follows the conventional Man and Biosphere structure comprised of Core areas (areas of higher level of protection within protected areas such as Arribes del Duero, Douro International, Montesinho and Lago de M.C. Iglesias et al. Journal of Environmental Management 322 (2022) 116045 3 Sanabria y Sierras de Segundera y Porto natural parks, as well as the Regional Natural Park Vale do Tua and the Sierra de la Culebra Site of Community Importance), Buffer zones (areas of lower level of protection inside protected areas and Natura 2000 sites) and Transition areas (the remaining areas) (Santamaría and Carril, 2018). Currently, 1,064 km 2 (9%) are allocated to the Core area, 4,203 km 2 (36%) to the Buffer zone, and 6,325 km 2 , (55%) to the Transition area (Palliwoda et al., 2021). The landscape of the BR is diverse and heterogeneous. Altitude Fig. 1. Map of the Meseta Ib´ erica Transboundary Biosphere Reserve showing the current distribution of management zones. Top-right corner shows the location of the study area within the Iberian Peninsula. Data supplied by ZASNET European Grouping of Territorial Cooperation. M.C. Iglesias et al. Journal of Environmental Management 322 (2022) 116045 4 ranges from 100 to over 2,127 m. a.s.l. (Xunta de Galicia, 2021). Climate is mostly Mediterranean with dry-warm summers and wet winters. Mean annual precipitation varies from 500 to 1200 mm following an altitudinal gradient and presents a strong seasonality (Deitch et al., 2017; Santos and Belo-Pereira, 2022). Land cover is mainly shrubland, farmland, and forests (Azevedo, 2012; Sil et al., 2017). The main types of forests are maritime pine (Pinus pinaster) plantations, deciduous woodlands dominated by Pyrenean oak (Quercus pyrenaica) and evergreen woodlands dominated by holm oak (Quercus ilex) and cork oak (Quercus suber) (Azevedo, 2012). The BR also harbours high levels of species richness, including a high number of invertebrates, plants and around 250 species of vertebrates, among which flagship species such as the Iberian wolf (Canis lupus signatus), Egyptian vulture (Neophron pernocpterus), and iberian endemisms such as Bocage’s wall lizard (Podarcis bocagei) or Seoane’s viper (Vipera seoanei) (UNESCO, 2017b). The area hosts a human population of around 300,000 inhabitants (UNESCO, 2017b). Depopulation and ageing cause high rates of land abandonment (Sil et al., 2016). The Mediterranean type of climate in the region, characterized by wet mild winters and dry warm summers, together with the landscape changes derived from land abandonment (vegetation encroachment and afforestation), has led to an increased risk for severe wildfires (Sil et al., 2019). 2.2. Conceptual framework The research followed a systematic conservation planning approach based on spatial modelling of species distribution and supply of ES (Fig. 2) under historical (hereafter 2005 scenario) and future conditions (2050). We considered four landscape management options and four climate models projected under two climate change scenarios (Representative Concentration Pathways RCP 4.5 and 8.5), which were used to evaluate the impact of climate change on species distribution and ES. Spatial projections of potential distributions for 207 species and four ES were jointly incorporated into Marxan with Zones (Watts et al., 2009), a decision support tool that has been successfully used in spatial conservation planning integrating biodiversity and ES (Adams et al., 2016; Barbosa et al., 2019; Lanzas et al., 2019). Marxan with Zones prioritises the allocation of management zones with user-defined roles (i.e., which features are allocated to a management zone and which are not) in a flexible way that allows conflicting uses to be managed separately in different zones. In addition, we used potential fire intensity (measured by fireline intensity) as a penalty in the spatial prioritisation exercise, thus favouring selection of planning units expected to burn with lower intensities to achieve biodiversity and ES supply targets. With the resulting distribution of management zones, we compared the performance of the different scenarios in regard to species and ES coverage. We also compared our zone distributions with the zone distribution in the existing management plan of the BR. All Marxan analyses were conducted based on data for features and penalties gathered in a grid of 1 km 2 cells (hereafter “Planning Units”) covering the BR area. Fig. 2. Schematic workflow of the study. M.C. Iglesias et al. Journal of Environmental Management 322 (2022) 116045 5 2.3. Landscape management scenarios We considered four landscape management scenarios projected for the study area for the year 2050 (obtained from Campos et al., 2022). These scenarios depict future states of the landscape based on the implementation of forest and agricultural management options (“Afforestation” and “BAU” scenarios) and fire-smart (“FarmReturn” and “AgroforestRe”) policies supported by climate. Scenarios were defined according to the main trends of landscape change identified in 5 periods (1990–2018; 1990–2000; 2000–2006; 2006–2012 and 2012–2018), selecting the most representative trend for the storyline of each scenario. The landscape management scenarios are defined according to the following storylines (see details in Campos et al., 2022): ●Afforestation. Defined by forest expansion resulting from afforestation as identified in the landscape change trends from 1990 to 2000. This scenario was used to account for changes associated with increasing wood and bioenergy demand as well as for climate change mitigation. Main transitions are a strong increase in forest areas through conversion of semi-natural areas (shrubland and grassland) and an increase of deciduous/broadleaved species through natural succession and active planting mainly in shrubland (Table S1.2). ● Business as usual (BAU). Defined by land abandonment, following historical and current trends of abandonment in the area (Azevedo et al., 2011), based on 1990–2018 trends. This scenario results in a landscape dominated by shrublands growing in former agro-pastoral areas. ●FarmReturn. Defined by the support of agricultural policies promoting sustainable low maintenance farming and reverting land abandonment tendencies (i.e., European Union’s Common Agricultural Policy), contributing to biodiversity conservation and developing low fire hazard landscapes (Moreira & Pe’er, 2018), based on 2006–2012 trends. It is characterized by an increase in farmland at the expense of semi-natural (shrubland and grassland) areas. ● AgroforestryReturn (AgroforestRe). A scenario where the support of agricultural and agroforestry policies will create a potentially more fire-resilient and fire-resistant landscape with lower fuel load and connectivity, based on 2006–2012 and 2012–2018 data. The main trends are a moderate replacement of semi-natural areas and coniferous forest by croplands, and a strong replacement of deciduous forest, shrubland and grassland by agroforestry areas (e.g. sweet chestnut groves). Projections of land cover under each landscape management scenario were built using the Scenario Generator of InVEST (Sharp et al., 2020) with CORINE Land Cover (CLC; Copernicus, 2020) data for 2018 as the baseline Land Cover. CLC data were grouped into 10 broader land cover classes for analysis, namely urban, agriculture, grassland, agroforestry, forests (deciduous, coniferous and mixed), shrubland, water and others (Table S1.1). The trends defined above were used to produce landscape transition matrices for all scenarios. To account for the stochasticity of landscape dynamics, 10 simulations were run for each scenario. Finally, land-use/cover data projected under each of these scenarios were used as input to predict changes in species distribution, ES supply and fire hazard. 2.4. Biodiversity data For biodiversity predictive mapping, two sets of species distribution models (SDMs) were used to account for the joint effects of climate and land-use change: (1) SDMs based only on climate predictors, obtained from Campos et al. (2021), and (2) SDMs based on Land Use/Land Cover (LULC) and topographic variables, obtained from Campos et al. (2022). Both sets were built from presence/absence data for 168 birds, 24 reptiles and 15 amphibians from national atlases at 10-km resolution for the whole Iberian Peninsula to characterise the ecological niche of the species (see Titeux et al., 2017). Individual projections were obtained using 6 modelling algorithms and 10 replicates to account for modelling stochasticity and were then used to compute ensemble models considering AUC values as model weights for each future management scenario. These ensemble models were then downscaled and projected at 1-km resolution to the extent of the BR (Bombi & D’Amen, 2012) and reclassified into binary presence/absence maps using ROC optimised thresholds (Thuiller et al., 2009). Habitat models’ projections were obtained for the 2005 scenario and the four landscape management scenarios (2050). To deal with the uncertainty of climate change, we considered four widely used models climate models (IPSL-IPSL-CM5A-MR, ICHEC-EC-EARTH, MPI-M-MPI-ESM-LR and CNRM-CERFACS-CNRM-CM5) from the European Coordinated Downscaling Experiment (EURO-CORDEX; Jacob et al., 2020) under two climate scenarios (RCPs 4.5 and 8.5). RCP 4.5 corresponds to an intermediate anthropogenic radiative forcing of the climate system, with a mid-century peak in greenhouse gas emissions and a subsequent decline thereafter. RCP 8.5 is a fossil-fuel emissions intensive scenario, commonly considered the worst-case scenario (van Vuuren et al., 2011). For climate projections, we used average predictions of the four climate models under each RCP scenario (see Campos et al., 2021). All the SDMs procedures were performed using the “biomod2” R package (Thuiller et al., 2009). Only locations where species presence was predicted by both climatic and habitat models were used for this study. Complete modelling details are available in Campos et al. (2022). 2.5. Ecosystem services Four ES were selected covering the three highest levels of the Common International Classification of Ecosystem Services version 5.1 (CICES v5.1) categories: regulation and maintenance, provisioning, and cultural. The selection was limited to those ES potentially affected by land-use changes, since the entire research framework relies on landscape change scenarios. 2.5.1. Provisioning ES: cultivated terrestrial plants We used the amount of agricultural surface as a surrogate for provisioning services that depend on this type of land cover. This represents ES in the “Cultivated terrestrial plants for nutrition, materials or energy” group in CICES v5.1 and it was chosen since CORINE land cover maps do not differentiate between particular end uses of crops (i.e., nutrition, materials, energy). We used 1-km resolution LULC maps to identify planning units classified as agriculture. 2.5.2. Regulation and maintenance ES: climate regulation We used the InVEST Carbon Storage and Sequestration module (Sharp et al., 2020) to assess the dynamics of the climate regulation ecosystem service (CRES) in the 2005 and 2050 landscape scenarios. CRES is the contribution of terrestrial systems to regulate the concentration of greenhouse gases in the atmosphere (Haines-Young and Potschin, 2018). We used carbon sequestration rate (Mg C ha −1 yr −1 ) as a proxy of the capacity of ecosystems and landscapes to supply CRES. The InVEST module was fed with data on carbon stocks based on previous studies in the area [see Campos et al. (2022) and Sil et al. (2017) for a complete description]. Carbon stocks were estimated for seven major land cover classes (agriculture, agroforestry, deciduous, coniferous and mixed forest, and semi-natural grassland and shrubland) and four carbon pools (aboveground, belowground biomass, soil organic carbon, and dead organic matter). The amount of carbon gain (sequestration) or loss (emission) was computed as the difference between stocks in each pixel on two consecutive dates of the two periods (1990–2020 and 2020–2050). Raster maps of carbon sequestered or emitted in each period were divided by the number of years in that period to obtain maps of annual estimates to be used as inputs in Marxan with Zones. M.C. Iglesias et al. Journal of Environmental Management 322 (2022) 116045 6 2.5.3. Regulation and maintenance ES: soil erosion control Soil erosion control was estimated according to Guerra et al. (2014), who measured avoided soil erosion due to the effect of vegetation, providing the actual ecosystem service. The approach is based on the Revised Universal Soil Loss Equation (RUSLE) which estimates annual soil loss through the product of rainfall erosivity (R), soil erodibility (K), slope length and steepness (LS), cover-management (C) and conservation practices (P) factors, the latter not considered due to the absence of spatial data (Eq. (1)): A=R×K×LS ×C(1) It differs from the traditional application of RUSLE in the computation of erosion made under two conditions: i) the structural impact, i.e., the erosion that would occur if vegetation was absent (Eq. (2)): S=R×K×LS (2) and ii) the actual soil loss (Eq. (1)). Soil erosion control ES was therefore estimated by subtracting the structural impact (Eq. (2)) from the actual soil loss (Eq. (1)). Control of soil erosion was calculated for the four proposed land management scenarios and respective replicates where the C factor, obtained from Pimenta (1998), was used in the reclassification of CORINE land cover map classes. Rainfall erosivity (R), soil erodibility (K) and slope length and steepness (LS) were obtained from the European Soil Data Centre (ESDAC) (Panagos et al., 2014; Panagos, Ballabio, et al., 2015; Panagos et al., 2015b). All rates were normalized in a 0 to 1 scale to be used as input in Marxan with Zones. 2.5.4. Cultural ES: recreation Recreation potential was modelled following the ESTIMAP model for nature-based recreation (NBR) (Zulian et al., 2013). This model uses advanced multiple layers lookup tables (advanced LUT) to assign ES scores to land units based on cross-tabulation from different input layers. NBR potential combines ecosystem-based potential to provide NBR and distance to NBR potential. Ecosystem-based potential combines three sources of information into a single layer: i) Suitability of each LULC class to support recreation based on a score from 0 to 1 representing the suitability of each LULC class to support these activities (Vallecillo et al., 2019) (Table S2.1); ii) Area-based conservation measures according to which we assigned an additional score to Protected Areas and Natura 2000 sites considering their attractiveness to people when deciding where to spend their free-time considering conservation areas (World Database of Protected Areas) differently - Natural and Regional parks were assigned a score of 1 while Natura 2000 sites (Special Areas for Conservation, Special Protection Areas, and Sites of Community Importance) were assigned a score of 0.8; and iii) Water masses, to which we assigned a score of 1 to include important fluvial beaches and other inland water elements used for recreation collected in the European environmental Agency (EEA) state of bathing water database. The 3 components represented in raster layers were summed up obtaining a 0–3 layer of scores, subsequently normalized to the 0–1 range, representing the Recreation Potential Index (RPI). RPI was then classified in “Low”, “Medium” and “High” classes using the 33 and 66 percentiles. Distance to NBR potential indicates accessibility and remoteness of areas with recreation potential. Both metrics are based on the Euclidean distance (in km) from roads (OSM contributors, 2021) and urban settlements, respectively. These measures were cross-tabulated to obtain the distance matrix (Table S2.2). The NBR provision layer was obtained by cross-tabulating ecosystem-based potential and distance components according to parameters in Table S2.3. Final NBR scores ranged from 1 to 9 (Table S2.4). Of these, we only used high recreation provision classes (7, 8 and 9) as inputs in Marxan with Zones. 2.6. Landscape fire hazard We applied the FlamMap module from the FlamMap5 (v5) fire mapping and analysis system model (Finney et al., 2015) to assess the effect of landscape change on fire behaviour in a spatially explicit manner, and derive information on the potential fire hazard in the study area in past (using CLC 2006) and in future landscape scenarios (2050). This information was used as costs in Marxan with Zones. We assumed fireline intensity (kW m −1 ) as the descriptor of potential fire hazard. To express the resistance to control of a wildfire, fireline intensity outputs were reclassified to be used as inputs in Marxan with Zones according to a standard fire danger classification (Alexander and Lanoville, 1989): Class 1: Low (<500 kW/m); Class 2: Moderate (500–2000 kW/m); Class 3: High (2000–4000 kW/m); Class 4: Very High (4000–10000 kW/m); and Class 5: Extreme (>10000 kW/m). In FlamMap, raster layers of fuels and topographic conditions and tabular data, and several built-in parameters were used to set fuel moisture and weather variables (Table S3.1). Fuel models were allocated based on the correspondence between land cover classes and custom fuel models for Portugal (Fernandes et al., 2009) (Table S3.2). Canopy cover data for each forest type was based on previous work within the study area (Azevedo et al., 2011) (Table S3.3). For canopy fuel variables (stand height, canopy base height and canopy bulk density) we used data from Botequim et al. (2019) collected data in a Mediterranean climate area of SW Spain for P. pinaster and Q. pyrenaica in pure and mixed stands. Fire behaviour was simulated under severe (dry and windy) weather conditions, expected to be more common under climate change (Table S3.4). The fuel moisture content of surface fuels (dead and live) and foliar moisture content (FMC) of canopy fuels was set based on typical conditions (Fernandes, 2009). Wind speed is representative of wind gusts in active crown fires (Cruz and Alexander, 2019). Alignment between wind and slope was assumed for all simulations to depict maximum fire behaviour potential. All raster files for fuels and terrain were prepared at 100-m spatial resolution using GIS functions. 2.7. Spatial prioritisation of management zones We used Marxan with Zones (Watts et al., 2009) to prioritise the spatial allocation of the three management zones within the BR. Marxan with Zones uses data on the spatial distribution of conservation features (in our case species and ES) and costs (in our case fire intensity), to identify the most suitable allocation of management zones that allow achieving user-defined representation targets for the features at a minimum cost. Marxan with Zones also allows specifying the spatial aggregation within management zones and the spatial relationship between management zones. The mathematical problem that we addressed was therefore: minimise ∑ m i=1 ∑ p k=1 cikxik +b∑ m i1=1 ∑ m i2=1 ∑ p k1=1 ∑ p k2=1 cvi1,i2,k1,k2xi1,k1xi2,k2(3) subject to ∑ m i=1 ∑ p k=1 aijxik ≥tjk ∀j(4) where, cik is the cost of planning unit i if allocated under zone k; xik is a control variable that determines whether planning unit i has been allocated under zone k (1) or not (0); cvi1,i2,k1,k2 is the connectivity penalty for including only one of the pair of planning units i1, i2; xi1,k1 and xi2,k2 are control variables that take values of 1 when the planning unit i1 or i2 is included in the solution or 0 otherwise; b, or boundary length modifier (BLM), is a weight applied to the connectivity penalty used to aggregate planning units in space or determine the spatial structure of zones; aij is the contribution of planning unit i to the achievement of targets for feature j; and tjk is the representation target desired for each j feature under their respective zone k. M.C. Iglesias et al. Journal of Environmental Management 322 (2022) 116045 7 2.7.1. Representation targets We set an overall representation target of 200 km 2 for each species under all planning scenarios and time horizons. We selected this representation target to ensure an adequate representation of the rarest species in the study area, most in need of conservation action, while avoiding over-representing the most common ones. The target we set represents the full distribution of the 20 (about 10% of all species considered) rarest species, while only a small proportion (2–3%) of the distribution of the most common species. For the rare species that do not reach 200 km 2 in the area, we set their total distribution as the target. In the future scenarios, we used exactly the same targets although expecting that this would lead to some representation targets being impossible to achieve in the case of species that heavily decline or disappear from the area in the future. Last, for species that appear in the area only in future scenarios, we set new targets following the criteria above. We believe that this is a good way to identify species turnover within Marxan with Zones in scenarios of uncertainty. Regarding ES, we aimed to explore the maintenance of high levels of ES in the future while avoiding conflicts. Since higher amounts of ES would make some of the targets impossible to achieve (i.e., asking for more agriculture than is present in a given scenario), we set the targets to 70% of 2005 supply. Targets for the future scenarios were kept constant according to the absolute amounts of ES demanded for the historic scenario (Table 1). We replicated the same types of zones of the UNESCO zoning scheme (Transition, Buffer and Core areas) in our analyses. To distribute the above-mentioned targets across these three management zones, we first evaluated the potential relationships between the different features to look for potential trade-offs or opportunities to foster co-benefits, following recommendations in Hermoso et al. (2018), Lanzas et al. (2019) and Sil et al. (2016). Based on knowledge from the authors in the study area, we identified cultivated terrestrial plants as a conflicting ES that can negatively impact carbon sequestration, soil erosion and habitat for some species, and thus tried to allocate these to different management zones (Table 2). For species, we distributed the overall targets above according to their distribution ranges to ensure that species with distributions below 200 occurrences achieve their targets within the Core area. For species above 200 occurrences, we split their targets between the Core area and Buffer zone and for open habitat and generalist species we allowed a portion of their targets to be met in the Transition area. For ecosystem services, the overall 70% target was partitioned in fractions of 10, 25 and 35% and distributed according to their potential impacts on conservation purposes (Table 1). We allowed a small proportion of ES targets to be met in zones where they might cause conflicts with other objectives because we assumed that adequate management of the BR can allow for small portions of incompatible ES to coexist within the same management zone. 2.7.2. Spatial configuration of management zones Marxan with Zones allows specifying the degree of spatial aggregation within management zones as well as the spatial arrangement among zones through weighting factors in the objective function: the Boundary Length Modifier (BLM) and the weights in the “zoneboundary” file. We used an overall BLM value of 1 and calibrated the “zoneboundary” file parameters following Serra et al. (2020) to ensure that the Buffer zone buffers the Core area, and the Transition and Core areas are not connected. 2.7.3. Spatial penalties In our research framework we penalised the selection of planning units with a high fire intensity risk, assuming that whenever fire suppression difficulty was rated high or very high, potential fire damage is higher e (high penalty), while in areas where fire suppression difficulty was rated low to moderate fire damage is lower (low penalty). This cost was equally applied to all zones. 2.7.4. Feature penalties Failing to achieve representation targets results in penalties. As such, target achievement is encouraged in the optimisation procedure. The feature penalties are weighted by a Feature Penalty Factor (FPF) in the Marxan objective function, so that high SPF results in all targets achieved, while low SPF can lead to some features not meeting their targets. To ensure that representation targets were always met, we used a FPF of 10 for all features except for rare species and ES, for which we used a FPF of 100. With the specifications detailed above, we ran Marxan with Zones 100 times (10 million iterations in each individual run) for each of the 81 simulations (4 landscape management scenarios X 10 replicates X 2 RCPs +1 for the “2005” scenario) and kept the best solution over those runs for subsequent comparative analyses across scenarios. 2.8. Analysis of Marxan with zones solutions We compared the solutions obtained under each scenario by recording the extent and the mean potential fire intensity within each zone. We also compared the amount of ES and species distributions covered within each management zone across landscape management scenarios and RCPs. Finally, we used the Jaccard index to compare the spatial allocation of management zones derived from Marxan with Zones with the configuration of the zoning currently implemented in Meseta Ib´ erica (Eq. (5)). The Jaccard index measures the spatial overlap between the distribution of a given management zone under two alternative conditions (current and best solution), ranging from 0 (no planning units in common) to 1 (all planning units in common). Jaccard =Best solutions ∩Current zones Best solutions ∪Currentzones (5) 3. Results 3.1. Targets and areas selected The areas selected by Marxan with Zones met representation targets Table 1 Target distribution for ES across management zones. Total targets accounted for 70% of the total supply in the 2005 scenario and were distributed according to their compatibility with agricultural practices. Zone Agriculture (ha) Carbon (Mg C ha −1 yr −1 ) Erosion control (normalized rates from t ha −1 yr −1 ) Recreation (number of PUs with high value) Transition 168386 25463 1079 361 Buffer 120275 63657 2699 1264 Core 48110 89121 3779 632 Total 336772 178243 7558 2529 Table 2 Target distribution for species (in number of occurrences) across management zones. Targets were set according to their total number of occurrences in the study area and distributed according to habitat preferences. Rare species are those with 200 or less occurrences across the study area. Common Species are those with more than 200 occurrences. Habitat Zone Rare Species Common Species Generalist & Open Habitat Transition 0 50 Buffer 0 50 Core All occurrences 100 Total All occurrences 200 Forest, Wetlands & Semi-open habitat Transition 0 0 Buffer 0 75 Core All occurrences 125 Total All occurrences 200 M.C. Iglesias et al. Journal of Environmental Management 322 (2022) 116045 8 for all species and ES under all management scenarios and RCPs. However, the number of planning units selected under each management zone differed among scenarios and RCPs (Fig. 3a). In the ‘FarmReturn’ and ‘AgroforestRe’ scenarios, and under both RCPs, the Transition area required around 200 km 2 less area than the 2005 scenario. ‘BAU’, especially under RCP 4.5, required about 200 km 2 more area than in 2005, and ‘Afforestation’ required a similar area under RCP 4.5, but less area under RCP 8.5. The Buffer zone remained fairly constant across future scenarios and RCPs, requiring an area only slightly higher in comparison to 2005. The same was observed for the Core area, but the increase in extent required in this case was higher. In all 3 zones, ‘BAU’ showed a degree of variability between replicates of the same scenario higher than the other scenarios. 3.2. Comparison with the current management plan The Core area in Marxan outputs was almost four times higher than in the actual planning of the BR (3872–3973 km 2 in Marxan vs. 1064 km 2 currently). In contrast, Buffer zones and Transition areas in Marxan covered less area compared to their actual extent (2777–2881 vs. 4203 km 2 for the Buffer zone and 2020–2517 vs. 6144 km 2 for the Transition area). Accordingly, the spatial overlap of the distribution of management zones in Marxan best solutions and the actual zoning of the BR was low across all management scenarios and RCPs, including the 2005 scenario, as evidenced by Jaccard index scores ranging between 0.08 and 0.26 (Fig. 3b). The Jaccard index was higher in the Transition and Buffer zones than in the Core area for all scenarios and RCPs. The Jaccard index for Core areas was highest under FarmReturn and RCP 8.5 and lowest under Afforestation and RCP 8.5. For Buffer zones, all landscape management scenarios showed Jaccard index scores lower than the 2005 scenario although variation was low. Jaccard values were higher for RCP 8.5 than for RCP 4.5. For the Transition area, only the ‘FarmReturn’ scenario (RCP 8.5) showed a distribution closer to the current zonation than the 2005 scenario while ‘BAU’ was the least similar under both RCPs. 3.3. Fire intensity within management zones Average potential fire intensity revealed important differences among management zones (Fig. 4). The Transition area generally presented lower fire intensity across management scenarios in comparison with Buffer and Core areas. Among management scenarios, ‘BAU’ showed higher fire intensities in the future compared to 2005 in all management zones but especially in the Core area. The other three scenarios showed fire intensities lower than in 2005, with ‘Afforestation’ showing the lowest fire intensities, except for the Core area under RCP 8.5 for which ‘AgroforestRe’ presented the lowest value. 3.4. Coverage of ecosystem services There were some scenarios where the representation of ES in areas aiming to secure different ES could lead to conflicts between objectives, such as the high representation of agricultural areas in the Core area under future scenarios (e.g., ‘Afforestation’ and ‘AgroforestRe’ and 2005 Fig. 3. a) Area (number of 1 km 2 planning units) allocated to management zones according to landscape management scenarios and RCPs. b) Similarity (represented by the Jaccard index) between the zonation of Marxan with Zones’ best solutions and the zonation currently implemented in the Meseta Ib´ erica Transboundary Biosphere Reserve. 2005 scenario is represented by a single line since we used one map only. Boxplots aggregate results of 10 runs made for each landscape management scenario. Lower and upper hinges of the boxplots correspond to the first and third quartiles (Q1 and Q3), while the vertical line inside the box represents the median. Lower whisker represents data at Q1 – 1.5* IQR and upper whisker represents data at Q3 +1.5 * IQR. Data beyond that range are called outliers and represented individually with points. M.C. Iglesias et al. Journal of Environmental Management 322 (2022) 116045 9 scenarios; Fig. 5). Regulating services had high representation in the Transition area under the Afforestation and BAU scenarios for Carbon sequestration and all scenarios for erosion control. Nature-based recreation was highly represented in the Transition and particularly in the Core area (Fig. 5). 3.5. Temporal turnover of species In the 2005 scenario, Marxan with Zones met representation targets for all species. However, under 2050 scenarios there were some missing targets due to the strong decline in distribution area or local extinction predicted for some species. The number of species whose presence is expected to decline below the representation target sought, ranged from 20 to 22 under RCP 4.5 and from 21 to 23 under RCP 8.5, which is around 10% of the total number of species in both cases (Table 3). In addition, there were 39 species (18.8%) under RCP 4.5 and 33 species (15.9%) under RCP 8.5 that were predicted to completely disappear from the study area by 2050 (Table 3). Contrastingly, there were species that were not initially in the area but are expected to be present in future scenarios: 10 species under RCP 8.5 and 11 species under RCP 4.5 representing (5% of total species) that always met their representation targets (Table 3). All other species met their targets both in 2005 and 2050 (“Persist” species in Table 3). 4. Discussion This study illustrates how a flexible approach based on splitting representation targets between different management zones could contribute to harmonising conservation with other management objectives, one of the major challenges of the Meseta Ib´ erica Biosphere Reserve as well as other conservation areas. This approach could also contribute to facilitating zoning implementation and securing larger targets (Lanzas et al., 2019). However, the implementation of this holistic approach would need careful coordination across all local stakeholders involved in land management to avoid conflicts between objectives and to identify best management practices (Abarca et al., 2022). Our simulations indicate that to conciliate biodiversity conservation with the sustainable supply of ES in the upcoming decades, changes in the distribution and extent of the zones of the BR would be required. To improve the role of Core areas, the current extent of this zone would need to be expanded threefold by 2050, with different spatial prioritisation depending on the land-use policy to be implemented in the future, as suggested for other protected areas of Mediterranean climate (Martinez-Harms et al., 2021; Lanzas et al., 2021; Regos et al., 2018). Core areas in future scenarios overlap, to a great extent, with existing Transition areas indicating that the required expansion of Core areas should Fig. 4. Mean potential fire intensity class per zone in each landscape management scenario and RCP. 2005 scenario is represented by a single line since we used one unique map. For landscape management scenarios, the 10 runs for each scenario are aggregated in boxplots. Lower and upper hinges of the boxplots correspond to the first and third quartiles (Q1 and Q3), while the vertical line inside the box represents the median. Lower whisker represents data at Q1 – 1.5* IQR and upper whisker represents data at Q3 + 1.5 * IQR. Data beyond that range are outliers and represented individually with points. M.C. Iglesias et al.