Received: 24 March 2021 | Accepted: 21 June 2021 DOI: 10.1111/1365-2745.13752 RECONCILING RESILIENCE ACROSS ECOLOGICAL SYSTEMS, SPECIES AND SUBDISCIPLINES Research Article Forest resilience to global warming is strongly modulated by local-­scale topographic, microclimatic and biotic conditions Jofre Carnicer1,2,3 | Maria Vives-­Ingla1,2 | Laura Blanquer1,2 | 1 1 1,2 Xavier Méndez-­Camps | Carme Rosell | Santi Sabaté | Emília Gutiérrez1 | Teresa Sauras1 | Josep Peñuelas2 | Adrià Barbeta1 1 BEECA-­UB, Department of Evolutionary Biology, Ecology and Environmental Abstract Sciences, University of Barcelona, 1. Resilience of endangered rear edge populations of cold-­adapted forests in the Barcelona, Spain 2 Mediterranean basin is increasingly altered by extreme heatwave and drought CREAF, E08193 Bellaterra (Cerdanyola del Vallès), Catalonia, Spain pressures. It remains unknown, however, whether microclimatic variation in these 3 IRBIO, Biodiversity Research Institute, isolated forests could ultimately result in large intra-­population variability in the University of Barcelona, Barcelona, Spain demographic responses, allowing the coexistence of contrasting declining and re- Correspondence silient trends across small topographic gradients. Multiple key drivers promoting Jofre Carnicer Email: jofre.carnicer@ub.edu spatial variability in the resilience of rear edge forests remain largely unassessed, including amplified and buffered thermal exposure induced by heatwaves along Funding information Spanish Government, Grant/Award topographic gradients, and increased herbivory pressure on tree saplings in de- Number: CGL2016-­78093-­R , PID2020-­ faunated areas lacking efficient apex predators. Here we analysed whether indica- 117636GB-­C21 tors of forest resilience to global warming are strongly modulated by local-­scale Handling Editor: Lorena Gomez Aparicio topographic, microclimatic and biotic conditions. 2. We studied a protected rear edge forest of sessile oak Quercus petraea, applying a suite of 20 indicators of resilience of tree secondary growth, including multidec- adal and short-­term indices. We also analysed sapling recruitment success, recruit/ adult ratios and sapling thermal exposure across topographic gradients. We found large within population variation in secondary growth resilience, in recruitment success and in thermal exposure of tree saplings to heatwaves, and this variability was spatially structured along small-­scale topographical gradients. Multidecadal resilience indices and curves provide useful descriptors of forest vulnerability to climate warming, complementing assessments based in the analysis of short-­term resilience indicators. Species-­specific associations of trees with microclimatic vari- ability are reported. 3. Biotic factors are key in determining long-­term resilience in climatically stressed rear edge forests, with strong limitation of sapling recruitment by increased roe Jofre Carnicer and Maria Vives-­Ingla contributed equally to this work and share first authorship. This is an open access article under the terms of the Creative Commons Attribution-­NonCommercial-­NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-­commercial and no modifications or adaptations are made. © 2021 The Authors. Journal of Ecology published by John Wiley & Sons Ltd on behalf of British Ecological Society Journal of Ecology. 2021;00:1–18.  wileyonlinelibrary.com/journal/jec | 1 2 | Journal of Ecology CARNICER et al. deer and wild boar herbivory. Our results also support non-­stationary effects of climate determining forest growth responses and resilience, showing increased negative effects of warming and drought over the last decades in declining stands. 4. Synthesis. Our findings do not support scenarios predicting spatially homogene- ous distributional shifts and limited resilience in rear edge populations, and are more supportive of scenarios including spatially heterogeneous responses, char- acterised with contrasting intra-­population trends of forest resilience. We con- clude that forest resilience responses to climate warming are strongly modulated by local-­scale microclimatic, topographic and biotic factors. Accurate predictions of forest responses to changes in climate would therefore largely benefit from the integration of local-­scale abiotic and biotic factors. KEYWORDS defaunation, drought, forest resilience, global warming, heatwave, land-­atmospheric feedbacks, mammal herbivory, microclimate 1 | I NTRO D U C TI O N Imbalanced networks of mammal herbivores stand as a potentially important, but currently unassessed, driver of reduced demographic Rear edge tree populations of cold-­adapted temperate species resilience in forest rear edge areas in the Mediterranean basin (Cairns are especially vulnerable to global change in the Mediterranean & Moen, 2004; Gómez & Hódar, 2008; Munier et al., 2010). Human-­ basin (Hampe & Petit, 2005; Pulido et al., 2019; Vilà-­C abrera & induced global defaunation processes can alter local interaction net- Jump, 2019). In Mediterranean rear edge forests, one of the most works and introduce drastic changes in the processes determining threatening components of global change is the increased frequency, forest demographic resilience (Estes et al., 2011; Young et al., 2016). intensity and duration of droughts and summer heatwaves (Carnicer, Studies in multiple land ecosystems indicate that the defaunation of Domingo-­Marimon, et al., 2019; Jump & Peñuelas, 2006; Spinoni top predators can facilitate an overabundance of large herbivorous et al., 2019; Turco et al., 2017; Vicente-­Serrano et al., 2014). Despite mammals and drastically change forest dynamics (Estes et al., 2011; increasing climatic impacts, recent studies report limited demo- Gómez & Hódar, 2008; McLaren & Peterson, 1994; Nogués-­Bravo graphic declines and retraction in the distribution of many temperate et al., 2016; Ripple & Beschta, 2007; Roemer et al., 2009; Smit trees in the Mediterranean basin, suggesting sustained demographic et al., 2015). For example, in the eastern Iberian Peninsula the his- resilience to current climatic impacts (Batllori et al., 2020; Dittmar torical defaunation of apex predators (wolves) has facilitated over et al., 2003; Granda et al., 2018; Hacket-­Pain & Friend, 2017; Tegel the last decades the expansion of consumers of tree saplings and et al., 2014; Vilà-­C abrera et al., 2019). However, the mechanisms en- tree roots such as the Roe deer Capreolus capreolus and the wild hancing forest resilience in spite of increased warming impacts in boar Sus scrofa (Giménez-­Anaya et al., 2020), increasing pressure rear edge areas still remain poorly understood for most dominant on tree saplings at the regional scale. Tree sapling recruitment is a tree species (Canadell & Jackson, 2021; Granda et al., 2014; Lloret & key process determining demographic resilience in rear edge areas Batllori, 2021; Lloret et al., 2012; Muffler et al., 2020; Vilà-­C abrera & of Mediterranean forests (Carnicer et al., 2014; Granda et al., 2014). Jump, 2019). Altitudinal and topographic gradients in rear edge areas This process could be largely impaired by increased herbivory in can determine locally heterogeneous demographic responses to in- areas affected by long-­term processes of defaunation of top pred- creased warming (Granda et al., 2018; Hacket-­Pain & Friend, 2017). ators and reduced human hunting pressure on large mammal her- However, the modulation by local-­scale topographic gradients of bivores (Côté et al., 2004; Mayer et al., 2000; Partl et al., 2002). In the impacts of extreme summer heatwaves in relict forests remains Mediterranean rear edge areas large population expansions of wild poorly quantified. Similarly, studies assessing the negative synergis- boar and roe deer are currently co-­occurring, potentially affecting tic effects of the multiple drivers of global change in the resilience of key rear edge areas (Barrios-­Garcia & Ballari, 2012; Giménez-­Anaya rear edge forests are also warranted. For example, we lack detailed et al., 2020; Massei & Genov, 2004; Saez-­Royuela & Telleria, 1986). studies jointly considering in rear edge areas the coupled effects of Management efforts to control increasing populations of wild boar increased drought, extreme summer heatwave events and altered and roe deer often obtain limited results, due to multiple logistic networks of mammal herbivores negatively affecting tree sapling limitations and social barriers impeding an optimal human hunting recruitment (Cairns & Moen, 2004; Castro et al., 1999, 2004; Den pressure (Giménez-­Anaya et al., 2020; González-­Crespo et al., 2018; Ouden et al., 2005; Munier et al., 2010). Massei et al., 2015; Vicente et al., 2019). CARNICER et al. Journal of Ecology | 3 Beyond the effects of altered networks of mammal herbivores, inland from the Mediterranean Sea. This protected area (13,694 ha) other local-­scale factors, such as local gradients of microclimatic and comprises mid-­altitude mountain ranges (300–­1104 m) covered by soil conditions can importantly determine the vulnerability of rear mixed Mediterranean oak and pine stands (Quercus ilex, Quercus hu- edge forests to global change (Graae et al., 2018; Granda et al., 2014; milis, Quercus petraea, Pinus halepensis, Pinus nigra and Pinus sylves- Herrera & Bazaga, 2008; Lenoir et al., 2017; Vilà-­C abrera et al., 2019; tris). Forests are dominated by mature successional stages (with ages Zellweger et al., 2020). Topographical features (e.g. the steepness in dominant trees >100 years), with abundant and well-­developed and aspect of slopes, elevation) and canopy cover can largely mod- holm and sessile oaks and limited pine regeneration dynamics ify local heat exchanges at fine scales, altering climatic exposure of (Carnicer et al., 2014). The current forest structure is the result of a saplings and adult trees and potentially increasing forest resilience multidecadal period of rural abandonment since the 1950s and the to climatic extremes (Bramer et al., 2018; Carnicer, Stefanescu, lack of active forest management in most of the protected areas of et al., 2019; Carnicer et al., 2017; Körner & Paulsen, 2004; Scherrer the park. Quercus petraea stands are restricted to isolated mesic sites & Körner, 2010). Previous studies highlight that heterogeneity of mi- in north-­facing hillsides, whereas the most xeric sites are covered croclimatic conditions at these scales can be especially relevant for by Mediterranean mixed stands of holm oak (Q. ilex), Aleppo pine seedling establishment and tree recruitment (Dobrowski et al., 2015; (P. halepensis) and xeric shrublands. In the study area, sessile oaks Lloret et al., 2005; Príncipe et al., 2019; von Arx et al., 2013). stands are distributed in a relatively narrow topographic band To effectively assess demographic resilience in forest rear edge (750–­890 m), dominating the northern side of the mountain range areas under pressure, detailed integrated measurements of microcli- and completely surrounded by holm oak forests. Sessile oak stands matic buffering effects, exposure of saplings to summer heatwaves are restricted to humid areas but sparse individuals also colonise and and mammal sapling predation are needed. To address this knowl- develop in the elevated areas of southern-­facing slopes (Figure 1B). edge gap, here we aim to evaluate the demographic resilience of a While Q. ilex was fairly abundant and present in all topographic rear edge population of sessile oak (Quercus petraea Liebl.), assess- areas sampled, sessile oak saplings and adults were absent or ing tree growth, recruitment rates and herbivore damage patterns rarely observed in the warmer southern-­oriented bands A5 and A6 along a topographic gradient in a protected Mediterranean moun- (Figure 1B). In contrast to these Mediterranean isolated rear edge tain range. For this, we will address the following specific objectives: areas, sessile oak is a very abundant tree at higher latitudes of the (a) to measure microclimatic exposure in oak tree saplings during an European continent, growing preferentially in humid areas in both extreme summer heatwave event attributed to climate change (van calcareous and siliceous soils (Figure 1A). Pine species in the sam- Oldenborgh et al., 2019), assessing the operation of local buffering pled area were represented by a few isolated and very scarce adult effects of the topographic gradient on the observed thermal expo- individuals of P. halepensis (located in the south-­facing slope areas, sure, (b) to describe the variability of tree sapling recruitment and i.e. the topographic bands A4–­A6 in Figure 1B). A few isolated indi- herbivore damage along the topographic gradient, (c) to analyse the viduals of P. sylvestris were also detected in elevated northern-­facing contrasting patterns of tree secondary growth and resilience be- areas (i.e. corresponding to the topographic band A3 in Figure 1B). tween climatically exposed and buffered areas of the topographic Due to this marginal presence and representation in the studied for- gradient. est stands, trends for pine species were subsequently not consid- ered in the analyses of forest resilience. 2 | M ATE R I A L S A N D M E TH O DS 2.1.1 | Extreme summer heatwave 2.1 | Study area event and quantification of microclimatic exposure in tree saplings More than 25 Quercus species inhabit the Mediterranean basin (Castagneri et al., 2020), including relict populations of deciduous In late June 2019 south-­western and central Europe experienced oak species such as the sessile oak (Q. petraea). This species shows a the hottest summer heatwave event recorded in the available in- higher susceptibility to drought in comparative ecophysiological stud- strumental records (Sousa et al., 2019; Zhao et al., 2020). The ies (Arend et al., 2011; Bréda & Granier, 1996; Corcuera et al., 2002; heatwave was related to a persistent atmospheric wave train propa- Shestakova et al., 2014). This study analyses demographic resilience gating eastward from the mid-­high latitudes of the North Atlantic in an isolated sessile oak population located at the southernmost to Eurasia (Zhao et al., 2020) and has been linked to climate change limit of the species distribution in the Iberian Peninsula (Figure 1A). in the available attribution studies (van Oldenborgh et al., 2019). This forest stand is remarkable for its large extension (4 ha), maturity Record-­breaking temperatures for the month of June were re- and conservation state, given that along rear edge Mediterranean corded at many locations in Europe, including France, Switzerland, areas well-­developed stands of Q. petraea are rather uncommon due Austria, Germany, the Czech Republic, Italy and Spain. For example, to intensive logging over the last century (Bou & Vilar, 2019). The France experienced temperatures in excess of 45℃ for the first time studied population is located within the Natural Park of Sant Llorenç in recorded history (van Oldenborgh et al., 2019). Along the stud- del Munt i l’Obac, in Catalonia (NE Iberian Peninsula), at 40 km ied topographic gradient, we measured the modulation of thermal 4 | Journal of Ecology CARNICER et al. (A) (B) (C) (D) (E) (F) F I G U R E 1 The study zone is a rear edge forest of Quercus petraea, which has been recently impacted by heatwave dynamics, and increasing mammal herbivore pressure dynamics. (A) The map illustrates the location of the isolated rear edge population of Q. petraea (star symbol), in relation to the continuous continental distribution of this species in western Europe. (B) A diagram of the sampling units defined along the topographic gradient in the rear edge population (topographic bands A1–­A6). Three topographic bands were distributed between 750 and 850 m [A1–­A6: 750–­783 m (elevational range: 33 m); A2–­A5: 784–­817 m; A3–­A4: 818–­850 m]. Black dots indicate the location of temperature–­humidity sensors allowing the measurement of microclimatic variability and tree sapling thermal exposure during an extreme summer heatwave (June 2019). As illustrated in the diagram, in the studied rear edge forest the drought-­sensitive Q. petraea is restricted to northern slopes and high elevation areas (i.e. ridge areas, 850 m a.s.l.), being progressively replaced by Quercus ilex stands in south facing slopes (topographic bands A5 and A6 in panel B). (C) Observed temporal dynamics of midday understorey temperatures at 25 cm height during the studied extreme summer heatwave period (coloured rectangle). (D) Observed temporal dynamics of air humidity (25 cm height) in the understorey during the extreme summer heatwave event. (E) Temporal dynamics of wild boar and roe deer densities in the studied forest rear edge area, showing a sustained increase over the last decade. Roe deer densities are expressed in percentage of detection in field censuses (Rosell et al., 2019). The orange line depicts the hunting effort in the study area (hunted wild boar individuals/km2). (F) A diagram summarising the multidecadal and short-­term resilience indices applied in the study of secondary growth [measured as Basal Area Increment (BAI) in declining (red line) and non-­declining areas (black line) of the rear edge forest] exposure in continuously monitored oak saplings in the peak of this bands (A1–­A6 in Figure 1B, 750–­850 m). These topographic bands unprecedented heatwave event. To assess how topographic vari- comprised the north and the southern slope of the studied rear edge ability, aspect and slope modified tree sapling thermal exposure and forest (Figure 1B). To continuously assess the spatial variation in cli- humidity conditions, the study zone was divided into six topographic matic exposure, six temperature–­humidity sensors (EL-­USB-­2, Lascar CARNICER et al. Journal of Ecology | 5 Electronics) were distributed in early spring across the mountain the late June extreme heatwave event, affecting soil energy balance range (Figure 1B, black dots), covering the whole altitudinal range of during the event (thermal amplification processes, see Table S1 for the sessile oak stand. The automatic sensors were located at 25 cm a precise definition). height above the soil surface using metal stakes, and were protected from direct solar radiation by a plastic envelope sustained by a wire-­ mesh cylinder (installed 5 cm above the sensor and thus precluding 2.1.2 | Altitudinal patterns of herbivore damage, the direct incidence of solar radiation), recording microclimate data sapling recruitment and forest structure with hourly frequency. To assess thermal exposure outside the for- est stand, an additional LASCAR temperature–­humidity sensor was To assess the abundance of recruits and its spatial variability we es- installed in an open bare ground area in the southern slope of the tablished a total of 36 linear transects covering the whole topographic range (800 m, in the middle of the topographic band A5, in a for- gradient (establishing six transect replicates per topographic band, per- est gap area). This additional LASCAR sensor was also installed at pendicular to the elevational gradient, see Tables S2 and S3; Figure S1). 25 cm height above the soil surface and protected from direct solar Tree sapling recruits and adult individuals were categorised in four size radiation. In addition, using portable thermocouples and thermo-­ categories following the Spanish Forest National Inventory classifi- hygrometer devices, we measured with monthly resolution in all cation (i.e. R1–­R4 for saplings and SC1–­SC4 for adult trees; Table S2; topographic bands (A1–­A6) spring and early summer midday soil hu- Carnicer et al., 2014; Coll et al., 2013). Changes in tree density (individ- midity (11:00 a.m. to 16:00 p.m.), soil temperature (1–­15 cm depth), uals/ha) across the topographic bands A1–­A6 were assessed for each soil surface temperature, tree sapling leaf temperatures and air tem- diametric class for holm and sessile oak, applying linear, loess and poly- perature at 2.5, 5, 10, 25, 50 and 100 cm height. Six measurement nomial models. Significant differences between topographic bands points were established in each topographic band (A1–­A6) and in were evaluated applying Tukey–­Kramer tests. To assess recruitment the control forest gap area (A5-­open area without canopy) and were success in the forest for Q. petraea and Q. ilex we calculated Recruit/ visited in April, May and June. We used temperatures measured at Adult indices (Table 1A) in all the topographic bands of the gradient. 1 m height as a reference for calculating positive and negative ther- Two independent indicators of mammal herbivore damage on sap- mal offset values at different height levels (see Table S1, De Frenne lings were evaluated [terminal twig damage by roe deer (TTD), and et al., 2021). Canopy cover (percentage) was visually estimated in all rooting damage by wild boar (RD), Table S4]. The patterns of herbivore 2 these sites in a zenithal area of 1 m above the focal point of ther- damage along the elevational gradient were evaluated (topographic mal measurement (see Carnicer, Stefanescu, et al., 2019; Vives-­Ingla bands A1–­A6). In each topographic band, six parallel transects were et al., 2020 for methodological details). established quantifying terminal twig damage (TTD) in all recorded In the six sampling points established in each topographic band saplings for R1 and R2 size classes. For the assessment of soil sur- monthly soil humidity measurements were performed using DELTA-­T face affected by wild boar rooting (RD), three plots (1 m2) were ran- SM150 (Delta-­T Devices Ltd), ambient humidity measures were domly measured in each transect (i.e. 18 plots per each topographic taken with Thermo-­hygrometer DC108 (Bonajay Technology Ltd). A band). Table S3 summarises the species-­specific transect length and thermocouple Hanna HI935005N type K (Hanna Instruments Ltd) width applied for each tree species and size category (Quercus petraea, was also used for the measurement of leaf, soil surface and air tem- Quercus ilex). Terminal twig damage (TTD) in tree saplings has been peratures. Soil underground temperature was measured employing previously applied as a reliable indicator of herbivore activity in roe Hanna HI98509 (Hanna Instruments Ltd). In each topographic band deer populations (Partl et al., 2002). In the case of wild boar popula- (A1–­A6), we randomly selected and marked three saplings of both tions in the study area, acorns and the underground parts of plants holm and sessile oak for continuous monitoring of temperature and (roots, bulbs and rhizomes) constitute an important fraction of the humidity conditions with monthly frequency using thermocouples. diet (Giménez-­Anaya et al., 2020). In this regard, indicators of rooting Saplings were <130 cm high, corresponding to the first and second activity (RD) in the soil provide a good proxy of the spatial variability sapling size categories of the Spanish National Forest Inventory (R1 of damage by wild boar populations (Bueno et al., 2009, Table S4). and R2, Table S1, Coll et al., 2013). All saplings selected were iden- tified with plastic labels and georeferenced. For all microclimatic thermal measurements, we calculated thermal warming and cooling 2.1.3 | Resilience indices, secondary tree growth and effects relative to reference measurements (i.e. local air tempera- relationships with climatic conditions tures measured at 1 m height without direct radiation incidence, see Carnicer et al., 2017; Carnicer, Stefanescu, et al., 2019). The heat- We applied a suite of forest resilience indices to contrast, along the wave event of late June of 2019 was preceded by an anomalous dry topographic gradient, the observed responses of secondary growth period. Standardised Precipitation Evapotranspiration drought indi- of Q. petraea in climatically exposed and climatically buffered areas ces (SPEI; extracted from SPEIbase v2.6; https://spei.csic.es/datab​ (Carnicer, Domingo-­Marimon, et al., 2019; Lloret et al., 2011). The in- ase.html) indicated a sustained preceding drought period, reporting dices were applied to dendrochronological samples of the forest (see −1.8 SD units in June SPEI 6 relative to the long-­term mean (1950–­ below, Table S5; Figure S2). As specified in Table 1B, we applied both 2010). The preceding period of soil drought possibly interacted with short-­term resilience indices (Lloret et al., 2011) and multidecadal 6 | Journal of Ecology CARNICER et al. TA B L E 1 A summary of the recruitment success indices (A) and the resilience, resistance and recovery indices (B) applied in the study of the trends of secondary growth in declining and non-­declining areas of the rear edge forest. A detailed mathematical description for each index is provided in Table S6 Index Type Description Reference (A) Recruitment success indices R sp/AT ratio (Recruit/Adult ratio) S Ratio of the density of recruits of a focal tree species (R sp) relative This paper to the density of multiple tree species (AT ) of the dominant trees inhabiting the canopy R sp/A sp ratio (Specific Recruit/ S Ratio of the density of recruits of a focal species (R sp) relative to This paper Adult ratio) the density of a single tree species (A sp) of the dominant trees inhabiting the canopy (B) Secondary growth resilience indices RS (Resilience) S Capacity to reach pre-­disturbance performance levels. Ratio Lloret et al. (2011) between the performance after and before disturbance MRS (Multidecadal resilience) M Resilience index relative to a reference period before the onset This paper of forest decline. The same reference period is used for multiple drought events MRCFW (Multidecadal M Multidecadal change in resistance and resilience. Beginning at the This paper. Carnicer, Domingo-­ resistance–­resilience curve drought year D, a sequence is calculated such that, for each year Marimon, et al. (2019) (fixed window)) y, the growth of a period w from y is divided by the growth of a reference period MRCIW (Multidecadal M Multidecadal change in resistance and resilience. Beginning at the This paper. Carnicer, Domingo-­ resistance–­resilience curve drought year D, a sequence is calculated such that, for each year Marimon, et al. (2019) (increasing window)) y, the growth from D to y is divided by the performance during a reference period RC (Recovery) S Ability to recover relative to the damage experienced during Lloret et al. (2011) disturbance. Ratio between performance after and during disturbance MRCFW (Multidecadal recovery M Recovery relative to the growth observed during a severe drought This paper. Carnicer, Domingo-­ curve (fixed window)) event. For each year y after the drought, the difference between Marimon, et al. (2019) the growth of a period w beginning at y and of the drought is divided by this last term MRCIW (Multidecadal recovery M Recovery relative to the growth observed during a severe drought This paper. Carnicer, Domingo-­ curve (increasing window)) event. For each year y after the drought, the difference between Marimon, et al. (2019) the growth from the first year to y and of the drought is divided by this last term MD (Multidecadal decline) M Observed actual growth relative to a reference period before the This paper decline RT (Resistance) S Reversal of the reduction in ecological performance during Lloret et al. (2011) disturbance. Ratio between the performance during and before the disturbance MRT (Multidecadal resistance) M Resistance index relative to a reference period before the onset This paper of forest decline GRT (Growth recovery time) S Duration (in years) starting in the drought year until reaching the Thurm et al. (2016) pre-­drought growth levels again. Recovery period MGRT (Multidecadal growth M Duration (years) starting in the drought year until reaching the Modified from Thurm et al. (2016) recovery time) levels of growth before the onset of forest decline TGR (Total growth reduction) S Cumulative growth reduction in the drought year plus all years in Thurm et al. (2016) the recovery period MTGR (Multidecadal total M Cumulative growth reduction in the drought year plus all years in Modified from Thurm et al. (2016) growth reduction) the multidecadal recovery period AGR (Average growth reduction) S Ratio of total growth reduction divided by the recovery period. Schwarz et al. (2020) Annual mean of growth reduction during the recovery period MAGR (Multidecadal average M Ratio of multidecadal total growth reduction divided by the Modified from Schwarz et al. (2020) growth reduction) length of the multidecadal recovery period (Continues) CARNICER et al. Journal of Ecology | 7 TA B L E 1 (Continued) Index Type Description Reference ARR (Average recovery rate) S Magnitude of growth reduction in the drought year (i.e. 1 − RT ) Modified from Schwarz et al. (2020) divided by the length of the recovery period. How much growth is averagely recovered each year of the recovery period MARR (Multidecadal average M Magnitude of growth reduction in the drought year (i.e. 1 − MRT ) Modified from Schwarz et al. (2020) recovery rate) divided by the length of the multidecadal recovery period PARR (Percental average S Average recovery rate relative to the growth reduction in the Modified from Schwarz et al. (2020) recovery rate) drought year. Percentage of the lost growth that is averagely recovered each year of the recovery period PMARR (Percental multidecadal M Multidecadal average recovery rate relative to the growth Modified from Schwarz et al. (2020) average recovery rate) reduction in the drought year. Percentage of the lost growth that is averagely recovered each year of the multidecadal recovery period resilience indices (Carnicer, Domingo-­Marimon, et al., 2019, obtained from the NERC Data Repository for Atmospheric Science Figure 1F). A total of 20 indices were calculated, including resistance, and Earth Observation (https://catal​ogue.ceda.ac.uk/). Standardised recovery and resilience indices (Table 1; Table S6). Note that while Precipitation Evapotranspiration drought indices (SPEI) were obtained most indices reported an averaged annual value of forest resilience, from the dataset SPEIbase v2.6 (https://spei.csic.es/datab​ase.html). some of them specifically provided an integrated value for a larger, The relationship of secondary growth (Basal Area Increment, BAI, interannual period of time (Table 1B). Table S6 provides a detailed dependent variable) with climatic variables was assessed using both mathematical description for each of the indices. An R code facilitat- GLM and GLMM approaches, analysing significant associations with ing the calculation of this integrated suite of indices is provided in the monthly temperatures and precipitation. In GLM analyses, we applied Supplementary file 1 (Table 1; Table S6). The analysis of microclimatic a model selection approach (AIC criterion) and selected the strong- time series (i.e. temperature–­humidity 25-­cm height LASCAR sensors est climatic predictors, using the functions glm and step from the stats data) indicated that the two lowland topographic bands of the forest package, and used as the response variable the group-­averaged BAI located in northern slopes (i.e. A1–­A2, valley area in Figure 1B) were (for valley and ridge areas). For GLMM analyses, we used the R pack- significantly buffered during extreme heatwave events. Over the age lme4, applying a Gaussian family with identity link. In this case, course of the heatwave event these lowland bands experienced sig- individual-­level BAI was the response variable and we added tree nificantly lower extreme summer temperatures and higher air humid- identity as a random factor. All statistical and dendrochronological ity values relative to higher topographic bands inhabited by Quercus analyses were done using R 4.0.3 (R Core Team, 2020). A period of petraea [R2 = 0.67, p < 0.0001 (lower midday June temperature); continued forest decline in secondary growth was detected (1980–­ R 2 = 0.69, p < 0.0001 (higher midday air humidity)]. This different ex- 2019), preceded by a previous period of more stable growth (1960–­ posure was assessed analysing LASCAR sensor data for June 2019. To 1979). We analysed non-­stationary effects of the climatic variables contrast forest growth trends in these climatically exposed and buff- on growth for these two contrasting time periods (e.g. Carnicer, ered topographic bands two complementary approaches were im- Domingo-­Marimon, et al., 2019). The analyses of resilience indices fo- plemented. First, for selected dominant trees of Q. petraea located in cused on the major periods of drought detected in each decade in the either valley areas (A1–­A2 bands) or ridge areas (A3–­A4) we obtained period of forest decline (1980–­2019), corresponding to the summers dendrochronological time series (see Table S5). We compared growth of 1982, 1999, 2006 and 2015 (Figure S3). In addition, we developed trends in these two areas exposed to significantly different thermal a detailed analysis of the relevant time-­scales for secondary growth regimes (Figure 1C,D). We applied the recommended optimal sam- of summer drought events, applying monthly SPEI indices of variable ple size effort indicated for Q. petraea in previous dendrochronologi- time length for all seasons (1–­24 months, summarised in Figure S4). cal studies (following Mérian & Lebourgeois, 2011), achieving robust Complementarily, the relationships of BAI with SPEI drought index EPS values (>0.96, Table S5). Two cores per adult tree were sampled at multiple monthly time-­scales were also assessed (1, 3, 6, 12 and using a 5 mm diameter increment borer, each corresponding to the 24 months, results not shown). eastern and western direction respectively. Samples were processed, dried, mounted and glued to wood supports. Cores were visually dated with a binocular microscope. Tree ring-­width series were cal- 3 | R E S U LT S culated for each forest area using the software WinDENDRO (Regent Instruments Inc. 2019) and dplR package in R (Bunn et al., 2020). 3.1 | Observed variability in microclimatic Cross-­dating was carried out using the software COFECHA (Holmes, conditions 1983). From the selected series, Basal Area Increment was calculated as a proxy of secondary tree growth. To analyse the climatic relation- During the heatwave event, we observed strong, significant differ- ships with tree growth trends, climatic data for the study site were ences in climatic exposure across microsites along the elevational 8 | Journal of Ecology CARNICER et al. June heatwave peak. 12 am June heatwave peak. 12 am Plant thermal exposure Sapling leaf temperature (º C) (A) (B) (C) Air relative humidity (%) 40 a ab bc c c bc c b ab a a ab 34 Air temperature (ºC) Q. ilex 35 33 Q. petraea 60 Northern Southern 30 aspect aspect 32 25 31 40 20 Northern Southern 30 Northern Southern aspect aspect aspect aspect 15 20 29 Elevation 750 850 750 750 850 750 750 850 750 (m a.s.l.) A1 A2 A3 A4 A5 A6 A1 A2 A3 A4 A5 A6 A1 A2 A3 A4 A5 A6 Topographic band Topographic band Topographic band F I G U R E 2 Climatic exposure of tree saplings during the extreme heatwave event of June 2019. (A) Midday temperature and (B) Midday relative humidity recorded by LASCAR EL-­USB-­2 temperature–­humidity sensors in the understorey. In the panels A and B different lower case letters indicate significantly different exposure between topographic bands (A1–­A6) in Tukey–­Kramer tests (p < 0.0001). (C) Observed thermal exposure in leaf saplings (midday values, 12–­14:00 hr) during the peak of the heatwave event (27 June). A spline fit is illustrated describing the observed variation in leaf temperature of the saplings of Quercus petraea (grey line) and Quercus ilex (black line) gradient, both in thermal and air humidity values (Figure 2A,B) and We observed a significantly negative recruitment trend in in leaf sapling temperature (Figure 2C). Figure S5 reports thermal Q. petraea, showing a strong decrease in the abundance of recruits exposure values observed at midnight. The observed temporal dy- between consecutive recruit size classes (Figures 4C and 5A). In namics of temperature and humidity conditions in the topographic contrast, Q. ilex showed sustained recruitment trends and non-­ bands A1–­A6 are described in Figure 1C,D. During the extreme significant differences were observed between size classes R1 and heatwave event, the thermal exposure of oak saplings below the R2 (Figure 5B). Recruitment decrease in Q. petraea saplings of size understorey was largely buffered when compared to control sites R2 was coherent with the observed high rates of predation of de- in open forested gaps completely exposed to direct sun radiation, ciduous oak saplings by roe deer (terminal twig damage, Figure 5C; which were not successfully colonised by oak saplings (Figure 3). Figure S8). A sustained increase in herbivore pressure due to signif- Significant effects of local thermal amplification and thermal buffer- icant increases in wild boar and roe deer abundances has been ob- ing processes were detected, determining sapling exposure during served in the area of study (Figure S9, Rosell et al., 2019). The rooting the heatwave event (Figure 3). No significant differences of canopy damage by wild boar was spatially concentrated and dominant factor cover were detected across the topographic bands dominated by in river bed and trail areas (Figure S10). Q. petraea in the period of the heatwave event (A1-­A4; Tukey–­Kramer test, F = 1.28; p = 0.28). Tree cover was elevated in most areas of the forest (>65%–­85%), and consequently understorey temperatures 3.3 | Contrasting patterns of tree during the heatwave were non-­strongly linked to minor variations growth and resilience between climatically in tree canopy cover percentage (Tables S7 and S8). However, a high exposed and buffered areas percentage of Q. ilex treelets in the understorey was associated with local thermal cooling patterns (Table S8; Figure S6). The time series of Basal Area Increment (BAI) of Q. petraea along the topographic gradient showed contrasting declining and non-­ declining trends in different areas of the rear edge (Figure 6A). Since 3.2 | Patterns of recruitment and herbivore damage the 1980s a divergent, heterogeneous spatial trend in secondary in Quercus ilex and Quercus petraea growth was observed, characterised by significant declining growth rates in the more thermally exposed mountain ridge areas [A3–­A4 The topographic gradient was associated with contrasting patterns bands] and positive growth trends in the valley areas [A1–­A 2 bands]. in the abundance and distribution of tree species, with adult indi- These divergent long-­term trends in secondary growth were more viduals of Q. petraea clearly dominating in northern topographic efficiently captured by multidecadal resilience indices (Figure 6B,C; bands characterised with milder thermal regimes. In contrast adult Table S9; Figures S11–­S14). As illustrated in Figure 6B, ridge areas individuals of Q. ilex were significantly more abundant in southern showed a significant reduction in multidecadal resilience values in topographic bands (Figure 4A). While tree saplings of Q. ilex were the last decades (i.e. 1999, 2006, 2015 drought events) when com- distributed across the whole altitudinal gradient with elevated den- pared to the observed impacts on resilience in the 1980s. Moreover, sities, saplings of Q. petraea were abundant only in the topographic the multidecadal resilience index was significantly lower in ridge bands A1–­A4 and were rarely recorded in the southern aspect areas relative to valley areas (asterisks in Figure 6B). Resilience curves (Figure 4B; Figure S7). (MRCFW ) also indicated a sustained decrease in resilience in the last CARNICER et al. Journal of Ecology | 9 (A) (B) (C) (D) F I G U R E 3 Contrasting midday thermal exposure in areas buffered by the forest canopy and a control area without canopy cover (forest open gap, topographic band A5). Thermal measurements were made in the southern aspect (A5) during the extreme heatwave event of June 2019 using portable thermocouples (panels A and B). A description of observed soil humidity dynamics across the whole topographic gradient is provided in panel C (April–­June period). (A) Observed thermal variability in A5 in two spatially adjacent areas [open forest gap (brown)/forested area (blue)] in relation to the distance to the ground. (B) Observed trends of thermal offset (see Table S1 for a definition), describing amplification and buffering effects relative to air temperatures measured at 1 m height above the ground in the same focal point. (C) Observed changes in soil humidity in the drought period preceding the heatwave event (April–­June period, for all topographic bands). The INSET graphic highlights observed significant differences (One-­way ANOVA, R 2 = 0.23, p < 0.01) in soil humidity of June between the control forest gap and forest understorey of the band A5 (i.e. southern aspect, 784–­817 m). Soil moisture in A5 was extremely low in June (showing a mean value of 1.18 ± 0.35% in gap areas, and 2.88 ± 0.52% in adjacent understorey areas covered by Quercus ilex trees). (D) A visual scheme summarising the dependence of sapling thermal exposure on the preceding soil drought period (panel C) and the cooling effects produced by the canopy structure (panels A and B). Blue circles represent canopy cooling effects reported in panels A and B. The orange circle represents the extreme thermal exposure observed in areas without canopy cover and significantly lower soil humidity (reported in panels A–­C). Red arrows schematically represent the possible interactions between the documented preceding drought period (panel C) and the late June heatwave event, determining thermal exposure on saplings due to vegetation and other land-­atmospheric feedbacks (reviewed in Table S1) decades in climatically exposed areas (Ridge, Figure 6C; Figure S13). summer drought (June rains and temperatures) were dominant and Short-­term resilience indices were not fully capturing the observed explained a larger fraction of growth variability (Table 2; Figure S15). declining trends (Figure S11; Table S9), but provided complementary The observed negative effects of summer drought and positive ef- information on short-­term recovery capacity. The observed correla- fects of winter and early spring temperatures on secondary growth tionships between all the short-­term and multidecadal indices listed are fully consistent with previous detailed ecophysiological studies in Table 1 are summarised in Figure 6D and Figure S14. The correla- of the determinants of secondary growth in this species (Bréda & tionships indicate that multidecadal and short-­term indices provide Granier, 1996, see the discussion for further details). qualitatively different and complementary information on forest re- silience (green arrow in Figure 6D). Variable non-­stationary effects of climatic variables were de- 4 | D I S CU S S I O N tected between the contrasting periods of forest growth [stable pe- riod (1960–­1979); declining period (1980–­2018)] in GLM and GLMM We have applied a diverse suite of 20 forest resilience indices to as- analyses (Table 2, Tables S10–­S12). In the first period positive ef- sess intra-­population trends in an endangered rear edge population fects on secondary growth of warm winter and spring temperatures of Q. petraea. As suggested by recent works (Capdevila et al., 2020; were the dominant factors. In contrast, during the warmer period Carnicer, Domingo-­Marimon, et al., 2019; Ingrisch & Bahn, 2018; (1980–­2018), the negative effects of climatic variables linked to Lloret et al., 2011; Schwarz et al., 2020; Thurm et al., 2016), our results 10 | Journal of Ecology CARNICER et al. (A) Abundance of adult trees (ind./ha) (B) (C) 3 Northern aspect Southern aspect Recruit/Adult ratio (R1) Recruit/Adult ratio (R2) Q. petraea Q. ilex 1,000 R2 OLS = 0.24* R2 OLS = 0.46**** Q. ilex 2 AICGLM = 77.84 2 AICGLM = 82.60 Q. petraea 500 R2 OLS = 0.60**** AICGLM = 5587.7 1 1 R OLS = 0.61**** 2 Q. ilex AICGLM = 277.86 R2 OLS = 0.17* R2 OLS = 0.20* AICGLM = 20.86 AICGLM = 75.45 0 Q. petraea Elevation 750 850 750 750 850 750 750 850 750 (m a.s.l.) A1 A2 A3 A4 A5 A6 A1 A2 A3 A4 A5 A6 A1 A2 A3 A4 A5 A6 Species R1 Consecutive recruit size classes R2 Quercus ilex Quercus petraea Topographic gradient F I G U R E 4 Effects of topographic heterogeneity on adult tree and sapling abundance of Quercus trees in the studied rear edge forest. (A) Adult tree abundance for Quercus petraea and Quercus ilex across the topographic gradient (including the larger size categories SC2, SC3 and SC4). Two-­and three-­order polynomial fits were applied, fitting generalised linear models (GLM with a Poisson distribution and log link) and ordinary least squares (OLS) models. GLM fits are illustrated, and the associated Akaike information criterion (AIC) is indicated. Complementarily, the variance (R 2) and significance (asterisks) reported by OLS polynomial fits are also indicated (*p < 0.05, **p < 0.01, ****p < 0.0001). (B) Observed changes in recruit/adult ratio (R sp/AT ratio, Table 1A) along the topographic gradient for the recruit size class R1. (C) Observed changes in recruit/adult ratio for the recruit size class R2. The green arrow highlights the sharp reduction on recruitment success in Q. petraea between consecutive size classes (R1 and R2) Sapling preference (A) Quercus petraea (B) Quercus ilex R = 0.15, p = 0.001 2 R = 0.02, p = 0.24 2 (C) 100 Sapling abundance (ind./ha) 3,000 Terminal twig damage (%) (herbivory by Capreolus) 4,000 75 2,000 50 2,000 Recruitment failure 1,000 25 R2 = 0.22, 0 0 p = 0.002 0 R1 R2 R1 R2 Quercus ilex Quercus petraea Tree species Consecutive sapling size classes F I G U R E 5 Observed trends in sapling recruitment success, and differences in herbivore damage by roe deer Capreolus capreolus. (A) Observed recruitment failure in R2 saplings of Quercus petraea across two consecutive size classes (R1–­R 2). (B) Sustained recruitment success in Quercus ilex. (C) Observed differences in terminal twig damage associated with roe deer herbivory between Q. petraea and Q. ilex. A one-­way analysis of variance (ANOVA) is provided in the three panels CARNICER et al. Journal of Ecology | 11 (A) Forest area (B) 16 Valley Resilient **** *** **** Multidecadal resilience Basal area increment 1.5 Forest area 12 index (MRS) (cm2 /year) Valley 1.0 Ridge 8 Decline 4 0.5 Ridge R2 = 0.32, p = 0.0002 1960 1980 2000 2020 1982 1999 2006 2015 Time (years) Drought event (year) (C) (D) PMARR MARR Multidecadal resistance- MAGR 1999 2006 2015 MTGR Pearson Forest area MGRT resilience curve correlation MRT 1.0 2 Valley MD Ridge MRS 0.5 PARR ARR 0.0 AGR −0.5 1 TGR GRT −1.0 RT RC RS 1990 2000 2010 T R R PA R RR MD MT RT MA R MA R PM RR R S C GRT S MG RT R R MR R TG AG AR G G AR M Time (years) Short-term and multidecadal resilience indices F I G U R E 6 Forest secondary growth and resilience trends. (A) Time series (1960–­2020) of secondary tree growth represented by Basal Area Increment (cm2/year) for the mountain ridge (violet line, declining strand) and valley areas (green line, resilient stand), showing the contrasting tendencies observed since the 1980s in climatically exposed and buffered areas of the rear edge. (B) Observed intra-­population variability in forest resilience indicators, applying in this case the multidecadal resilience index (MRs, Table 1) for the selected drought events. MRs is calculated for the tree population at the ridge (violet box plots) and valley areas (green, resilient areas of the forest). A one-­ way ANOVA test is provided comparing MRs mean between ridge and valley areas (***p < 0.001; ****p < 0.0001). In the illustrated results reference time period considered for the resilience analysis is a 2-­year period (1979–­1980). For other reference periods the results were qualitatively the same (see Supporting Information). (C) Observed trends in the multidecadal resistance–­resilience curve (MRCFW ), with a fixed window of 2 years, illustrating the diverging trends between the valley and ridge areas (green and violet lines). Thin lines illustrate the trends for individual trees in both areas. (D) Observed Pearson correlationships between all the analysed resilience indices TA B L E 2 General Linear models of secondary tree growth, testing for the effect of climatic variables. BAI (Basal Area Increment) was used as a measure of secondary tree growth. Forest area indicates the location of the stand in the topographic gradient t Time period Forest area Predictor Estimate SE value Pr (>|t|) R2 p-­value 1960–­1979 Ridge TApr 24.87 15.75 1.58 0.11 0.13 0.001 TJan 37.77 11.47 3.29 0.001 PreJun 0.81 0.50 1.61 0.11 Valley TApr 106.56 24.94 4.27 <0.001 0.30 <0.001 TJan 79.73 18.17 4.39 <0.001 PreJun 1.45 0.79 1.83 0.069 1980–­2018 Ridge TJun −35.83 8.67 −4.12 <0.001 0.64 <0.001 PreJun 1.57 0.34 4.66 <0.001 Valley PreJun 3.49 0.70 4.98 <0.001 0.40 <0.001 Bold letters highlight the most significant predictors detected. indicate that resilience trends are best captured applying multiple indi- multidecadal resilience indices and the application of multiannual re- cators. In our case study, we focused on capturing and integrating both silience curves (Figure 6; Table S9) provide useful descriptors of forest multidecadal and short-­term variability in resistance, resilience and re- vulnerability to climate warming in its intra-­population variability, com- covery capacity (Table 1, Table S6). We observed that in forest stands plementing in this way assessments based in the analysis of short-­term affected by long-­term declines and non-­stationary climatic effects, resilience indicators (Gazol et al., 2017, 2018; Lloret et al., 2011). 12 | Journal of Ecology CARNICER et al. We observed that forest resilience indices differed in areas af- Dawson, 2003). However, in our case study we found that mammal fected by contrasting levels of thermal exposure during extreme herbivory (i.e. a biotic stressor) largely determined the fine-­scale heatwave events. Thermal buffering and amplification effects patterns of a key life-­history stage (i.e. sapling recruitment). Due to during extreme heatwave events remain still poorly described in these biotic factors, recruitment success largely differed between Mediterranean ecosystems and elsewhere (Carnicer, Stefanescu, these competing tree species (Figure 5), characterised by contrast- et al., 2019; Frey et al., 2016; Schwaab et al., 2020). In this paper ing climatic niches, and by spatially overlapped distributions in a we report a first detailed description of the modulation of the narrow forest ecotone. This highlights that in order to predict the thermal exposure in tree saplings during an unprecedented ex- performance and assess the resilience of rear edge populations to treme summer heatwave event, based on detailed field measure- global change, integrative studies considering multiple drivers are ments (van Oldenborgh et al., 2019; Zhao et al., 2020). Variability more likely to capture all relevant demographic processes. in thermal exposure during the heatwave was clearly linked to We also found high spatial heterogeneity in the response to cli- local topography and vegetation cover effects (Figures 1C, 2 and matic factors in the rear edge, with positive growth responses in cli- 3). In a relatively short altitudinal distance, we measured con- matically buffered sections of the topographic gradient (valley areas), trasting microclimatic conditions. Particularly, during the heat- which contrasted with the negative trends observed in climatically wave event the highest air temperature in the understorey and exposed areas of the gradient (Figure 6). The response of most resil- the lowest air relative humidity coincided with the forest sections ience indices in Q. petraea was strikingly different across a gradient closest to the ridge (Figure 2A–­C). Microclimate variability along of tens of metres along the ecotone. Within the same population the gradient was further intensified during a summer heatwave and inside a rear edge area, climatically buffered microsites were relative to previous colder periods of the summer (Figure 1C,D). associated with greater resilience and recovery capacity (Figure 6; Because an increase in summer heatwave frequency is very likely Figure S13). Our results show that within-­population variability in for the Mediterranean region in the coming decades (Guiot & growth and resilience can be really high and that not all trees might Cramer, 2016; IPCC, 2018), the effect of topography might ex- be equally able to cope with future warmer and drier conditions. acerbate within-­site variability and modulate impacts of extreme The difference in growth and resilience between ridge and valley events on tree demographic performance (Rita et al., 2020). Q. petraea trees has been progressively exacerbated since the 1980 Climate buffering effects of forest canopies are well-­k nown (Figure 6). Trees in the valley area maintained a positive growth (De Frenne et al., 2021; Zellweger et al., 2020), although they trend, achieving growth rates also reported for other non-­stressed depend on water availability and thus can vary intra-­ and inter-­ stands at the same latitude, but presenting lower growth rates than annually (Carnicer, Stefanescu, et al., 2019; Davis et al., 2019; von those observed in northern populations of the same species (Bréda Arx et al., 2013). Microclimatic variability seems to be relevant in & Granier, 1996; Martínez-­Sancho et al., 2018). In sharp contrast, the studied rear edge forest, since Q. petraea recruitment success trees in the ridge area have been affected by an increasing water was restricted to mesic, northern slope microsites and secondary deficit, possibly due to a post-­1980 climatic shift to warmer and growth and resilience were also largely affected by topographical drier conditions reported in the Eastern Iberian Peninsula (Carnicer, variability (Figure 6). We detected greater vulnerability and lim- Domingo-­Marimon, et al., 2019). This can be concluded from the ited resilience of the tree stand situated in the higher and more climate-­growth model (Table 2), which clearly indicates more sig- exposed part of the topographic gradient, in line with previous nificant negative effects of summer drought and summer tempera- studies highlighting the key role of local-­s cale topography in deter- tures after 1980. Climate impacts on tree growth have shifted from mining thermal exposure (Bramer et al., 2018; Graae et al., 2018; positive effects of warm late winter and early spring temperatures Lenoir et al., 2013; Scherrer & Körner, 2010). In addition, we de- (1960–­1970s) to strong limitations by warm and dry summer ef- tected species-­specific cooling effects of some species (Q. ilex) in fects over the last decades (Table 2). These climate-­induced trends microclimatic measurements (Table S8; Figure S6). This suggests are fully consistent with previous studies of secondary drought in that the current expansion of treelets of Q. ilex and other simi- Q. petraea reporting negative effects of summer drought stress and lar oak species in the understorey of most Mediterranean forests, positive effects of late winter and early spring temperatures (Bréda which is linked to the widespread advance of forest succession & Granier, 1996; Michelot, Bréda, et al., 2012; Michelot, Simard, over extensive areas (Carnicer et al., 2014), might be altering the et al., 2012). The secondary growth in this species starts prior to patterns of microclimatic exposure at the regional scale. budburst and leaf expansion (Michelot, Bréda, et al., 2012; Michelot, Our results indicate that rear edge tree populations are vulner- Simard, et al., 2012); and almost half of the annual ring can be pro- able to the combined effect of both global change-­induced abiotic duced in the spring period (Bréda & Granier, 1996). Although warmer and biotic stressors. Rear edge populations often occur in climate temperatures can promote cambium cell division (Pérez-­de-­Lis refugia, but their persistence can be often restricted or facilitated by et al., 2017), this process is also limited by water availability (Eilmann biotic interactions (Hampe, 2004; Hampe & Jump, 2011). Among abi- et al., 2006). Increased resilience observed in valley trees could be otic factors, climate is often recognised as the most relevant factor facilitated by topography, promoting soil water storage in lowland limiting species distributions, including trailing and rear edge areas areas, and dampening the negative effects of drought for the val- (Harsch & HilleRisLambers, 2016; Normand et al., 2009; Pearson & ley trees (Dymond et al., 2017). Similarly, shallow soils in ridge areas CARNICER et al. Journal of Ecology | 13 cannot store as much water as soil further down the valley, promot- populations is a major concern and challenging task in the NE Iberian ing limited resilience and recovery capacity in these stands. Resilient Peninsula (Giménez-­Anaya et al., 2020), and increasing hunting pres- trees in the valley appear to be less sensitive to climate than other sure in protected areas often fails to effectively regulate population populations located at the core of the species distribution (Martínez-­ size (Figure S9). Similarly, in the rest of Europe, wild boar popula- Sancho et al., 2018). At least a part of the studied rear edge forest tions have grown at rates not matched by human hunting pressure appears to benefit from buffering environmental conditions that de- (González-­Crespo et al., 2018; Massei et al., 2015). As a result, mul- couple this stand from macroclimatic fluctuations, matching the defi- tiple ecosystem alterations are currently caused by increasing wild nition of a climatic and hydrological refugia (Fan et al., 2020; Lenoir boar and roe deer impacts (Barrios-­Garcia & Ballari, 2012). Previous et al., 2017; McLaughlin et al., 2017). In this sense, our findings are studies have also shown a limiting effect of wild boar damage on in line with those of recent studies reporting higher resistance to the recruitment of Quercus saplings (Gómez et al., 2003; Gómez & drought in rear edge populations, compared to core distribution Hódar, 2008). populations (Muffler et al., 2020; Vilà-­C abrera et al., 2019). Figure 7 summarises possible future scenarios for Mediterranean Our findings do not support scenarios predicting rapid, spatially rear edge forest resilience in the next decades, accounting for the homogeneous distributional shifts and limited resilience in rear edge mixed effects of climate change warming impacts, altered mammal populations of the Mediterranean biome. These uniform and simple herbivory and microclimatic variability. The current scenario in our scenarios of massive forest retreat can be often predicted applying study case is dominated low wild predator pressure and suboptimal climate-­envelope species distribution models (SDMs, as discussed in human hunting pressure, which ultimately result in a high mammal Hampe, 2004; Hampe & Jump, 2011). Our results instead are more herbivory pressure and recruitment failure (Figure 7A). These trends supportive of more complex scenarios including spatially heteroge- co-­occur with still moderate but increasing impacts of warming, neous responses, characterised with contrasting intra-­population which affect the forest resilience in non-­buffered areas of the forest. trends of forest resilience. Such heterogeneity would be mediated We suggest that in the next decades, resilient scenarios will possibly by local thermal buffering and amplification effects, and variable require effective mitigation and adaptation policies to limit climate drought stress linked to topographic gradients, soil and land cover warming impacts on forests (IPCC, 2018), combined with regional variability (Carnicer, Stefanescu, et al., 2019). We provide evidence plans for wild mammal management (Figure 7B, violet arrows). In the for a significant within population variation in resilience, emphasising eastern Iberian Peninsula wolve populations are expected to natu- in our case study the relevance of microclimatic variability in the dy- rally recover in the following decades, from present very low den- namics of rear edge populations. Furthermore, our results also sup- sities in the area [in the order of 0.0004 ind./km2, which contrast port the existence of non-­stationary effects of climate determining with current high regional densities of wild boar (3–­15 ind./km2) and forest growth responses (Table 2; Astigarraga et al., 2020; Carnicer, roe deer (0.5–­10 ind./km2); Rosell et al., 2019]. Ineffective imple- Domingo-­Marimon, et al., 2019; Peltier & Ogle, 2020). Finally, as pre- mentation of climate policies and regional wild fauna management viously mentioned, our results also emphasise an important role of plans would possibly lead to low resilience scenarios for protected other interacting global change drivers such as imbalanced recruit- and endangered rear edge forests, increasing negative interactions ment due to increased herbivore pressure in locally altered, defau- between unregulated wild mammal populations and warming im- nated trophic networks, that we discuss in detail below. pacts on protected areas and on multiple land sectors in adjacent Significant species-­specific differences were found in the pat- areas (Figure 7C; Badano et al. 2015; Havlík et al., 2014; Kolström terns of tree sapling abundance and recruitment success. Although et al., 2011; Kuijper et al., 2019; Martin et al., 2020; Nardone both Q. petraea and Q. ilex presented similar recruitment success for et al., 2010; Rojas-­Downing et al., 2017). According to these sce- the first size class (R1), a larger survival rate in the transient phase narios, we suggest that forest resilience studies would benefit from from the first to the second class was observed in the evergreen combining the analysis of growth resilience relationships and indi- Q. ilex species (R2, Figures 4 and 5). This was possibly related to the cators of demographic resilience in different tree life-­history stages strong and selective herbivore pressure by roe deer specifically af- (reviewed in Capdevila et al., 2020). fecting broadleaved oak saplings of Q. petraea (Figure 5C), which is We conclude that forest resilience and heterogeneous demo- consistent with previous studies of selective consumption of tree graphic responses to climate warming are ultimately modulated by saplings and dietary preferences in this mammal species (e.g. Partl the influence of the local environment on individual trees (Levine et al., 2002). et al., 2016). Accurate predictions of forest responses to changes in We also observed a spatially clustered pattern of damage by wild climate would largely benefit, thus, from the integration of local-­scale boars, which was mainly concentrated nearby river bed and forest variation of key factors such as microclimatic conditions and biotic trail areas (Figure S10). In these specific areas, the effects of root- interactions, as our results and previous studies support (Lembrechts ing damage by wild boar were dominant and often completely pre- et al., 2019; Vilà-­Cabrera & Jump, 2019; Zellweger et al., 2020). cluded recruitment success in sessile oaks. Wild boar and roe deer Macroclimate should not be used alone to predict recruitment suc- are the main ungulate species inhabiting the studied area (Rosell cess and forest resilience in future scenarios because fine-­scale mi- et al., 2019; Torres et al., 2016). Controlling roe deer and wild boar croclimate and biotic interactions can be also important factors. 14 | Journal of Ecology CARNICER et al. Future scenarios More resilient Emerging (A) Current scenario (B) Active-management (D) Tree growth Resilient tree growth and Impacts on carbon defaunation and Historical carbon uptake Recovered carbon uptake Mitigation/ stocks, protected forest Low wild High mammal predator Low mammal Adaptation regeneration dynamics predator pressure pressure herbivore failure herbivore policies and ecosystem services pressure pressure 1 ºC 1.5–3 ºC Moderate Increased global wild Increased impacts global warming warming impacts (1 ºC) RCP 2.6/4.5 W/m 2 fauna in food Low hunting Optimal management agrosystems pressure hunting variability pressure variability Impacts on forest Reduced tree Active wild fauna management growth and plantations carbon uptake Recruitment failure Recovered Variable, predator herbivore pressure pressure cycles Impacts on livestock apex predators 3–5 ºC High global Natural dynamics, without warming human-induced population RCP 6/8.5 W/m2 Suboptimal Interactions and parallel control hunting variability impacts of wild mammal pressure Climate-change impacts fauna and global warming in multiple sectors (C) management and mitigation F I G U R E 7 A synthetic summary of multiple scenarios for forest resilience, explicitly considering the complex interactions between global warming impacts, wild fauna management and local microclimatic variability. (A) A synthetic diagram of the current scenario in the studied Quercus petraea rear edge forest. (B, C) Future scenarios ranging from low to high atmospheric concentrations of greenhouse gases (Representative Concentration Pathways, IPCC, 2018). The width of the green and orange arrows in panels A, B and C represents the relative magnitude of the expected effects of global warming and microclimatic buffering processes. (D) Emerging combined impacts of global warming and unregulated wild mammal populations on protected forests and adjacent land areas. The arrows in panel D indicate the key ecological and climatic processes involved, according to the grey panel on the left AC K N OW L E D G E M E N T S ORCID Lorena Gómez-­Aparicio, Jonathan Lenoir, Iain Stott and two anony- Jofre Carnicer https://orcid.org/0000-0001-7454-8296 mous reviewers provided comments and scientific advice that largely Maria Vives-­Ingla https://orcid.org/0000-0003-4887-8392 improved the manuscript. Thanks to Angel Miño, Vicenç Bros, Javi Laura Blanquer https://orcid.org/0000-0002-7459-2075 Sala, Daniel Pons, Sonia Llobet, and Parc Natural de Sant Llorenç del Santi Sabaté https://orcid.org/0000-0003-1854-0761 Munt for technical and logistic support. This research was funded by Emília Gutiérrez https://orcid.org/0000-0002-6085-5700 PID2020-­117636GB-­C21 and GGL-­2016-­78093-­R grants. Teresa Sauras https://orcid.org/0000-0003-1294-6410 Josep Peñuelas https://orcid.org/0000-0002-7215-0150 AU T H O R S ' C O N T R I B U T I O N S Adrià Barbeta https://orcid.org/0000-0002-8357-1719 J.C. designed the research; X.M.-­C., J.C., L.B., A.B., M.V.-­I., C.R., S.S., E.G., J.P. and T.S. collected the field data; M.V.-­I. and J.C. designed and REFERENCES implemented the resilience calculations; J.C., M.V.-­I. and A.B. contrib- Arend, M., Kuster, T., Günthardt-­Goerg, M. S., & Dobbertin, M. (2011). uted new analytic tools and analysed the data; J.C. wrote an initial Provenance-­specific growth responses to drought and air warming in three European oak species (Quercus robur, Q. petraea and Q. pubes- draft of the paper and all authors contributed to the final document. cens). Tree Physiology, 31(3), 287–­297. https://doi.org/10.1093/treep​ hys/tpr004 PEER REVIEW Astigarraga, J., Andivia, E., Zavala, M. A., Gazol, A., Cruz-­A lonso, V., The peer review history for this article is available at https://publo​ Vicente-­S errano, S. M., & Ruiz-­B enito, P. (2020). Evidence of non-­ stationary relationships between climate and forest responses: ns.com/publo​n/10.1111/1365-­2745.13752. Increased sensitivity to climate change in Iberian forests. Global Change Biology, 26(9), 5063–­5 076. https://doi.org/10.1111/ DATA AVA I L A B I L I T Y S TAT E M E N T gcb.15198 Data is available from Dryad Digital Repository: https://doi. Badano, E. I., Samour-­Nieva, O. R., Flores, J., & Douterlungne, D. (2015). Microclimate and seeding predation as drivers of tree recruitment in org/10.5061/dryad.vq83b​k3t5 (Carnicer, 2021). CARNICER et al. Journal of Ecology | 15 human-­disturbed oak forests. Forest Ecology and Management, 356, Castagneri, D., Carrer, M., Regev, L., & Boaretto, E. (2020). Precipitation 93–­100. https://doi.org/10.1016/j.foreco.2015.07.031 variability differently affects radial growth, xylem traits and ring po- Barrios-­Garcia, M. N., & Ballari, S. A. (2012). Impact of wild boar (Sus rosity of three Mediterranean oak species at xeric and mesic sites. scrofa) in its introduced and native range: A review. Biological Invasions, Science of the Total Environment, 699, 134285. 14(11), 2283–­2300. https://doi.org/10.1007/s1053​0 -­012-­0229-­6 Castro, J., Gómez, J. M., García, D., Zamora, R., & Hódar, J. A. (1999). Seed Batllori, E., Lloret, F., Aakala, T., Anderegg, W. R. L., Aynekulu, E., predation and dispersal in relict Scots pine forests in southern Spain. Plant Bendixsen, D. P., Bentouati, A., Bigler, C., Burk, C. J., Camarero, Ecology, 145(1), 115–­123. https://doi.org/10.1023/A:10098​65703906 J. J., Colangelo, M., Coop, J. D., Fensham, R., Floyd, M. L., Galiano, Castro, J., Zamora, R., Hódar, J. A., & Gómez, J. M. (2004). Seedling L., Ganey, J. L., Gonzalez, P., Jacobsen, A. L., Kane, J. M., … Zeeman, establishment of a boreal tree species (Pinus sylvestris) at its south- B. (2020). Forest and woodland replacement patterns following ernmost distribution limit: Consequences of being in a marginal drought-­related mortality. Proceedings of the National Academy of Mediterranean habitat. Journal of Ecology, 266–­277. https://doi. Sciences of the United States of America, 117, 29720–­29729. https:// org/10.1111/j.0022-­0 477.2004.00870.x doi.org/10.1073/pnas.20023​14117 Coll, M., Peñuelas, J., Ninyerola, M., Pons, X., & Carnicer, J. (2013). Bou, J., & Vilar, L. (2019). Sessile oak forest plant community changes Multivariate effect gradients driving forest demographic responses on the NE Iberian Peninsula over recent decades. Journal of Plant in the Iberian Peninsula. Forest Ecology and Management, 303, 195–­ Ecology, 12(5), 894–­906. https://doi.org/10.1093/jpe/rtz029 209. https://doi.org/10.1016/j.foreco.2013.04.010 Bramer, I., Anderson, B. J., Bennie, J. J., Bladon, A. J., De Frenne, P., Corcuera, L., Camarero, J. J., & Gil-­Pelegrín, E. (2002). Functional groups Hemming, D. L., Hill, R., Kearney, M. R., Körner, C., Korstjens, A. H., in Quercus species derived from the analysis of pressure-­volume Lenoir, J., Maclean, I., Marsh, C. D., Morecroft, M. D., Ohlemüller, R., curves. Trees -­Structure and Function, 16(7), 465–­472. https://doi. Slater, H. D., Suggitt, A. J., Zellweger, F., & Gillingham, P. K. (2018). org/10.1007/s0046​8-­0 02-­0187-­1 Advances in monitoring and modelling climate at ecologically rele- Côté, S. D., Rooney, T. P., Tremblay, J. P., Dussault, C., & Waller, D. M. vant scales. Advances in Ecological Research, 58, 101–­161. https://doi. (2004). Ecological impacts of deer overabundance. Annual Review org/10.1016/BS.AECR.2017.12.005 of Ecology, Evolution, and Systematics, 35, 113–­147. https://doi. Bréda, N., & Granier, A. (1996). Intra-­and interannual variations of tran- org/10.1146/annur​ev.ecols​ys.35.021103.105725 spiration, leaf area index and radial growth of a sessile oak stand Davis, K. T., Dobrowski, S. Z., Holden, Z. A., Higuera, P. E., & Abatzoglou, J. (Quercus petraea). In Annales des sciences forestières (Vol. 53, No. 2–­3, T. (2019). Microclimatic buffering in forests of the future: The role of pp. 521–­536). EDP Sciences. local water balance. Ecography, 42(1), 1–­11. https://doi.org/10.1111/ Bueno, C. G., Alados, C. L., Gómez-­García, D., Barrio, I. C., & García-­ ecog.03836 González, R. (2009). Understanding the main factors in the extent and De Frenne, P., Lenoir, J., Luoto, M., Scheffers, B. R., Zellweger, F., Aalto, distribution of wild boar rooting on alpine grasslands. Journal of Zoology, J., Ashcroft, M. B., Christiansen, D. M., Decocq, G., De Pauw, K., 279(2), 195–­202. https://doi.org/10.1111/j.1469-­7998.2009.00607.x Govaert, S., Greiser, C., Gril, E., Hampe, A., Jucker, T., Klinges, D. H., Bunn, A. G., Korpela, M., Biondi, F., Campelo, F., Mérian, P., Qeadan, F., Koelemeijer, I. A., Lembrechts, J. J., Marrec, R., Meeussen, C., Ogée, Zang, C., Buras, A., Cecile, J., Mudelsee, M., & Schulz, M. (2020). J., Tyystjärvi, V., Vangansbeke, P., & Hylander, K. (2021). Forest Package ‘dplR’. R Package Version 1.7.2. microclimates and climate change: Importance, drivers and future Cairns, D. M., & Moen, J. (2004). Herbivory influences tree lines. Journal research agenda. Global Change Biology, 27(11), 2279–­2297. https:// of Ecology, 92(6), 1019–­1024. https://doi.org/10.1111/j.1365-­2745.​ doi.org/10.1111/gcb.15569 2004.00945.x Den Ouden, J., Jansen, P. A., & Smit, R. (2005). Jays, mice and oaks: preda- Canadell, J. G., & Jackson, R. B. (Eds.). (2021). Ecosystem collapse and tion and dispersal of Quercus robur and Q. petraea in North-­western climate change: An introduction. In Ecosystem collapse and climate Europe. Seed Fate: Predation, Dispersal, and Seedling Establishment, 13, change (pp. 1–­9). Springer. 223–­239. Capdevila, P., Stott, I., Beger, M., & Salguero-­Gómez, R. (2020). Towards a Dittmar, C., Zech, W., & Elling, W. (2003). Growth variations of Common comparative framework of demographic resilience. Trends in Ecology & beech (Fagus sylvatica L.) under different climatic and environmen- Evolution, 35, 776–­785. https://doi.org/10.1016/j.tree.2020.05.001 tal conditions in Europe –­ A dendroecological study. Forest Ecology Carnicer, J., Coll, M., Pons, X., Ninyerola, M., Vayreda, J., & Peñuelas, and Management, 173(1–­3), 63–­78. https://doi.org/10.1016/S0378​ J. (2014). Large-­scale recruitment limitation in Mediterranean pines: -­1127(01)00816​-­7 The role of Quercus ilex and forest successional advance as key re- Dobrowski, S. Z., Swanson, A. K., Abatzoglou, J. T., Holden, Z. A., Safford, gional drivers. Global Ecology and Biogeography, 23(3), 371–­384. H. D., Schwartz, M. K., & Gavin, D. G. (2015). Forest structure and https://doi.org/10.1111/geb.12111 species traits mediate projected recruitment declines in western Carnicer, J., Domingo-­Marimon, C., Ninyerola, M., Camarero, J. J., Bastos, US tree species. Global Ecology and Biogeography, 24(8), 917–­927. A., López-­Parages, J., Blanquer, L., Rodríguez-­Fonseca, B., Lenton, https://doi.org/10.1111/geb.12302 T. M., Dakos, V., Ribas, M., Gutiérrez, E., Peñuelas, J., & Pons, X. (2019). Dymond, S. F., Bradford, J. B., Bolstad, P. V., Kolka, R. K., Sebestyen, S. Regime shifts of Mediterranean forest carbon uptake and reduced D., & DeSutter, T. M. (2017). Topographic, edaphic, and vegetative resilience driven by multidecadal ocean surface temperatures. Global controls on plant-­available water. Ecohydrology, 10(8), 1–­12. https:// Change Biology, 25(8), 2825–­2840. https://doi.org/10.1111/gcb.14664 doi.org/10.1002/eco.1897 Carnicer, J., Stefanescu, C., Vives-­Ingla, M., López, C., Cortizas, S., Eilmann, B., Weber, P., Rigling, A., & Eckstein, D. (2006). Growth reac- Wheat, C., Vila, R., Llusià, J., & Peñuelas, J. (2019). Phenotypic bio- tions of Pinus sylvestris L. and Quercus pubescens Willd. to drought markers of climatic impacts on declining insect populations: A key years at a xeric site in Valais, Switzerland. Dendrochronologia, 23(3), role for decadal drought, thermal buffering and amplification effects 121–­132. https://doi.org/10.1016/j.dendro.2005.10.002 and host plant dynamics. Journal of Animal Ecology, 88(3), 376–­391. Estes, J. A., Terborgh, J., Brashares, J. S., Power, M. E., Berger, J., Bond, https://doi.org/10.1111/1365-­2656.12933 W. J., Carpenter, S. R., Essington, T. E., Holt, R. D., Jackson, J. B. C., Carnicer, J., Wheat, C., Vives-­Ingla, M., Ubach, A., Domingo-­Marimón, C., Marquis, R. J., Oksanen, L., Oksanen, T., Paine, R. T., Pikitch, E. K., Nylin, S., Stefanescu, C., Vila, R., Wiklund, C., & Peñuelas, J. (2017). Ripple, W. J., Sandin, S. A., Scheffer, M., Schoener, T. W., … Wardle, D. Evolutionary responses of invertebrates to global climate change: A. (2011). Trophic downgrading of planet earth. Science, 333(6040), The role of life-­history trade-­offs and multidecadal climate shifts. In 301–­3 06. https://doi.org/10.1126/scien​ce.1205106 S. N. Johnson & T. H. Jones (Eds.), Global climate change and terrestrial Fan, B., Tao, W., Qin, G., Hopkins, I., Zhang, Y., Wang, Q., Lin, H., & Guo, invertebrates (pp. 317–­3 48). Wiley. L. (2020). Soil micro-­climate variation in relation to slope aspect, 16 | Journal of Ecology CARNICER et al. position, and curvature in a forested catchment. Agricultural and Hampe, A., & Petit, R. J. (2005). Conserving biodiversity under climate Forest Meteorology, 290. https://doi.org/10.1016/j.agrfo​rmet.2020.​ change: The rear edge matters. Ecology Letters, 8(5), 461–­467. https:// 107999 doi.org/10.1111/j.1461-­0248.2005.00739.x Frey, S. J. K., Hadley, A. S., Johnson, S. L., Schulze, M., Jones, J. A., & Harsch, M. A., & HilleRisLambers, J. (2016). Climate warming and sea- Betts, M. G. (2016). Spatial models reveal the microclimatic buffering sonal precipitation change interact to limit species distribution shifts capacity of old-­growth forests. Science Advances, 2(4). https://doi. across Western North America. PLoS ONE, 11(7), 1–­17. https://doi. org/10.1126/sciadv.1501392 org/10.1371/journ​al.pone.0159184 Gazol, A., Camarero, J. J., Anderegg, W. R. L., & Vicente-­Serrano, S. M. Havlík, P., Valin, H., Herrero, M., Obersteiner, M., Schmid, E., Rufino, M. (2017). Impacts of droughts on the growth resilience of Northern C., Mosnier, A., Thornton, P. K., Böttcher, H., Conant, R. T., Frank, Hemisphere forests. Global Ecology and Biogeography, 26(2), 166–­176. S., Fritz, S., Fuss, S., Kraxner, F., & Notenbaert, A. (2014). Climate https://doi.org/10.1111/geb.12526 change mitigation through livestock system transitions. Proceedings Gazol, A., Camarero, J. J., Vicente-­Serrano, S. M., Sánchez-­Salguero, R., of the National Academy of Sciences, 111(10), 3709–­3714. https://doi. Gutiérrez, E., de Luis, M., Sangüesa-­Barreda, G., Novak, K., Rozas, org/10.1073/pnas.13080​4 4111 V., Tíscar, P. A., Linares, J. C., Martín-­Hernández, N., Martínez del Herrera, C. M., & Bazaga, P. (2008). Adding a third dimension to the edge Castillo, E., Ribas, M., García-­González, I., Silla, F., Camisón, A., of a species' range: Altitude and genetic structuring in mountainous Génova, M., Olano, J. M., … Galván, J. D. (2018). Forest resilience landscapes. Heredity, 100(3), 275–­285. https://doi.org/10.1038/ to drought varies across biomes. Global Change Biology, 24(5), 2143–­ sj.hdy.6801072 2158. https://doi.org/10.1111/gcb.14082 Holmes, R. L. (1983). Computer-­assisted quality controlin tree-­ring dat- Giménez-­Anaya, A., Bueno, C. G., Fernández-­Llario, P., Fonseca, C., ing and measurement. Tree-­Ring Bulletin, 43(1), 69–­78. García-­González, R., Herrero, J., & Rosell, C. (2020). What do Ingrisch, J., & Bahn, M. (2018). Towards a comparable quantification of we know about wild boar in Iberia? In F. M. Angelici & L. Rossi resilience. Trends in Ecology & Evolution, 33(4), 251–­259. https://doi. (Eds.), Problematic wildlife II (pp. 251–­271). Springer. https://doi. org/10.1016/j.tree.2018.01.013 org/10.1007/978-­3-­030-­42335​-­3 IPCC. (2018). Global Warming of 1.5°C. An IPCC Special Report on the im- Gómez, J. M., García, D., & Zamora, R. (2003). Impact of vertebrate pacts of global warming of 1.5°C above pre-­industrial levels and related acorn-­and seedling-­predators on a Mediterranean Quercus pyrenaica global greenhouse gas emission pathways, in the context of strength- forest. Forest Ecology and Management, 180(1–­3), 125–­134. https:// ening the global response to the threat of climate change, sustainable doi.org/10.1016/S0378​-­1127(02)00608​- ­4 development, and efforts to eradicate poverty [V. Masson-­Delmotte, Gómez, J. M., & Hódar, J. A. (2008). Wild boars (Sus scrofa) affect the P. Zhai, H.-­O. Pörtner, D. Roberts, J. Skea, P. R. Shukla, A. Pirani, W. recruitment rate and spatial distribution of holm oak (Quercus ilex). Moufouma-­Okia, C. Péan, R. Pidcock, S. Connors, J. B. R. Matthews, Forest Ecology and Management, 256(6), 1384–­1389. https://doi. Y. Chen, X. Zhou, M. I. Gomis, E. Lonnoy, T. Maycock, M. Tignor, & org/10.1016/j.foreco.2008.06.045 T. Waterfield (Eds.)]. World Meteorological Organization Technical González-­Crespo, C., Serrano, E., Cahill, S., Castillo-­Contreras, R., Document. Cabañeros, L., López-­Martín, J. M., Roldán, J., Lavín, S., & López-­ Jump, A. S., & Peñuelas, J. (2006). Genetic effects of chronic habitat Olvera, J. R. (2018). Stochastic assessment of management strate- fragmentation in a wind-­pollinated tree. Proceedings of the National gies for a Mediterranean peri-­urban wild boar population. PLoS ONE, Academy of Sciences, 103(21), 8096–­8100. 13(8), 1–­19. https://doi.org/10.1371/journ​al.pone.0202289 Kolström, M., Lindner, M., Vilén, T., Maroschek, M., Seidl, R., Lexer, M. Graae, B. J., Vandvik, V., Armbruster, W. S., Eiserhardt, W. L., Svenning, J., Netherer, S., Kremer, A., Delzon, S., Barbati, A., Marchetti, M., & J.-­C ., Hylander, K., Ehrlén, J., Speed, J. D. M., Klanderud, K., Bråthen, Corona, P. (2011). Reviewing the science and implementation of cli- K. A., Milbau, A., Opedal, Ø. H., Alsos, I. G., Ejrnæs, R., Bruun, H. H., mate change adaptation measures in European forestry. Forests, 2(4), Birks, H. J. B., Westergaard, K. B., Birks, H. H., & Lenoir, J. (2018). 961–­982. https://doi.org/10.3390/f2040961 Stay or go–­how topographic complexity influences alpine plant pop- Körner, C., & Paulsen, J. (2004). A world-­wide study of high altitude ulation and community responses to climate change. Perspectives treeline temperatures. Journal of Biogeography, 31(5), 713–­732. in Plant Ecology, Evolution and Systematics, 30, 41–­50. https://doi. https://doi.org/10.1111/j.1365-­2699.2003.01043.x org/10.1016/j.ppees.2017.09.008 Kuijper, D. P. J., Churski, M., Trouwborst, A., Heurich, M., Smit, C., Kerley, Granda, E., Alla, A. Q., Laskurain, N. A., Loidi, J., Sánchez-­Lorenzo, A., & G. I. H., & Cromsigt, J. P. G. M. (2019). Keep the wolf from the door: Camarero, J. J. (2018). Coexisting oak species, including rear edge How to conserve wolves in Europe's human-­dominated landscapes? populations, buffer climate stress through xylem adjustments. Tree Biological Conservation, 235, 102–­111. https://doi.org/10.1016/j. Physiology, 38(2), 159–­172. https://doi.org/10.1093/treep​hys/tpx157 biocon.2019.04.004 Granda, E., Escudero, A., & Valladares, F. (2014). More than just Lembrechts, J. J., Nijs, I., & Lenoir, J. (2019). Incorporating microcli- drought: Complexity of recruitment patterns in Mediterranean for- mate into species distribution models. Ecography, 42(7), 1267–­1279. ests. Oecologia, 176(4), 997–­1007. https://doi.org/10.1007/s0044​ https://doi.org/10.1111/ecog.03947 2-­014-­3 064-­x Lenoir, J., Graae, B. J., Aarrestad, P. A., Alsos, I. G., Armbruster, W. S., Guiot, J., & Cramer, W. (2016). Climate change: The 2015 Paris Agreement Austrheim, G., Bergendorff, C., Birks, H. J., Bråthen, K. A., Brunet, J., thresholds and Mediterranean basin ecosystems. Science, 354(6311), Bruun, H. H., Dahlberg, C. J., Decocq, G., Diekmann, M., Dynesius, M., 465–­468. https://doi.org/10.1126/scien​ce.aah5015 Ejrnaes, R., Grytnes, J. A., Hylander, K., Klanderud, K., … Svenning, J. Hacket-­Pain, A. J., & Friend, A. D. (2017). Increased growth and reduced C. (2013). Local temperatures inferred from plant communities sug- summer drought limitation at the southern limit of Fagus sylvatica L., gest strong spatial buffering of climate warming across Northern despite regionally warmer and drier conditions. Dendrochronologia, Europe. Global Change Biology, 19(5), 1470–­1481. 44, 22–­3 0. https://doi.org/10.1016/j.dendro.2017.02.005 Lenoir, J., Hattab, T., & Pierre, G. (2017). Climatic microrefugia under an- Hampe, A. (2004). Bioclimate envelope models: What they detect and thropogenic climate change: Implications for species redistribution. what they hide. Global Ecology and Biogeography, 13(5), 469–­471. Ecography, 40, 253–­266. https://doi.org/10.1111/ecog.02788 https://doi.org/10.1111/j.1466-­822X.2004.00090.x Levine, N. M., Zhang, K., Longo, M., Baccini, A., Phillips, O. L., Lewis, S. Hampe, A., & Jump, A. S. (2011). Climate relicts: Past, present, future. L., Alvarez-­Dávila, E., De Andrade, A. C. S., Brienen, R. J. W., Erwin, Annual Review of Ecology, Evolution, and Systematics, 42(1), 313–­333. T. L., Feldpausch, T. R., Mendoza, A. L. M., Vargas, P. N., Prieto, A., https://doi.org/10.1146/annur​ev-­ecols​ys-­10271​0 -­145015 Silva-­Espejo, J. E., Malhi, Y., & Moorcroft, P. R. (2016). Ecosystem CARNICER et al. Journal of Ecology | 17 heterogeneity determines the ecological resilience of the Amazon to towards the dry distribution margin of European beech. Journal of climate change. Proceedings of the National Academy of Sciences of the Biogeography, 47(9), 1910–­1921. https://doi.org/10.1111/jbi.13884 United States of America, 113(3), 793–­797. https://doi.org/10.1073/ Munier, A., Hermanutz, L., Jacobs, J. D., & Lewis, K. (2010). The inter- pnas.15113​4 4112 acting effects of temperature, ground disturbance, and herbivory on Lloret, F., & Batllori, E. (2021). Climate-­induced global forest shifts due to seedling establishment: Implications for treeline advance with cli- heatwave-­drought. In J. G. Canadell & R. B. Jackson (Eds.), Ecosystem mate warming. Plant Ecology, 210(1), 19–­3 0. https://doi.org/10.1007/ collapse and climate change (pp. 155–­186). Springer. s1125​8-­010-­9724-­y Lloret, F., Escudero, A., Iriondo, J. M., Martínez-­V ilalta, J., & Valladares, Nardone, A., Ronchi, B., Lacetera, N., Ranieri, M. S., & Bernabucci, U. F. (2012). Extreme climatic events and vegetation: The role of stabi- (2010). Effects of climate changes on animal production and sus- lizing processes. Global Change Biology, 18(3), 797–­8 05. https://doi. tainability of livestock systems. Livestock Science, 130(1–­3), 57–­69. org/10.1111/j.1365-­2486.2011.02624.x https://doi.org/10.1016/j.livsci.2010.02.011 Lloret, F., Keeling, E. G., & Sala, A. (2011). Components of tree re- Nogués-­Bravo, D., Simberloff, D., Rahbek, C., & Sanders, N. J. (2016). silience: Effects of successive low-­growth episodes in old pon- Rewilding is the new Pandora's box in conservation. Current Biology, derosa pine forests. Oikos, 120(12), 1909–­1920. https://doi. 26(3), R87–­R91. https://doi.org/10.1016/j.cub.2015.12.044 org/10.1111/j.1600-­0706.2011.19372.x Normand, S., Treier, U. A., Randin, C., Vittoz, P., Guisan, A., & Svenning, J. Lloret, F., Peñuelas, J., & Estiarte, M. (2005). Effects of vegetation canopy C. (2009). Importance of abiotic stress as a range-­limit determinant and climate on seedling establishment in Mediterranean shrubland. for European plants: Insights from species responses to climatic gra- Journal of Vegetation Science, 16(1), 67–­76. https://doi.org/10.1111/ dients. Global Ecology and Biogeography, 18(4), 437–­4 49. https://doi. j.1654-­1103.2005.tb023​39.x org/10.1111/j.1466-­8238.2009.00451.x Martin, J. L., Chamaillé-­Jammes, S., & Waller, D. M. (2020). Deer, Partl, E., Szinovatz, V., Reimoser, F., & Schweiger-­Adler, J. (2002). Forest wolves, and people: Costs, benefits and challenges of living to- restoration and browsing impact by roe deer. Forest Ecology and gether. Biological Reviews, 95(3), 782–­8 01. https://doi.org/10.1111/ Management, 159(1–­2), 87–­100. https://doi.org/10.1016/S0378​ brv.12587 -­1127(01)00712​-­5 Martínez-­Sancho, E., Dorado-­Liñán, I., Gutiérrez Merino, E., Matiu, M., Pearson, R. G., & Dawson, T. P. (2003). Predicting the impacts of climate Helle, G., Heinrich, I., & Menzel, A. (2018). Increased water-­use effi- change on the distribution of species: Are bioclimate envelope mod- ciency translates into contrasting growth patterns of Scots pine and els useful? Global Ecology and Biogeography, 12(5), 361–­371. https:// sessile oak at their southern distribution limits. Global Change Biology, doi.org/10.1046/j.1466-­822X.2003.00042.x 24(3), 1012–­1028. https://doi.org/10.1111/gcb.13937 Peltier, D. M., & Ogle, K. (2020). Tree growth sensitivity to climate is Massei, G., & Genov, P. V. (2004). The environmental impact of wild boar. temporally variable. Ecology Letters, 23(11), 1561–­1572. https://doi. Galemys, 16(1), 135–­145. org/10.1111/ele.13575 Massei, G., Kindberg, J., Licoppe, A., Gačić, D., Šprem, N., Kamler, J., Pérez-­de-­Lis, G., Olano, J. M., Rozas, V., Rossi, S., Vázquez-­Ruiz, R. A., Baubet, E., Hohmann, U., Monaco, A., Ozoliņš, J., Cellina, S., Podgórski, & García-­González, I. (2017). Environmental conditions and vas- T., Fonseca, C., Markov, N., Pokorny, B., Rosell, C., & Náhlik, A. (2015). cular cambium regulate carbon allocation to xylem growth in Wild boar populations up, numbers of hunters down? A review of deciduous oaks. Functional Ecology, 31(3), 592–­603. https://doi. trends and implications for Europe. Pest Management Science, 71(4), org/10.1111/1365-­2435.12789 492–­500. https://doi.org/10.1002/ps.3965 Príncipe, A., Matos, P., Sarris, D., Gaiola, G., do Rosário, L., Correia, O., & Mayer, J. J., Nelson, E. A., & Wike, L. D. (2000). Selective depredation Branquinho, C. (2019). In Mediterranean drylands microclimate af- of planted hardwood seedlings by wild pigs in a wetland resto- fects more tree seedlings than adult trees. Ecological Indicators, 106, ration area. Ecological Engineering, 15(Suppl. 1), 79–­85. https://doi. https://doi.org/10.1016/j.ecoli​nd.2019.105476. 105476 org/10.1016/S0925​-­8574(99)00074​-­9 Pulido, F., Castagneyrol, B., Rodríguez-­Sánchez, F., Cáceres, Y., Pardo, McLaren, B. E., & Peterson, R. O. (1994). Wolves, moose, and tree rings A., Moracho, E., Kollmann, J., Valladares, F., Ehrlén, J., Jump, A. S., on Isle Royale. Science, 266, 1555–­1558. https://doi.org/10.1126/ Svenning, J.-­C ., & Hampe, A. (2019). Asymmetry in marginal popu- scien​ce.266.5190.1555 lation performance foreshadows widespread species range shifts. McLaughlin, B. C., Ackerly, D. D., Klos, P. Z., Natali, J., Dawson, T. E., BioRxiv, 529560. https://doi.org/10.1101/529560 & Thompson, S. E. (2017). Hydrologic refugia, plants, and climate R Core Team. (2020). R: A language and environment for statistical comput- change. Global Change Biology, 23(8), 2941–­2961. ing. R Foundation for Statistical Computing. Mérian, P., & Lebourgeois, F. (2011). Consequences of decreasing the Ripple, W. J., & Beschta, R. L. (2007). Hardwood tree decline following number of cored trees per plot on chronology statistics and climate-­ large carnivore loss on the Great Plains, USA. Frontiers in Ecology growth relationships: A multispecies analysis in a temperate climate. and the Environment, 5(5), 241–­246. https://doi.org/10.1890/1540-­ Canadian Journal of Forest Research, 41(12), 2413–­2422. https://doi. 9295(2007)5%5B241​:HTDFL​C%5D2.0.CO;2 org/10.1139/X11-­149 Rita, A., Camarero, J. J., Nolè, A., Borghetti, M., Brunetti, M., Pergola, N., Michelot, A., Bréda, N., Damesin, C., & Dufrêne, E. (2012). Differing Serio, C., Vicente-­Serrano, S. M., Tramutoli, V., & Ripullone, F. (2020). growth responses to climatic variations and soil water deficits of The impact of drought spells on forests depends on site conditions: Fagus sylvatica, Quercus petraea and Pinus sylvestris in a temperate The case of 2017 summer heatwave in southern Europe. Global forest. Forest Ecology and Management, 265, 161–­171. https://doi. Change Biology, 26(2), 851–­863. https://doi.org/10.1111/gcb.14825 org/10.1016/j.foreco.2011.10.024 Roemer, G. W., Gompper, M. E., & van Valkenburgh, B. (2009). The eco- Michelot, A., Simard, S., Rathgeber, C., Dufrêne, E., & Damesin, C. (2012). logical role of the mammalian mesocarnivore. BioScience, 59(2), 165–­ Comparing the intra-­annual wood formation of three European 173. https://doi.org/10.1525/bio.2009.59.2.9 species (Fagus sylvatica, Quercus petraea and Pinus sylvestris) as re- Rojas-­Downing, M. M., Nejadhashemi, A. P., Harrigan, T., & Woznicki, lated to leaf phenology and non-­structural carbohydrate dynamics. S. A. (2017). Climate change and livestock: Impacts, adaptation, Tree Physiology, 32, 1033–­1045. https://doi.org/10.1093/treep​hys/ and mitigation. Climate Risk Management, 16, 145–­163. https://doi. tps052 org/10.1016/j.crm.2017.02.001 Muffler, L., Weigel, R., Hacket-­Pain, A. J., Klisz, M., van der Maaten, Rosell, C., Pericas, B., Colomer, J., & Navàs, F. (2019). Program for moni- E., Wilmking, M., Kreyling, J., & van der Maaten-­Theunissen, M. toring populations of boar in Catalonia. 2018-­2019 season. Generalitat (2020). Lowest drought sensitivity and decreasing growth synchrony de Catalunya. 18 | Journal of Ecology CARNICER et al. Saez-­Royuela, C., & Telleria, J. L. (1986). The increased population of M., Villanúa, D., Rosell, C., … Acevedo, P. (2019). Harmonization of the Wild Boar (Sus scrofa L.) in Europe. Mammal Review, 16, 97–­101. the use of hunting statistics for wild boar density estimation in dif- https://doi.org/10.1111/j.1365-­2907.1986.tb000​27.x ferent study areas. EFSA Supporting Publications, 16(9). https://doi. Scherrer, D., & Körner, C. (2010). Infra-­red thermometry of alpine land- org/10.2903/sp.efsa.2019.en-­1706 scapes challenges climatic warming projections. Global Change Biology, Vicente-­Serrano, S. M., Lopez-­Moreno, J. I., Beguería, S., Lorenzo-­L acruz, 16(9), 2602–­2613. https://doi.org/10.1111/j.1365-­2486.2009.02122.x J., Sanchez-­Lorenzo, A., García-­Ruiz, J. M., Azorin-­Molina, C., Morán-­ Schwaab, J., Davin, E. L., Bebi, P., Duguay-­Tetzlaff, A., Waser, L. T., Haeni, Tejeda, E., Revuelto, J., Trigo, R., Coelho, F., & Espejo, F. (2014). M., & Meier, R. (2020). Increasing the broad-­leaved tree fraction Evidence of increasing drought severity caused by temperature rise in European forests mitigates hot temperature extremes. Scientific in southern Europe. Environmental Research Letters, 9(4), https://doi. Reports, 1–­9, https://doi.org/10.1038/s4159​8-­020-­71055​-­1 org/10.1088/1748-­9326/9/4/044001 Schwarz, J., Skiadaresis, G., Kohler, M., Kunz, J., Schnabel, F., Vitali, Vilà-­C abrera, A., & Jump, A. S. (2019). Greater growth stability of trees V., & Bauhus, J. (2020). Quantifying growth responses of trees to in marginal habitats suggests a patchy pattern of population loss and drought—­A critique of commonly used resilience indices and recom- retention in response to increased drought at the rear edge. Ecology mendations for future studies. Current Forestry Reports. https://doi. Letters, 22(9), 1439–­1448. https://doi.org/10.1111/ele.13329 org/10.1007/s4072​5-­020-­0 0119​-­2 Vilà-­C abrera, A., Premoli, A. C., & Jump, A. S. (2019). Refining predictions Shestakova, T. A., Aguilera, M., Ferrio, J. P., Gutiérrez, E., & Voltas, J. of population decline at species' rear edges. Global Change Biology, (2014). Unravelling spatiotemporal tree-­ring signals in Mediterranean 25(5), 1549–­1560. https://doi.org/10.1111/gcb.14597 oaks: A variance-­covariance modelling approach of carbon and ox- Vives-­Ingla, M., Sala-­Garcia, J., Stefanescu, C., Peñuelas, J., & Carnicer, ygen isotope ratios. Tree Physiology, 34(8), 819–­838. https://doi. J. (2020). Varying thermal exposure, host-­plant traits and ovipo- org/10.1093/treep​hys/tpu037 sition behaviour across vegetation ecotones. bioRxiv. https://doi. Smit, C., Ruifrok, J. L., van Klink, R., & Olff, H. (2015). Rewilding with large org/10.1101/2020.02.11.944439 herbivores: The importance of grazing refuges for sapling establish- Von Arx, G., Graf Pannatier, E., Thimonier, A., & Rebetez, M. (2013). ment and wood-­pasture formation. Biological Conservation, 182, 134–­ Microclimate in forests with varying leaf area index and soil mois- 142. https://doi.org/10.1016/j.biocon.2014.11.047 ture: Potential implications for seedling establishment in a chang- Sousa, P. M., Barriopedro, D., Ramos, A. M., García-­Herrera, R., Espírito-­ ing climate. Journal of Ecology, 101(5), 1201–­1213. https://doi. Santo, F., & Trigo, R. M. (2019). Saharan air intrusions as a relevant org/10.1111/1365-­2745.12121 mechanism for Iberian heatwaves: The record breaking events Young, H. S., McCauley, D. J., Galetti, M., & Dirzo, R. (2016). Patterns, of August 2018 and June 2019. Weather and Climate Extremes, 26, causes, and consequences of anthropocene defaunation. Annual https://doi.org/10.1016/j.wace.2019.100224 Review of Ecology, Evolution, and Systematics, 47(1), 333–­358. https:// Spinoni, J., Barbosa, P., De Jager, A., McCormick, N., Naumann, G., Vogt, doi.org/10.1146/annur​ev-­ecols​ys-­11241​4-­054142 J. V., Magni, D., Masante, D., & Mazzeschi, M. (2019). A new global Zellweger, F., De Frenne, P., Lenoir, J., Vangansbeke, P., Verheyen, K., database of meteorological drought events from 1951 to 2016. Bernhardt-­Römermann, M., Baeten, L., Hédl, R., Berki, I., Brunet, J., Journal of Hydrology: Regional Studies, 22, https://doi.org/10.1016/j. Van Calster, H., Chudomelová, M., Decocq, G., Dirnböck, T., Durak, ejrh.2019.100593 T., Heinken, T., Jaroszewicz, B., Kopecký, M., Máliš, F., … Coomes, D. Tegel, W., Seim, A., Hakelberg, D., Hoffmann, S., Panev, M., Westphal, T., (2020). Forest microclimate dynamics drive plant responses to warming. & Büntgen, U. (2014). A recent growth increase of European beech Science, 368(6492), 772–­775. https://doi.org/10.1126/scien​ce.aba6880 (Fagus sylvatica L.) at its Mediterranean distribution limit contradicts Zhao, W., Zhou, N., & Chen, S. (2020). The record-­breaking high tem- drought stress. European Journal of Forest Research, 133(1), 61–­71. perature over Europe in June of 2019. Atmosphere, 11(5), 16–­20. https://doi.org/10.1007/s1034​2-­013-­0737-­7 https://doi.org/10.3390/atmos​11050524 Thurm, E. A., Uhl, E., & Pretzsch, H. (2016). Mixture reduces climate sen- sitivity of Douglas-­fir stem growth. Forest Ecology and Management, 376, 205–­220. https://doi.org/10.1016/j.foreco.2016.06.020 S U P P O R T I N G I N FO R M AT I O N Torres, R. T., Carvalho, J., Fonseca, C., Serrano, E., & López-­Martín, J. Additional Supporting Information may be found online in the M. (2016). Longterm assessment of roe deer reintroductions in Supporting Information section. NorthEast Spain: A case of success. Mammalian Biology, 81(4), 415–­ 422. https://doi.org/10.1016/j.mambio.2016.05.002 Turco, M., Von Hardenberg, J., AghaKouchak, A., Llasat, M. C., Provenzale, A., & Trigo, R. M. (2017). On the key role of droughts in the dynamics How to cite this article: Carnicer, J., Vives-­Ingla, M., of summer fires in Mediterranean Europe. Scientific Reports, 7(1), 1–­ Blanquer, L., Méndez-­C amps, X., Rosell, C., Sabaté, S., 10. https://doi.org/10.1038/s4159​8-­017-­0 0116​-­9 Gutiérrez, E., Sauras, T., Peñuelas, J., & Barbeta, A. (2021). Van Oldenborgh, G. J., Philip, S., Kew, S., Vautard, R., Boucher, O., Otto, Forest resilience to global warming is strongly modulated by F., Haustein, K., Soubeyroux, J.-­M., Ribes, A., & Robin, Y. (2019). Human contribution to the record-­ breaking June 2019 heatwave in local-­scale topographic, microclimatic and biotic conditions. France (p. 32). World Weather Attribution (WWA). Journal of Ecology, 00, 1–­18. https://doi. Vicente, J., Palencia, P., Plhal, R., Blanco-­A guiar, J. A., Laguna, E., Soriguer, org/10.1111/1365-­2745.13752 R., López, J. F., Podgórski, T., Petrović, K., Apollonio, M., Scandura, M., Ferroglio, E., Zanet, S., Brivio, F., Keuling, O., Smith, G. C., Guibert,