PNAS PLUS

Climate extremes and predicted warming threaten Mediterranean Holocene firs forests refugia Raúl Sánchez-Salgueroa,b,1, J. Julio Camareroa, Marco Carrerc, Emilia Gutiérrezd, Arben Q. Allae, Laia Andreu-Haylesf,g, Andrea Heviah, Athanasios Koutavasi, Elisabet Martínez-Sanchoj, Paola Nolak, Andreas Papadopoulosl, Edmond Pashoe, Ervin Toromanie, José A. Carreiram, and Juan C. Linaresb a Instituto Pirenaico de Ecología–Consejo Superior de Investigaciones Científicas, 50192 Zaragoza, Spain; bDepartamento Sistemas Físicos, Químicos y Naturales, Universidad Pablo de Olavide, 41013 Sevilla, Spain; cDipartimento Territorio e Sistemi Agro-Forestali, Università degli Studi di Padova, Legnaro 35020, Italy; dDepartamento Biologia Evolutiva, Ecologia i Ciències Ambientals, University of Barcelona, 08028 Barcelona, Spain; eFakulteti i Shkencave Pyjore, Universiteti Bujqësor i Tiranës, 1029 Tirana, Albania; fTree-Ring Laboratory, Lamont-Doherty Earth Observatory, Palisades, NY 10964; gInstitut Català de Ciències del Clima, 08005 Barcelona, Spain; hForest and Wood Technology Research Centre, 33936 Siero, Asturias, Spain; iDepartment of Engineering Science and Physics, College of Staten Island, City University of New York, Staten Island, NY 10314; jDepartment of Ecology and Ecosystem Management, Technische Universität München, 85354 Freising, Germany; kDipartimento Scienze della Terra e dell’Ambiente, Università degli Studi di Pavia, 27100 Pavia, Italy; l Department of Forestry and Natural Environment Management, Technological Educational Institute of Stereas Elladas, 36100 Karpenissi, Greece; and mCentro de Estudios Avanzados de la Tierra, Universidad de Jaén, 23071 Jaén, Spain

|

Abies spp. climate change forward growth model

warming (11). The intense habitat fragmentation and degradation as a result of anthropogenic activities has reduced their population sizes in the Modern Era and increased the extinction risk for some of these taxa (12–14). More severe droughts could reduce the present ranges of Circum-Mediterranean firs, particularly the dry (southernmost or lowermost) distribution limits, but also in xeric sites (15). This calls for a better assessment of the CMFF vulnerability to projected climate warming and extreme climatic events. Comparing forest growth responses to observed and forecast climate can provide the right avenue to quantify CMFF vulnerability. The susceptibility of trees to extreme climatic events such as droughts is usually expressed by low radial-growth rates beyond a critical minimum threshold (4). Many investigations confirmed that extensive growth reduction can be considered as an early warning signal of stand vulnerability to episodes of drought- and/or heattriggered dieback (5, 16–20) with, in the most critical circumstances, cascade effects on local contraction of the species distribution area (21–24). Time series of tree-ring width indices (TRWi) coupled Significance

| dendroecology | emission scenarios |

Climate extremes are major drivers of long-term forest growth trends, but we still lack appropriate knowledge to anticipate their effects. Here, we apply a conceptual framework to assess the vulnerability of Circum-Mediterranean Abies refugia in response to climate warming, droughts, and heat waves. Using a tree-ring network and a process-based model, we assess the future vulnerability of Mediterranean Abies forests. Models anticipate abrupt growth reductions for the late 21st century when climatic conditions will be analogous to the most severe dry/ heat spells causing forest die-off in the past decades. However, growth would increase in moist refugia. Circum-Mediterranean fir forests currently subjected to warm and dry conditions will be the most vulnerable according to the climate model predictions for the late 21st century.

C

limate is a major driver of long-term changes and shifts in tree species distribution (1). Forecasts for the 21st century predict an increase of up to 6 °C in mean temperature, as well as more severe and lasting dry spells and heat waves (2). Climate extremes can play a key role, albeit little explored, in affecting the functioning and vulnerability of drought-prone relict forests (3) and in triggering forest die-off (4) with cascading effects on biodiversity, carbon, water, and nutrient cycling, and ultimately on the provisioning of forest ecosystem services (5). The persistence of tree species under new climatic conditions will depend upon their ability to acclimatize or buffer climate extremes (6). Special attention needs to be paid to tree-species refugia, i.e., sites where the species persisted during glaciations and currently form relict or isolated stands (7). Relict fir populations located in Circum-Mediterranean fir forests (CMFF), where water availability is the major forest growth constraint, could be threatened by the forecasted warmer and drier conditions (8). In addition, drought stress could also be amplified by forecasted hot temperature extremes (9). Understanding how these populations would respond to the predicted climate warming and amplified drought stress is relevant to forecasting how and where biodiversity refugia would endure adverse climates extremes during the 21st century (10). Circum-Mediterranean firs (Abies species) constitute an example of endangered forests and an ecosystem highly vulnerable to climate

www.pnas.org/cgi/doi/10.1073/pnas.1708109114

Author contributions: R.S.-S., J.J.C., and J.C.L. designed research; R.S.-S., J.J.C., M.C., E.G., A.H., and J.C.L. performed research; R.S.-S., J.J.C., A.H., and J.C.L. contributed new reagents/analytic tools; R.S.-S., J.J.C., M.C., E.G., A.Q.A., L.A.-H., A.K., E.M.-S., P.N., A.P., E.P., E.T., J.A.C., and J.C.L. contributed to field data; R.S.-S., J.J.C., A.H., and J.C.L. analyzed data; and R.S.-S., J.J.C., M.C., E.G., A.Q.A., L.A.-H., A.H., A.K., E.M.-S., P.N., A.P., E.P., E.T., J.A.C., and J.C.L. wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. Published under the PNAS license. Data deposition: The datasets reported in this paper have been deposited in the Pangaea repository, https://doi.org/10.1594/PANGAEA.882101. 1

To whom correspondence should be addressed. Email: [email protected].

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10. 1073/pnas.1708109114/-/DCSupplemental.

PNAS Early Edition | 1 of 9

ENVIRONMENTAL SCIENCES

Warmer and drier climatic conditions are projected for the 21st century; however, the role played by extreme climatic events on forest vulnerability is still little understood. For example, more severe droughts and heat waves could threaten quaternary relict tree refugia such as Circum-Mediterranean fir forests (CMFF). Using tree-ring data and a process-based model, we characterized the major climate constraints of recent (1950–2010) CMFF growth to project their vulnerability to 21st-century climate. Simulations predict a 30% growth reduction in some fir species with the 2050s businessas-usual emission scenario, whereas growth would increase in moist refugia due to a longer and warmer growing season. Fir populations currently subjected to warm and dry conditions will be the most vulnerable in the late 21st century when climatic conditions will be analogous to the most severe dry/heat spells causing dieback in the late 20th century. Quantification of growth trends based on climate scenarios could allow defining vulnerability thresholds in tree populations. The presented predictions call for conservation strategies to safeguard relict tree populations and anticipate how many refugia could be threatened by 21st-century dry spells.

ECOLOGY

Edited by William H. Schlesinger, Cary Institute of Ecosystem Studies, Millbrook, NY, and approved September 29, 2017 (received for review May 16, 2017)

with process-based forward models of growth offer a valuable tool for understanding forest growth responses to climate change (25). We projected the radial growth of CMFF (SI Appendix, Fig. S1 and Table S1) as a function of different 21st-century climate forecasts under two contrasting representative emission scenarios (Representative Concentration Pathways, or RCPs) (2, 26), applying the Vaganov–Shashkin–Lite model (hereafter VS-Lite model) (27, 28). We analyzed climate-growth associations paying special attention to the stability of Mediterranean fir refugia under climate extremes (see the conceptual framework in Fig. 1).

Assessment of Extremes Using Process-Based Growth Models. The VS-Lite model accurately tracked the year-to-year growth variability (Fig. 3 and SI Appendix, Figs. S6 and S7 and Tables S3 and S4). The growth response to temperature (gT) peaked during the growing season (May to September) in strictly Mediterranean firs (15), while the growth response to moisture (gM) dropped during summer in response to dry conditions (Fig. 3). Tree growth was limited by low temperatures (gT < gM) at the beginning and end of the growing season and by soil moisture availability (gM < gT) from late spring to autumn.

Results and Discussion We detected a generally positive influence on tree growth by previous late-summer wet and cool climatic conditions but a high site-to-site variability in the climate-growth relationships within each species (SI Appendix, Fig. S2). Years with significant growth reduction (negative pointer years) usually corresponded to heat waves and dry spells (SI Appendix, Figs. S3 and S4). Forecasted warmer and drier conditions corresponding to the RCP 8.5 scenario projected a growth decline at lower latitudes and elevations (Fig. 2 and SI Appendix, Fig. S3 and Table S2), although still with a high inter- and intraspecific variability.

Changes in the Climatic Thresholds of Circum-Mediterranean Firs. Under the warmest climate scenario (RCP 8.5), several Mediterranean firs (Figs. 4 and 5) are predicted to be increasingly constrained by soil moisture availability, i.e., more drought stress throughout the growing season (Figs. 4 and 5). Some fir species will experience persistent dry conditions at levels similar to those during 20thcentury extreme dry spells (Figs. 3–5 and SI Appendix, Fig. S8). On the contrary, the wettest sites are predicted to have an extended growing season (SI Appendix, Fig. S9). We provide empirically based constraints on modeled climateinduced changes in growth for CMFF and present an approach

Fig. 1. Conceptual framework for assessing forest vulnerability to climate change based on growth responses to extreme climate events. The figure shows the framework stages (Left column) and the procedures and tools used to fulfill this approach (Right column). (A) We quantified, for each site, how growth responds to climatic conditions and extreme climatic events using correlation coefficients between the mean series of TRWi and mean temperature (red) or total precipitation (blue) (1). Temporal window spans depend on the species and location [e.g., from previous (t−1) September up to current (t) October]. Significant (P < 0.05) values are located outside the dashed lines. We analyzed the impact of critical climate variables on growth using extreme pointer year analyses (2). Extreme negative growth years (very low TRWi values) corresponding to >50% raw growth reduction in more than 50% of trees on a site were selected (red squares). (B) To define the vulnerability threshold, we simulated TRWi using the VS-Lite process-based growth model (3) in which several parameters are estimated: temperature (T1) or soil moisture (M1) thresholds below which growth will not occur, optimal temperature (T2), or soil moisture (M2) above which growth is not limited (Materials and Methods; B, 3 is adapted from ref. 76). Growth depends on temperature (gT, red lines and red areas correspond to the mean and SD, respectively) and soil moisture limitations (gM, blue lines and blue areas correspond to the mean and SD, respectively). Note that high-gT or low-gM values indicate low and high growth limitations, respectively. Observed years with extremely low growth (e.g., 1995) based on observed drought- or heat-induced dieback episodes are indicated as vulnerability threshold (Vi) (critical tipping thresholds of growth stability are indicated by black dashed lines) (4). To model growth responses as a function of climate extremes and to define vulnerability thresholds, we compared observed and simulated (VS-Lite) TRWi series (5). (C) We compared observed and projected TRWi based on stepwise multiple linear regressions using emission and climate scenarios as predictors (6). Simulated (VS-Lite) growth responses were compared with extreme years under future climate scenarios (2050–2100) (7). The observed gT and gM values defined the observed vulnerability thresholds for each forest refugia comparing projected growth responses with forecasted climate scenarios (8).

2 of 9 | www.pnas.org/cgi/doi/10.1073/pnas.1708109114

Sánchez-Salguero et al.

Fig. 2. Projected ring-width indices (TRWi) for each fir species and site considering the warmest emission scenario (RCP 8.5). Species climate-based models means are shown as black lines, and gray areas represent the maximum and minimum TRWi values considering all sites for each species. The lower plots of each panel show the assessment of the impact on growth stability (sensitivity) of the occurrence of extreme forecasted 21st-century events by TRWi reductions (GR) computed considering the 5-y moving average values from 2011 to 2100. The gray-scale bars indicate TRWi reductions while red bars indicate selected climatically extreme years causing reduced growth in each species during the 21st century. Bold lines frame extreme values (TRWi reduction >25%). See SI Appendix, Fig. S5, for RCP 4.5 scenario and SI Appendix, Table S1, for site codes.

Sánchez-Salguero et al.

PNAS Early Edition | 3 of 9

PNAS PLUS ECOLOGY

drought resistance in some species (38) but not necessarily with a corresponding effect on growth (39, 40). Following a drought avoidance strategy (15), A. cephalonica copes better with long-term droughts than A. borisii-regis and can thus be found on a greater variety of substrates, including dry compact limestones. In contrast, A. borisii-regis occurs in more humid sites and at higher latitudes (41) (SI Appendix, Table S5). Our analyses suggest that for A. borisii-regis higher spring temperatures will negatively affect tree growth characteristics (Fig. 5). Something similar is currently observed in A. alba, where increased solar radiation accompanied by low rainfall and warmer temperatures had a negative effect on growth as less sunshine seemed to be sufficient for an adequate photosynthetic carbon assimilation in this shade-tolerant species (42). In the case of A. cilicica and A. pinsapo under a hypothetical moderate temperature increase (RCP 4.5 scenario) and if rainfall remained stable, the species’ growth responses would not be greatly disrupted (SI Appendix, Figs. S5 and S8). In contrast, a temperature increase in the previous autumn and winter may make these species more sensitive to cold spells in winter (15) (SI Appendix, Table S6). Warm autumns may also lead to a depletion of reserves as a result of the prevalence of respiration over photosynthesis having a negative impact on growth in the next year (43). Mediterranean firs show adaptive features to increasing drought stress by adjusting their phenology and hydraulic architecture (e.g., leaf area to sapwood area ratio) to local climatic conditions, particularly the most widespread species (44). The wide intersite variability in modeled VS-Lite parameters (Fig. 5 and SI Appendix, Table S4) may arise from parallel intersite differences in environmental conditions or from contrasting interspecies drought tolerance caused by local genetic adaptation and phenotypic plasticity (45). Anatomical and physiological studies suggest that, for some of these species, there is an early

ENVIRONMENTAL SCIENCES

that accounts for intraspecific traits and geographic shifts in climate response to forecasted 21st-century extreme events (Fig. 1). This represents an important advance in quantifying ecological responses to future climatic constraints since conservation of marginal tree populations will depend not only on climate forcing, but also on adequately considering vulnerability thresholds (29– 31). Our projections question the role of long-term evolutionary buffering effects on some relict Circum-Mediterranean fir species and highlight the impacts of forecasted 21st-century climate variability on the stability of climatic tree refugia (32) (Fig. 5). We projected stable growth conditions for some species (e.g., A. alba and A. borisii-regis) in the northern Mediterranean Basin under the moderate RCP 4.5 warming scenario, but a substantial decline in southern populations under the significantly warmer RCP 8.5 (Figs. 4 and 5 and SI Appendix, Fig. S8). On the contrary, growth enhancement due to a longer growing season is predicted at high-elevation (wettest) sites (33, 34) (Fig. 2 and SI Appendix, Table S2). We also expect that some species (e.g., A. cephalonica, A. pinsapo, and A. cilicica) growing in dry and low-elevation sites may be strongly and negatively influenced by warming-induced drought stress and will be the most sensitive to climate extremes under the RCP 8.5 scenario (Fig. 3 and SI Appendix, Fig. S3). This implies that some forest stands will be unlikely to keep pace with the enhanced decline forecast for the late 21st century, especially given that several episodes of long-lasting reductions in growth, linked to subsequent tree mortality (17, 20, 35), have already been observed in different dry regions of southern and central Europe (19, 20, 36). Our findings support the predictions of range contractions, local extinctions, and species composition changes in many fir stands due to warming temperatures (23) and the increasing frequency of drought-triggered dieback (5, 21, 22, 24, 37). Within this scenario, the increase in water-use efficiency related to rising atmospheric CO2 concentrations could ameliorate

Fig. 3. Simulated monthly growth response curves (gT, gM) using the VS-Lite model for the period 1950–2010. The growth responses consider temperature (gT, red lines and red areas correspond to the mean and SD, respectively) and soil moisture limitations (gM, blue lines and blue areas correspond to the mean and SD, respectively) for each species. See Fig. 1 and SI Appendix, Table S1, for site codes. Note that low-gM values indicate growth limitations. Selected extremely low growth years are indicated by dashed lines in different colors: green (previous year) and black (current year). The map displays site-level Pearson correlation coefficients computed by the VS-Lite growth model between observed and predicted series of ring-width indices. See SI Appendix, Fig. S7, for more locations.

closure of stomata that provides an efficient reduction of water loss in the case of water shortage (46, 47). However, as already observed, the loss in hydraulic performance during extreme events can translate into growth decline (19), especially given that high temperatures in combination with strong radiation can also negatively influence soil water balance (48). Unexpectedly, the highest vulnerability of fir species is recorded on soils with higher waterholding capacity, showing that soil features need to be explicitly considered to properly assess tree vulnerability and growth potential (49). For example, the strict water-stress avoidance strategy of A. pinsapo and A. cephalonica to limit water losses is reflected in rapid stomatal closure at rather positive stem water potentials and relatively high values of soil water content (15, 47). Stability of Climate Refugia Under Different Warming Scenarios.

Under the strongest temperature increases, the VS-lite growth model projected population shifts from moisture-limited to also warming-limited conditions during the growing season (Figs. 3 and 4). This would occur mainly in lower and dry sites where spring– summer conditions are presently close to the species’ limits (Fig. 5) and where the future climate will be similar to current extreme (drought/heat spells) events (Figs. 3 and 5). In general, climate projections predict shorter growing seasons in southern A. alba and lower elevation A. borisii-regis, A. cephalonica, A. cilicica, and A. pinsapo sites, but longer ones in the wettest and highelevation A. alba sites (25) (Figs. 3 and 4). This would imply a higher growth sensitivity to climate during the early growing season, suggesting a prominent role of spring water deficit as already observed in xylogenesis investigations (50). On the contrary, the minimum temperature threshold for growth (T1) and the temperature at which growth is not limited (T2) will increase for all CMFF (SI Appendix, Tables S4 and S7), with a potential increased frequency of evapotranspiration stress during the predicted shorter growing season (Fig. 5) (4). Thus, the forecasted drier and hotter growing 4 of 9 | www.pnas.org/cgi/doi/10.1073/pnas.1708109114

season would negatively affect the growth of CMFF stands at lower elevation and xeric sites by exceeding their functional thresholds for optimal growth responses (i.e., gT and gM) (Figs. 4 and 5). Since all of the CMFF experienced a fairly stable thermal range during most of the Holocene (51), an increase in water and thermal stress during the second half of the 21st century may have dreadful effects on the relict populations at the warm (lower) edges of the species’ ranges (Fig. 5). Assuming a low genetic variability within these species (see ref. 52), few fir taxa may potentially withstand a temperature increase of about 4.8 °C in xeric sites (Fig. 4). Our results thus suggest that some of the quaternary populations’ refugia at the warmer end of the species distribution (e.g., A. pinsapo, A. tazaotana, or low elevation A. cephalonica and A. cilicica) may be threatened within a climate warming scenario (Fig. 5 and SI Appendix, Table S7). Evaluating the climatic limits of species distributions over time is key to understanding the species responses to future climate change. Space–time extended datasets, such as those derived from tree rings, are needed to reconstruct species resilience to past extreme events and to predict future dieback processes (cf. 18 and 53). Our methodological framework is robust and provides results that are in accordance with what has been observed in field studies on forest vulnerability and climate-driven growth thresholds (15, 19). Nevertheless, we captured only part of the key drivers governing forest growth, and future simulations should also include information on shoot and leaf growth (54), demographic processes (55), and tree acclimation and disturbances (56). Conclusions. Our modeling approach helps to highlight locally resistant tree populations, providing drought-heat tolerance that can be potentially preadapted to future climatic conditions expected in more northerly locations. Having these sources for assisted migration can be highly relevant for the application of conservation strategies under global change or habitat fragmentation. Finally, Sánchez-Salguero et al.

PNAS PLUS ECOLOGY

growth analyses should be ideally supplemented by the genetic background of the species when stakeholders have to decide about management but also to confer resilience to forests. Process-based models allowed us to define the vulnerability thresholds of CMFF refugia by quantifying growth responses and identifying vulnerability thresholds to forecasted climate warming and more extreme dry spells (Fig. 1). We disentangled the respective roles of the different species- and site-specific factors that will predispose to growth decline and local extinction processes. The predictions under the business-as-usual scenario (RCP 8.5) forecasted a growth reduction in some southern forests. This knowledge can support forest managers to better prepare specific strategies to cope with forecasted climatic challenges. Materials and Methods Study Sites and Sampling Design. The study area includes the main ranges of most Circum-Mediterranean fir (Abies) species (CMFF) (SI Appendix, Fig. S1 and Table S1) (11). Fir species appear in the wettest sites of the Mediterranean Basin usually dominating in N-oriented sites with deep soils (57). They can be found on different parent materials, calcareous or noncalcareous, but grow best on deeper acid soils with high water reserves (15). These species develop in a subhumid climate where annual precipitation is between 600 and 1,000 mm, and mean annual temperature ranges from 6.6 to 17.2 °C (SI Appendix, Table S5). Thirty forests that were sampled across the Mediterranean region are distributed as follows: A. alba Miller, six sites in Italy (58), six sites in Spain (25); A. borisii-regis Mattfeld (five sites) (41, 48); A. cephalonica Loudon (three sites) (38, 41); A. cilicica (Antoine and Kotschy) Carrière (five sites); A. pinsapo Boissier (three sites); A. pinsapo var. maroccana (Trabut) Ceballos and Bolaños (one site); and A. pinsapo var. tazaotana (Cózar ex Huguet del Villar) Pourtet (one site) (SI Appendix, Fig. S1 and Table S1). We sampled these sites following dendrochronological methods and covering most of the geographical distribution of Abies species in the Mediterranean, thus capturing most of the ecological variability experienced by these species (SI Appendix, Fig. S1) (15). Most sampled sites are in protected areas (e.g., national parks), which guarantee that trees have been less

Sánchez-Salguero et al.

exposed to disturbances (e.g., logging or fire) than in nonprotected areas at least for the past 50–60 y. At each site, at least 12 dominant or codominant trees (when they were available) separated by at least 20 m from each other were cored at 1.3 m using Pressler increment borers on the cross-slope sides of the trunk whenever possible. For each study site, latitude, longitude, and mean elevation as well as topographic (slope and aspect) variables were recorded (SI Appendix, Table S1). Radial growth was measured in two to three radial cores per tree. Wood samples were sanded until rings were visible and then visually cross-dated. Once dated, treering widths were measured to the nearest 0.01 mm using a binocular microscope and a LINTAB measuring device (Rinntech). The accuracy of visual cross-dating and measurements was checked with the COFECHA program, which calculates moving correlations between each individual tree-ring series and the mean site series (59). Climate-Growth Associations. The homogenized and quality-checked CRU T.S. 3.23 dataset (www.cru.uea.ac.uk/data) was used for the period 1950–2010, providing a reliable climate data source across all of the study sites (60). This dataset contains monthly mean temperature and precipitation data gridded at a 0.5° spatial resolution that have been checked for homogeneity. To quantify climate-growth associations, we calculated mean tree-ring width series at the site scale (site chronology) (SI Appendix, Table S1). TRWi were calculated, adjusting negative exponential or linear functions and 30-y-long splines to obtain detrended growth series and corrected age– size trends. These splines allowed high-frequency (annual to decadal) growth variability to be preserved. We applied autoregressive models to model and eliminate most of the temporal (usually of first order) autocorrelation. We obtained residual or prewhitened TRWi series for each tree as ratios between the measurement and the fitted curve. Finally, we averaged the individual growth-index series into site-level chronologies following a hierarchical approach from tree to site. Among the developed chronologies, we considered only those covering the period 1950–2010, which corresponded to the period of most reliable climate data in the study area. To characterize these site TRWi chronologies, we also calculated several dendrochronological statistics (SI Appendix, Table S1) (61). The relationships between monthly climate data (mean temperatures and precipitation) and TRWi were assessed by calculating bootstrapped Pearson

PNAS Early Edition | 5 of 9

ENVIRONMENTAL SCIENCES

Fig. 4. Projected monthly growth response curves using VS-Lite model for the 2050–2100 period under the RCP 8.5 scenario. The response curves consider temperature (gT, red lines and red areas correspond to the mean and SD, respectively) and soil moisture limitations (gM, blue lines and blue areas correspond to the mean and SD, respectively) for each site (Fig. 1). Selected extreme low growth years are indicated by dashed lines in different colors: green (previous year) and black (current year) (Fig. 2 and SI Appendix, Fig. S5). The map symbols show projected growth trends (Pearson correlation coefficients of site mean TRWi series) considering the 2011–2100 period following growth projections (SI Appendix, Table S2). Correlation values higher than j0.25j are significant at P < 0.05. See SI Appendix, Table S1 for site codes and SI Appendix, Fig. S8, for more locations.

Fig. 5. Simulated (VS-Lite model) mean growth responses for the growing season and climate extreme years. The selected periods represent spring (March to May) temperature (gT, red circles indicate the mean and SD) and summer (June to September) soil moisture limitations (gM, blue triangles indicate the mean and SD) for each Abies species. See SI Appendix, Fig. S9 and Table S1 for site codes. Values were fitted for one observed (1950–2010) and two projected 21stcentury periods (2011–2049 and 2050–2100) based on RCP 8.5 emission scenario. Selected observed growth responses during extreme events (1950– 2010 period; see Figs. 1–3) with extremely low growth rates are indicated in different colors on the Left for each site. See also Figs. 2 and 3. These values defined the observed thresholds of vulnerability for each climate refugia compared with the forecasted growth response for the RCP 8.5 climate scenario. Different letters indicate significant (P < 0.05) differences between observed and projected periods based on Tukey’s HSD post hoc tests (lowercase and uppercase letters correspond to the 2011–2049 and 2050–2100 periods, respectively).

correlation coefficients for the common period 1950–2010. The temporal window of growth-climate comparisons included from the previous September to the current October. All stages of chronologies building and further analyses were performed using the R statistical software (62) and the package dplR (63). This first step of our conceptual framework defines the critical growing period at site and species levels and which climate variables (i.e., monthly temperature and precipitation) affect growth (Fig. 1). Climate Extremes and Pointer Years. As the term “climate extreme” usually defines both “extreme weather” and “extreme climate” events (cf. 63), we considered 27 climate indices (11 for precipitation and 16 for temperature)

6 of 9 | www.pnas.org/cgi/doi/10.1073/pnas.1708109114

from daily observations of precipitation and maximum and minimum temperature (ETCCDI database) calculated for each site by using the RClimDex package (see etccdi.pacificclimate.org/indices_def.shtml and SI Appendix, Table S8). RClimDex was used to compute climate indices from the observed and forecasted ensemble climate data (65–67). Spearman correlation coefficients were calculated for each site between TRWi and annual climate indices considering the previous and current year of tree-ring formation. This allowed the presence of any carryover effect of climate extremes to the following tree-ring formation to be assessed (SI Appendix, Table S9) (68). Vulnerability to climate change during these extreme climate events was assessed using extreme pointer year analysis. This was performed to determine

Sánchez-Salguero et al.

Growth Projections Based on Process-Based Modeling of Tree Growth. We used the VS-Lite model and a Bayesian parameter estimation approach to simulate TRWi as a function of climate (Fig. 1) (76–78). The model uses the Leaky Bucket Model of hydrology provided by the National Oceanic and Atmospheric Administration Climate Prediction Center (79) to estimate monthly soil moisture from temperature and total precipitation data. Snow dynamics are not explicitly considered in the model, and thus all precipitation is assumed to be liquid. For each year, the model simulates standardized tree-ring width anomalies from the minimum of the monthly growth responses to temperature (gT) and moisture (gM), modulated by insolation (gE). Day length is determined from site latitude and does not vary from year to year. The growth response functions for temperature (gT) and moisture (gM) in VS-Lite involve only two parameters. The first represents the temperature (T1) or moisture (M1) threshold below which growth will not occur, while the second is the optimal temperature (T2) or moisture (M2) above which growth is not limited by climate (Fig. 1). The growth function parameters were estimated for each site via Bayesian calibration. This scheme assumes uniform priors for the growth response parameters and independent, normally distributed errors for the modeled TRWi values. The posterior median for each parameter was used to obtain the “calibrated” growth response for a given site. Finally, the model was run over the entire period 1950–2010 using the calibrated parameters for each site to produce a simulated tree-ring chronology (TRWiVSL) that represents an estimate of the site climate signal of forest growth. A more detailed description of the approach can be found in ref. 78.

1. Davis MB, Shaw RG (2001) Range shifts and adaptive responses to Quaternary climate change. Science 292:673–679. 2. IPCC (2013) Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, eds Stocker TF, et al. (Cambridge Univ Press, Cambridge, UK). 3. Hampe A, Jump AS (2011) Climate relicts: Past, present, future. Annu Rev Ecol Evol Syst 42:313–333.

Sánchez-Salguero et al.

Vulnerability of Mediterranean Firs. To explore the stability of Mediterranean Abies refugia under different forcing climate scenarios, we compared the observed and forecasted growth responses to temperature (gT) and soil moisture (gM) for the 2011–2049 and 2050–2100 periods under the RCP scenarios (Fig. 1). We did this by analyzing the growth projections of each site climate-based model including Extreme Climate Indices (86). All forecasted values exceeding at least two SDs of the observed mean growth response (gT, gM) during the growing season were defined as vulnerable climate refugia; i.e., refugia exposed to a highly elevated climate instability (87). We chose the growth responses to observed extreme climate events to declare vulnerability since drought- or heat-induced forest dieback episodes have already been documented in response to heat waves and dry spells in several Mediterranean fir species, particularly in the case of Silver fir (A. alba) and Spanish fir (A. pinsapo) (13, 20, 36, 40, 83). ACKNOWLEDGMENTS. We thank the Climatic Research Unit for providing the climate databases used in this study; Dr. Geert Jan van Oldenborgh for assistance; the contributors to the International Tree-Ring Data Bank; José B. López-Quintanilla, José L. Sánchez Vallejo, Francisco Jarillo, José Ramón Guzmán Álvarez, and Fernando Ríos (Consejería de Medio Ambiente y Ordenación del Territorio, Junta de Andalucía); Miguel Sánchez González and Pedro Salguero Fernández for their support during fieldwork; and useful comments provided by K. Seftigen. A special thanks for the great efforts made by our late colleague, Fernando Molina Vázquez. R.S.-S. is supported by Postdoctoral Grant IJCI-201525845 (Fondo Europeo de Desarrollo Regional funds). This study was supported by projects 387/2011 (Organismo Autónomo de Parques Nacionales, Spanish Ministry of Environment), CGL2011-26654; FunDiver (CGL2015-69186C2-1-R); CoMo-ReAdapt (CGL2013-48843-C2-1-R) (Spanish Ministry of Economy); and Cambio Climático y Adaptación de los Bosques del Pirineo (CANOPEE) (Interreg V-A POCTEFA 2014-2020-Fondo Europeo de Desarrollo Regional funds).

4. Allen CD, Breshears DD, McDowell NG (2015) On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene. Ecosphere 6:1–55. 5. Anderegg WRL, Kane J, Anderegg LDL (2013) Consequences of widespread tree mortality triggered by drought and temperature stress. Nat Clim Chang 3:30–36. 6. Nicotra AB, et al. (2010) Plant phenotypic plasticity in a changing climate. Trends Plant Sci 15:684–692.

PNAS Early Edition | 7 of 9

PNAS PLUS ECOLOGY

Climate Projections. Only those climate variables highly correlated with TRWi (r > j0.30j, P < 0.05) were considered in the climate-based models and TRWi projections under different climate scenarios. The climate data projected for the 21st century were downloaded and downscaled (cf. 72) at a 0.5° spatial resolution from the fifth phase of the Coupled Model Intercomparison Project (ensemble CMIP5) (74). We used data for the scenario (RCP 8.5) that most closely tracked recent historical emissions (75), and one lower-emission scenario (RCP 4.5) in which the increase in annual emissions is more gradual during the early 21st century and declines after the mid-21st century. These scenarios result in 1.7–4.8 °C global warming by the year 2100, relative to the late–20th-century baseline. The RCP 4.5 scenario represents a situation where radiative forcing peaks at 4.5 W·m−2 after 2100 with temperature increases ranging between 0.9 and 2.6 °C during the 21st century, whereas in the RCP 8.5 scenario radiative forcing continuously rises to reach 8.5 W·m−2 in 2100 with a warming increase ranging between 1.4 and 4.8 °C (26).

Temperature (Ti) and soil moisture (Mi) growth parameters were sampled uniformly across intervals, and the growth parameter set producing the simulation that correlated most significantly with the corresponding observed TRWi series for each site was then used in the simulations. In addition, other parameters (e.g., soil moisture, runoff, root depth) were taken from published studies (25, 76–82). The model was evaluated 10,000 times for each site using three parallel Markov Chain Monte Carlo chains with uniform prior distribution for each parameter and a white Gaussian noise model error (78). To compute annual TRWi values, we integrated the overall simulated growth rates (i.e., the point-wise minimum of monthly gT, gM, and gE) over the time window from September of the year before growth to October of the year of tree-ring formation. This period was determined following previous xylogenesis and dendroecological studies performed on these species (15, 41, 50, 83). To evaluate the temporal stability of the calibrated growth response functions, we divided the period 1950–2010 into two 30-y intervals (1950–1980, 1980–2010) and withheld the second half for validation of the parameters estimated in the first half. Climate-growth relationships were reexamined by applying stepwise multiple linear regressions. This allowed the effects of climate and extreme indices in the observed TRWi data to be identified and to project TRWi through individual-site climate-based equations under future scenarios (Fig. 1). All continuous predictor variables were standardized to give them the same weight in the fitted models (i.e., the mean was subtracted from each value and divided by the SD), enabling the interactions to be tested and compared (84). In addition, we evaluated the existence of multicollinearity among explanatory variables by calculating the variance inflation factor, which, being lower than 2, confirmed no redundancy problems with the data. We used the function “step” of the R package “stats” (62) using the lowest Akaike Information Criterion to select the final regression equations. The models were fitted using Generalized Least-Squares estimation and the R package “nlme” (85). The selected models were run to forecast the TRWi of each site (hereafter, TRWip) for the 2011–2049 and 2050–2100 periods under the two selected RCP scenarios. Finally, we ran VS-Lite models on the TRWip series over the same periods to forecast growth responses (i.e., gT, gM) and growth function parameters (T1, T2 and M1, M2) under future climate projections (see above).

ENVIRONMENTAL SCIENCES

whether climatic factors were responsible for conspicuously narrower or wider tree rings as described in ref. 69 (Fig. 1). Species- and site-specific pointer years were calculated on raw annual ring-width values for each tree and were then transformed into normalized values (Cropper index). This method z-transforms tree growth in year i within a symmetric moving window of n years (5 y in our case), thereby providing the number of SDs of tree growth in individual years from the window average (70). To calculate the intensity of interannual growth anomalies, we defined the probability density functions of the standardized normal distribution of data (see ref. 71) and identified the strong negative growth anomalies during extreme events. The resulting time series of normalized values allowed for the interpretation of the growth characteristics under extreme events in terms of yearly SD units for each tree, site, and species. Cropper values were then normalized to a mean of 0 and an SD of 1. To evaluate the probability of a strong negative growth anomaly, threshold values were defined as follows: “extreme”: > 1.64; “strong”: > 1.28; and “weak” > 1; values between −1 and +1 were considered as “normal” (see ref. 71). The weak threshold value was used to maximize the number of years for further analyses. When more than 50% of all trees within a chronology exceeded the defined threshold, the year was considered a pointer year (69). To assess the impact on growth (sensitivity) of the occurrence of pointer years, we also calculated the annual percentage of growth change (GC) of each tree by calculating the ratio of tree growth in year i and the average growth of the 5 preceding years. To assess the tree response to 21st-century impact on forecasted growth trends, the TRWi reduction (GR: the ratio of TRWi in year i and the average TRWi of 5 preceding years) was also calculated for each tree. To calculate the pointer years and the GC components, we used the pointRes package (72).

7. Tzedakis PC, Lawson IT, Frogley MR, Hewitt GM, Preece RC (2002) Buffered tree population changes in a quaternary refugium: Evolutionary implications. Science 297: 2044–2047. 8. Linares JC, Delgado-Huertas A, Carreira JA (2011) Climatic trends and different drought adaptive capacity and vulnerability in a mixed Abies pinsapo - Pinus halepensis forest. Clim Change 105:67–90. 9. Seneviratne SI, Donat MG, Pitman AJ, Knutti R, Wilby RL (2016) Allowable CO2 emissions based on regional and impact-related climate targets. Nature 529:477–483. 10. Hampe A, Rodríguez-Sánchez F, Dobrowski S, Hu FS, Gavin DG (2013) Climate refugia: From the last glacial maximum to the twenty-first century. New Phytol 197:16–18. 11. Linares JC (2011) Biogeography and evolution of Abies (Pinaceae) in the Mediterranean Basin. The roles of long-term climatic changes and glacial refugia. J Biogeogr 38: 619–630. 12. Linares JC, Carreira JA (2009) Temperate-like stand dynamics in relict Mediterraneanfir (Abies pinsapo Boiss.) forests from Southern Spain. Ann For Sci 66:610. 13. Linares JC, Carreira JA, Ochoa V (2011) Human impacts drive forest structure and diversity. Insights from Mediterranean mountain forest dominated by Abies pinsapo. Eur J For Res 130:533–542. 14. Tinner W, et al. (2013) The past ecology of Abies alba provides new perspectives on future responses of silver fir forests to global warming. Ecol Monogr 83:419–439. 15. Aussenac G (2002) Ecology and ecophysiology of circum-Mediterranean firs in the context of climate change. Ann For Sci 59:823–832. 16. Bigler C, Bugmann H (2004) Predicting the time of tree death using dendrochronological data. Ecol Appl 14:902–914. 17. Bigler C, Bräker OU, Bugmann H, Dobbertin M, Rigling A (2006) Drought as an inciting mortality factor in Scots pine stands of the Valais, Switzerland. Ecosystems (N Y) 9:330–343. 18. Williams A, et al. (2013) Temperature as a potent driver of regional forest drought stress and tree mortality. Nat Clim Chang 3:292–297. 19. Camarero JJ, Gazol A, Sangüesa-Barreda G, Oliva J, Vicente-Serrano SM (2015) To die or not to die: Early-warning signals of dieback in response to a severe drought. J Ecol 103:44–57. 20. Cailleret M, et al. (2017) A synthesis of radial growth patterns preceding tree mortality. Glob Change Biol 23:1675–1690. 21. García-Valdés R, Zavala MA, Araújo MB, Purves DW (2013) Chasing a moving target: Projecting climate change-induced shifts in non-equilibrial tree species distributions. J Ecol 101:441–453. 22. García-Valdés R, Gotelli NJ, Zavala MA, Purves DW, Araújo MB (2015) Effects of climate, species interactions, and dispersal on decadal colonization and extinction rates of Iberian tree species. Ecol Modell 309:118–127. 23. Talluto MV, Boulangeat I, Vissault S, Thuiller W, Gravel D (2017) Extinction debt and colonization credit delay range shifts of eastern North American trees. Nat Ecol Evol 1:0182. 24. Zhang J, et al. (2017) Extinction risk of North American seed plants elevated by climate and land-use. J Appl Ecol 54:303–312. 25. Sánchez-Salguero R, et al. (2017) Assessing forest vulnerability to climate warming using a process-based model of tree growth: Bad prospects for rear-edges. Glob Change Biol 23:2705–2719. 26. van Vuuren DP, et al. (2011) The representative concentration pathways: An overview. Clim Change 109:5–31. 27. Vaganov EA, Hughes MK, Shashkin AV (2006) Growth Dynamics of Conifer Tree Rings (Springer, Berlin). 28. Vaganov EA, Anchukaitis KJ, Evans MN (2011) How well understood are the processes that create dendroclimatic records? A mechanistic model of the climatic control on conifer tree-ring growth dynamics. Dendroclimatology, Developments in Paleoenvironmental Research, eds Hughes MK, Swetnam TW, Diaz HF (Springer, Dordrecht, The Netherlands), pp 37–75. 29. Hampe A, Petit RJ (2005) Conserving biodiversity under climate change: The rear edge matters. Ecol Lett 8:461–467. 30. Dawson TP, Jackson ST, House JI, Prentice IC, Mace GM (2011) Beyond predictions: Biodiversity conservation in a changing climate. Science 332:53–58. 31. Rogers BM, Jantz P, Goetz SJ (2017) Vulnerability of eastern US tree species to climate change. Glob Change Biol 23:3302–3320. 32. Ruosch M, et al. (2016) Past and future evolution of Abies alba forests in Europe comparison of a dynamic vegetation model with palaeo data and observations. Glob Change Biol 22:727–740. 33. Büntgen U, et al. (2014) Placing unprecedented recent fir growth in a European-wide and Holocene-long context. Front Ecol Environ 12:100–106. 34. Gazol A, et al. (2015) Distinct effects of climate warming on populations of silver fir (Abies alba) across Europe. J Biogeogr 42:1150–1162. 35. McDowell N, et al. (2008) Mechanisms of plant survival and mortality during drought: Why do some plants survive while others succumb to drought? New Phytol 178: 719–739. 36. Camarero JJ, Bigler C, Linares JC, Gil-Pelegrin E (2011) Synergistic effects of past historical logging and drought on the decline of Pyrenean silver fir forests. For Ecol Manage 262:759–769. 37. Anderegg WRL, et al. (2012) The roles of hydraulic and carbon stress in a widespread climate-induced forest die-off. Proc Natl Acad Sci USA 109:233–237. 38. Koutavas A (2013) CO2 fertilization and enhanced drought resistance in Greek firs from Cephalonia Island, Greece. Glob Change Biol 19:529–539. 39. Peñuelas J, Hunt JM, Ogaya R, Jump AS (2008) Twentieth century changes of tree-ring delta C-13 at the southern range-edge of Fagus sylvatica: Increasing water-use efficiency does not avoid the growth decline induced by warming at low altitudes. Glob Change Biol 14:1076–1088.

8 of 9 | www.pnas.org/cgi/doi/10.1073/pnas.1708109114

40. Linares JC, Camarero JJ (2012) From pattern to process: Linking intrinsic water-use efficiency to drought-induced forest decline. Glob Change Biol 18:1000–1015. 41. Papadopoulos AM (2016) Tree-ring patterns and climate response of Mediterranean fir populations in Central Greece. Dendrochronologia 40:17–25. 42. Rolland C (1993) Tree-ring and climate relationships for Abies alba in the internal Alps. Tree-Ring Bull 53:1–11. 43. Adams HD, et al. (2009) Temperature sensitivity of drought-induced tree mortality portends increased regional die-off under global-change-type drought. Proc Natl Acad Sci USA 106:7063–7066. 44. George JP, et al. (2015) Inter- and intra-specific variation in drought sensitivity in Abies spec. and its relation to wood density and growth traits. Agric Meteorol 214215:430–443. 45. Valladares F, et al. (2014) The effects of phenotypic plasticity and local adaptation on forecasts of species range shifts under climate change. Ecol Lett 17:1351–1364. 46. Peguero-Pina JJ, et al. (2011) Hydraulic traits are associated with the distribution range of two closely related Mediterranean firs, Abies alba Mill. and Abies pinsapo Boiss. Tree Physiol 31:1067–1075. 47. Sánchez-Salguero R, et al. (2015) Regulation of water use in the southernmost European fir (Abies pinsapo Boiss.): Drought avoidance matters. Forests 6:2241–2260. 48. Pasho E, Toromani E, Allan AQ (2014) Climatic impact on tree-ring widths in Abies borisii-regis forests from South-East Albania. Dendrochronologia 32:237–244. 49. Nourtier M, et al. (2014) Transpiration of silver Fir (Abies alba Mill.) during and after drought in relation to soil properties in a Mediterranean mountain area. Ann Sci 71: 683–695. 50. Cuny HE, Rathgeber CBK, Lebourgeois F, Fortin M, Fournier M (2012) Life strategies in intra-annual dynamics of wood formation: Example of three conifer species in a temperate forest in north-east France. Tree Physiol 32:612–625. 51. Cheddadi R, et al. (2016) Temperature range shifts for three European tree species over the last 10,000 years. Front Plant Sci 7:1581. 52. Sánchez-Robles JM, et al. (2014) Phylogeography of SW Mediterranean firs: Different European origins for the North African Abies species. Mol Phylogenet Evol 79:42–53. 53. Pecl GT, et al. (2017) Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being. Science 355:1389. 54. Vitasse Y, et al. (2011) Assessing the effects of climate change on the phenology of European temperate trees. Agric Meteorol 151:969–980. 55. Matías L, Jump AS (2015) Asymmetric changes of growth and reproductive investment herald altitudinal and latitudinal range shifts of two woody species. Glob Change Biol 21:882–896. 56. Millar CI, Stephenson NL (2015) Temperate forest health in an era of emerging megadisturbance. Science 349:823–826. 57. Caudullo G, Tinner W (2016) Abies–Circum-Mediterranean firs in Europe: Distribution, habitat, usage and threats. European Atlas of Forest Tree Species, eds San-MiguelAyanz J, de Rigo D, Caudullo G, Houston Durrant T, Mauri A (Publications Office of the European Union, Luxembourg). 58. Carrer M, Nola P, Motta R, Urbinati R (2010) Contrasting tree-ring growth to climate responses of Abies alba toward the southern limit of its distribution area. Oikos 119: 1515–1525. 59. Grissino-Mayer HD (2001) Evaluating crossdating accuracy: A manual and tutorial for the computer program COFECHA. Tree-Ring Res 57:205–221. 60. Harris I, Jones PD, Osborn TJ, Lister DH (2014) Updated high-resolution grids of monthly climatic observations: The CRU TS3.10 Dataset. Int J Climatol 34:623–642. 61. Wigley TML, Briffa KR, Jones PD (1984) On the average value of correlated time series, with applications in dendroclimatology and hydrometeorology. J Clim 23:201–213. 62. R Development Core Team (2017) R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna), Version 3.3.3. Available at www.r-project.org/. Accessed March 30, 2017. 63. Bunn AG (2010) Statistical and visual crossdating in R using the dplR library. Dendrochronologia 28:251–258. 64. Frank D, et al. (2015) Effects of climate extremes on the terrestrial carbon cycle: Concepts, processes and potential future impacts. Glob Change Biol 21:2861–2880. 65. Donat MG, et al. (2013) Updated analyses of temperature and precipitation extreme indices since the beginning of the twentieth century: The HadEX2 dataset. J Geophys Res Atmos 118:2098–2118. 66. Sillmann J, Kharin VV, Zhang X, Zwiers FW, Bronaugh D (2013) Climate extremes indices in the CMIP5 multimodel ensemble: I. Model evaluation in the present climate. J Geophys Res 118:1716–1733. 67. Sillmann J, Kharin VV, Zwiers FW, Zhang X, Bronaugh D (2013) Climate extremes indices in the CMIP5 multimodel ensemble: II. Future climate projections. J Geophys Res 118:2473–2493. 68. Cook ER, Kairiukstis LA (1990) Methods of Dendrochronology: Applications in the Environmental Sciences (Kluwer Academic Publishers, International Institute for Applied Systems Analysis, Dordrecht, The Netherlands). 69. Schweingruber FH, Eckstein D, Bachet S, Bräker OU (1990) Identification, presentation and interpretation of event years and pointer years in dendrochronology. Dendrochronologia 8:9–38. 70. Cropper JP (1979) Tree ring skeleton plotting by computer. Tree-Ring Bull 39:47–60. 71. Neuwirth B, Schweingruber FH, Winiger M (2007) Spatial patterns of central European pointer years from 1901 to 1971. Dendrochronologia 24:79–89. 72. van der Maaten-Theunissen M, van der Maaten E, Bouriaud O (2015) pointRes: An R package to analyze pointer years and components of resilience. Dendrochronologia 35:34–38. 73. Moreno A, Hasenauer H (2016) Spatial downscaling of European climate data. Int J Climatol 36:1444–1458. 74. Taylor KE, Stouffer RJ, Meehl GA (2012) An overview of CMIP5 and the experiment design. Bull Am Meteorol Soc 93:485–498.

Sánchez-Salguero et al.

PNAS PLUS

81. Touchan R, Shishov VV, Meko DM, Nouiri I, Grachev A (2012) Process based model sheds light on climate sensitivity of Mediterranean tree-ring width. Biogeosciences 9:965–972. 82. Mina M, Martin-Benito D, Bugmann H, Cailleret M (2016) Forward modeling of tree-ring width improves simulation of forest growth responses to drought. Agric Meteorol 221:13–33. 83. Linares JC, Camarero JJ, Carreira JA (2009) Plastic responses of Abies pinsapo xylogenesis to drought and competition. Tree Physiol 29:1525–1536. 84. Zuur AF, Ieno EN, Walker N, Saveliev AA, Smith GM (2009) Mixed Effects Models and Extensions in Ecology with R (Springer, New York). 85. Pinheiro J, Bates D, DebRoy S, Sarkar D (2009) nlme: Linear and Nonlinear Mixed Effects Models. R package Version 3.1-96. Available at cran.r-project.org/web/packages/ nlme/index.html. Accessed April 15, 2017. 86. Sass-Klaassen U, et al. (2016) A tree-centered approach to assess impacts of extreme climatic events on forests. Front Plant Sci 7:1069. 87. Cavin L, Jump AS (2017) Highest drought sensitivity and lowest resistance to growth suppression are found in the range core of the tree Fagus sylvatica L. not the equatorial range edge. Glob Change Biol 23:362–379.

ENVIRONMENTAL SCIENCES

ECOLOGY

75. Lelieveld J, et al. (2016) Strongly increasing heat extremes in the Middle East and North Africa (MENA) in the 21st century. Clim Change 137:245–260. 76. Tolwinski-Ward SE, Tingley MP, Evans MN, Hughes MK, Nychka DW (2015) Probabilistic reconstructions of local temperature and soil moisture from tree-ring data with potentially time-varying climatic response. Clim Dyn 44:791–806. 77. Tolwinski-Ward SE, Evans MN, Hughes MK, Anchukaitis KJ (2011) An efficient forward model of the climate controls on interannual variation in tree-ring width. Clim Dyn 36:2419–2439. 78. Tolwinski-Ward SE, Anchukaitis KJ, Evans MN (2013) Bayesian parameter estimation and interpretation for an intermediate model of tree-ring width. Clim Past 9: 1481–1493. 79. Huang J, van den Dool HM, Georgakakos KP (1996) Analysis of model-calculated soil moisture over the United States (1931–1993) and applications to long-range temperature forecasts. J Clim 9:1350–1362. 80. Evans MN, et al. (2006) A forward modeling approach to paleoclimatic interpretation of tree-ring data. J Geophys Res 111:G03008.

Sánchez-Salguero et al.

PNAS Early Edition | 9 of 9

SUPPORTING INFORMATION (SI Appendix)

Climate extremes and predicted warming threaten Mediterranean Holocene firs forests refugia Raúl Sánchez-Salguero, J. Julio Camarero, Marco Carrer, Emilia Gutiérrez, Arben Q. Alla, Laia Andreu-Hayles, Andrea Hevia, Athanasios Koutavas, Elisabet Martínez-Sancho, Paola Nola, Andreas Papadopoulos, Edmond Pasho, Ervin Toromani, José A. Carreira and Juan C. Linares

The purpose of this supplementary material is to provide information that is of less central importance to the paper and that cannot be included in the main body of the text because of space limitations. The Supplementary Information contains 9 tables and 9 figures. The sequence of the supplementary figures and tables follows the citation order in the main text.

1

Supplementary Tables Table S1. Main characteristics of sampled sites (codes are as in SI Appendix, Fig. S1). Dbh: diameter at breast height. Number of trees and radii per site. Tree-ring width (TRW) was determined from cores taken at 1.3 m. Rbt: is the mean correlation between trees of indexed tree-ring width series; AC: is the first order autocorrelation, and MS: is the mean sensitivity considering the common 1930-2010 period. The EPS (expressed population signal) was > 0.85 for all sites in the 1950-2010 period (61).

Tree species

A. alba

A. borisii-regis

A. cephalonica

A. cilicica

A. pinsapo A. maroccana A. tazaotana

Site code

Site name

Country

Latitude (N)

Longitude (-W, +E)

Elevation (m a.s.l.)

Aspect

Slope

Dbh (cm)

Nº trees /radii

TRW ± SD (mm)

Rbt

AC

MS

Time span (years)

ROSE PESC IANN CUGN PANT GARA GUAR PORH PORL SJPE BOUM MONT ABRH ABRM ABRL KARP1 PERT AINO KARP2 PARN leb2* syr2* tu24* DEAC IBAC APSB APSG APSN ABMC ABTZ

Rosello Pescopennataro Iannace (Pollino) Cugno dell’Acero (Pollino) Pantano (Pollino) Gariglione Guara Peña Oroel-High Peña Oroel-Low San Juan de la Peña Boumort Montseny Bredhi i Hotovës High Bredhi i Hotovës Middle Bredhi i Hotovës Low Karpenissi (Evrytania) Pertouli (Trikala) Ainos Mountain Karpenissi (Evrytania) Parnitha (Attica) Wadi Balat Bedayat Al Khandak Al Tawil Gazipasa Forest, Antalya Derebukak Ibradi Civarinda Sierra Bermeja Sierra de Grazalema Sierra de las Nieves Talassemtane Tazaot

Italy Italy Italy Italy Italy Italy Spain Spain Spain Spain Spain Spain Albania Albania Albania Greece Greece Greece Greece Greece Lebanon Syria Turkey Turkey Turkey Spain Spain Spain Morocco Morocco

42.50 42.48 40.57 40.57 40.57 39.20 42.30 42.52 42.52 42.52 42.20 41.77 40.34 40.35 40.35 38.88 39.55 38.14 38.87 38.17 34.47 35.57 36.45 37.32 37.21 36.48 36.77 36.72 35.13 35.23

14.57 14.50 16.35 16.37 16.40 17.07 -0.12 -0.32 -0.32 -0.41 1.20 2.43 20.38 20.39 20.38 21.86 21.47 20.67 21.86 23.74 36.23 36.20 32.52 31.46 31.49 -5.21 -5.42 -5.10 -5.12 -5.08

1086 1370 1400 1350 1500 1680 1428 1604 1587 1393 1583 1550 1361 1144 1058 1120 1270 1450 1260 1223 1175 1450 1770 1600 1500 1300 1275 1350 1600 1685

E E N-W N-W N-E S-E N-NW N-NW N N-NE NE SW NW SE N E E-SE NE NE E NW W SW SO SO N N N-NE NO NE

35 28 20 25 25 28 47 36 24 22 35 22 25 22 19 48 23 35 35 30 40 33 27 27 31

89.7 67.1 92.0 76.6 111.8 83.4 52.5 59.1 46.5 46.0 64.2 68.0 46.3 45.4 39.1 47.0 49.0 68.4 45.3 44.0 34.3 23.9 49.1 30.1 28.5 37.3 32.1 31.8 29.5 45.1

15/32 15/35 14/31 16/35 16/34 13/27 23/46 12/24 23/46 14/28 17/40 21/50 22/31 36/54 18/28 12/23 12/21 23/54 15/25 12/19 16/19 33/35 7/7 6/10 5/10 23/26 64/85 69/97 14/23 14/18

2.45 ± 0.55 1.61 ± 0.50 1.38 ± 0.95 1.29 ± 0.59 1.39 ± 0.80 1.59 ± 0.82 3.11 ± 0.79 2.57 ± 0.79 2.41 ± 0.70 2.37 ± 0.49 1.79 ± 1.24 1.22 ± 0.78 2.30 ± 1.03 2.43 ± 1.26 1.83 ± 0.68 1.85 ± 0.85 2.46 ± 0.34 1.07 ± 0.66 1.47 ± 0.37 1.46 ± 0.24 1.15 ± 0.33 0.79 ± 0.26 1.59 ± 0.41 2.48 ± 0.72 2.28 ± 0.61 1.17 ± 0.43 2.25 ± 0.59 1.65 ± 0.97 1.45 ± 0.19 2.91 ± 0.73

0.30 0.42 0.21 0.34 0.29 0.31 0.48 0.71 0.67 0.57 0.36 0.37 0.26 0.3 0.31 0.61 0.57 0.54 0.44 0.49 0.46 0.52 0.50 0.59 0.56 0.57 0.50 0.52 0.42 0.48

0.09 0.32 0.13 0.08 0.22 0.23 0.09 0.78 0.76 0.85 0.16 0.26 0.77 0.76 0.72 0.73 0.75 0.68 0.80 0.77 0.82 0.77 0.83 0.79 0.62 0.76 0.70 0.81 0.73 0.82

0.15 0.16 0.18 0.17 0.14 0.15 0.20 0.20 0.18 0.23 0.22 0.17 0.22 0.26 0.2 0.19 0.19 0.23 0.19 0.21 0.25 0.32 0.22 0.23 0.22 0.30 0.26 0.24 0.22 0.23

1848–1998 1823–1998 1720–1998 1703–1998 1697–1998 1719–1999 1887–2011 1892–2000 1889–2000 1855–1999 1804–2002 1587–2010 1900–2010 1874–2010 1893–2010 1890–2012 1880-2013 1769-2007 1857–2013 1868–2012 1722-2001 1795-2001 1797-2003 1869-2009 1943-2009 1795-2014 1764-2015 1710-2015 1798-2006 1689-2009

* ITRDB database chronologies (https://www.ncdc.noaa.gov/data-access/paleoclimatology-data)

2

Table S2. Modelled percentage change in average tree-ring width (TRWi) for projected climate change trends (TRWip) considering the ensemble IPCC AR5 emission scenarios (RCP 4.5 and RCP 8.5). Values are calculated such as the projected 2011-2049 and 2050-2100 annual mean minus measured 1950–2010 annual mean. Trends (τ, Pearson correlation coefficient) of site mean tree-ring width indices series considering the 2011–2100 period following growth projections based on two IPCC AR5 emission scenarios (see Fig. 4 and SI Appendix Fig. S8). 2011-2049 RCP 4.5 Tree species

A. alba

A. borisii-regis

A. cephalonica

A. cilicica

A. pinsapo A. maroccana A. tazaotana

Site code ROSE PESC IANN CUGN PANT GARA GUAR PORH PORL SJPE BOUM MONT ABRH ABRM ABRL KARP1 PERT AINO KARP2 PARN leb2 syr2 tu24 DEAC IBAC APSB APSG APSN ABMC ABTZ

Mean 2.4 7.1 -5.6 -9.1 -2.3 -3.2 -3.9 -18.8 -12.1 -31.6 3.7 45.9 18.9 -9.3 -10.9 -2.5 -9.8 -1.3 -7.8 13.9 -9.3 -7.8 12.5 30.6 4.5 15.2 2.2 8.2 -3.6 -9.6

RCP 8.5 Mean 42.6 37.4 -7.8 -13.3 -10.6 -4.9 -5.8 -8.3 -10.4 -5.8 -2.6 35.1 15.2 -8.7 -6.8 -3.2 -7.8 3.2 -7.5 7.4 -11.3 6.3 15.1 20.4 6.8 14.5 15.1 14.2 -1.6 -10.1

2050-2100 RCP 4.5 Mean -5.2 -7.2 -8.7 -17.9 -11.2 -6.5 -14.4 -35.2 -18.8 -68.2 34.4 65.4 27.1 -14.6 -38.1 -6.8 -48.3 -10.8 -24.7 4.8 -36.9 -16.1 32.1 53.1 6.7 68.3 -5.7 17.3 -12.2 -25.7

RCP 8.5 Mean 68.1 61.0 -19.3 -35.6 -13.3 -7.5 -18.7 -42.3 -32.7 -59.4 35.7 59.8 30.1 -41.8 -21.6 -8.7 -58.6 -11.5 -33.3 5.2 -38.2 -21.3 35.4 49.8 6.9 47.3 -10.6 14.7 -16.5 -30.0

2011-2100 Trend (τ Pearson) RCP 4.5 -0.15 -0.34 -0.25 -0.23 -0.44 0.12 -0.33 -0.73 -0.50 -0.81 0.78 0.91 0.62 -0.62 -0.83 -0.33 -0.77 -0.27 -0.69 -0.30 -0.79 -0.29 0.55 0.79 0.11 0.68 -0.41 0.17 -0.53 -0.64

RCP 8.5 0.85 0.58 -0.79 -0.71 -0.63 -0.26 -0.77 -0.93 -0.83 -0.92 0.85 0.91 0.85 -0.97 -0.87 -0.51 -0.92 -0.58 -0.87 0.10 -0.90 -0.75 0.80 0.88 0.07 0.79 -0.77 -0.10 -0.84 -0.86

3

Table S3. Pearson correlation coefficients calculated between observed site series of tree-ring width indices (TRWi) and VS-Lite indices (TRWiVSL) for the calibration period (1950-2010); and for the sub-periods 1950–1980 and 1980–2010 (see SI Appendix, Fig. S6), also considering projected site series of tree-ring width indices (TRWip) and VS-Lite future indices (TRWip-VSL) for the periods 2010-2050 and 20502100 and two IPCC AR5 emission scenarios (RCP 4.5 and RCP 8.5). Correlation values higher than 0.25 are significant at P < 0.05.

IPCC emission scenarios Tree species

A. alba

A. borisii-regis

A. cephalonica

A. cilicica

A. pinsapo A. maroccana A. tazaotana

Site code ROSE PESC IANN CUGN PANT GARA GUAR PORH PORL SJPE BOUM MONT ABRH ABRM ABRL KARP1 PERT AINO KARP2 PARN leb2 syr2 tu24 DEAC IBAC APSB APSG APSN ABMC ABTZ

1950-1980

Period 1980-2010

1950-2010

0.41 0.33 0.45 0.47 0.35 0.41 0.59 0.45 0.49 0.52 0.41 0.55 0.40 0.43 0.41 0.36 0.33 0.18 0.31 0.28 0.55 0.20 0.37 0.48 0.45 0.37 0.45 0.38 0.49 0.47

0.33 0.22 0.51 0.30 0.24 0.48 0.62 0.32 0.38 0.38 0.53 0.41 0.32 0.35 0.34 0.37 0.42 0.15 0.32 0.38 0.47 0.15 0.27 0.49 0.31 0.37 0.47 0.52 0.35 0.41

0.36 0.28 0.42 0.31 0.26 0.39 0.67 0.34 0.44 0.44 0.40 0.44 0.31 0.41 0.39 0.38 0.28 0.26 0.36 0.24 0.51 0.19 0.30 0.43 0.36 0.37 0.44 0.46 0.38 0.41

RCP 4.5 2010-2050 2050-2100 0.34 0.27 0.46 0.35 0.24 0.38 0.47 0.33 0.43 0.44 0.65 0.82 0.41 0.63 0.57 0.29 0.45 0.25 0.24 0.26 0.79 0.29 0.50 0.74 0.28 0.27 0.43 0.29 0.35 0.64

0.52 0.36 0.38 0.34 0.25 0.55 0.40 0.41 0.33 0.50 0.59 0.77 0.35 0.58 0.41 0.43 0.49 0.36 0.56 0.24 0.39 0.29 0.40 0.58 0.37 0.25 0.24 0.26 0.37 0.42

RCP 8.5 2010-2050 2050-2100 0.69 0.29 0.57 0.42 0.34 0.49 0.50 0.53 0.55 0.62 0.58 0.87 0.63 0.61 0.85 0.40 0.58 0.34 0.38 0.31 0.77 0.34 0.63 0.73 0.34 0.36 0.49 0.42 0.71 0.80

0.68 0.30 0.70 0.44 0.54 0.34 0.59 0.85 0.64 0.86 0.65 0.90 0.69 0.79 0.59 0.45 0.89 0.44 0.67 0.33 0.82 0.72 0.74 0.34 0.29 0.51 0.59 0.35 0.82 0.81

4

Table S4. Statistics of the Bayesian estimation of VS-Lite growth response parameters (1950-2010). Statistics of the Bayesian estimation of the site-by-site tuned VS-Lite growth response parameters (T1, T2, M1, and M2) for the period 1950-2010 and two sub-periods (1950-1980 and 19802010).

Tree species

Abies alba-Italy

Abies alba-Spain

Abies borisii-regis

Abies cephalonica

Parameter (unit)

1950-2010

1950-1980

1980-2010

Mean

Min

Max

SD

Mean

Min

Max

SD

Mean

Min

Max

SD

T1 (ºC)

7.893

6.950

8.323

0.495

8.254

8.032

8.638

0.206

7.294

6.758

8.193

0.575

T2 (ºC)

12.839

11.179

13.811

0.995

12.986

11.130

15.292

1.381

12.659

10.453

15.613

2.033

M1 (v/v)

0.024

0.001

0.083

0.030

0.025

0.007

0.068

0.022

0.034

0.002

0.064

0.024

M2 (v/v)

0.163

0.102

0.319

0.085

0.184

0.150

0.210

0.026

0.127

0.100

0.178

0.030

T1 (ºC)

2.573

1.589

5.707

1.564

4.284

2.358

8.246

2.097

3.111

1.696

6.548

1.938

T2 (ºC)

13.843

10.840

23.020

4.693

13.443

10.650

18.449

3.187

12.993

10.455

17.033

2.797

M1 (v/v)

0.057

0.020

0.090

0.031

0.031

0.001

0.071

0.028

0.040

0.026

0.055

0.011

M2 (v/v)

0.375

0.261

0.469

0.091

0.368

0.267

0.482

0.080

0.387

0.260

0.483

0.081

T1 (ºC)

3.932

2.911

6.963

1.711

3.378

2.100

5.715

1.371

5.381

3.074

7.891

2.175

T2 (ºC)

12.225

10.585

14.023

1.461

12.568

11.007

14.420

1.245

11.697

10.917

12.890

0.814

M1 (v/v)

0.036

0.003

0.081

0.034

0.040

0.006

0.088

0.030

0.060

0.025

0.087

0.025

M2 (v/v)

0.389

0.177

0.489

0.122

0.395

0.190

0.487

0.119

0.341

0.164

0.470

0.148

T1 (ºC)

7.340

5.438

8.563

1.669

7.285

5.157

8.570

1.856

8.181

8.058

8.406

0.195

T2 (ºC)

13.014

11.238

14.822

1.792

12.300

10.570

15.488

2.764

13.958

11.000

18.012

3.632

M1 (v/v)

0.035

0.008

0.079

0.038

0.025

0.005

0.048

0.022

0.054

0.053

0.055

0.001

M2 (v/v)

0.135

0.109

0.177

0.037

0.144

0.113

0.199

0.048

0.163

0.106

0.195

0.050

5

Table S4. continued Tree species

Abies cilicica

Abies pinsapo

Abies maroccana

Abies tazaotana

Parameter (unit)

1950-2010

1950-1980

1980-2010

Mean

Min

Max

SD

Mean

Min

Max

SD

Mean

Min

Max

SD

T1 (ºC)

7.183

6.020

8.390

1.068

7.714

5.600

8.600

1.225

6.803

5.600

8.215

1.131

T2 (ºC)

13.209

11.079

17.847

2.699

14.790

10.530

18.680

3.551

15.088

10.566

23.404

5.661

M1 (v/v)

0.033

0.006

0.053

0.021

0.049

0.002

0.093

0.038

0.041

0.003

0.073

0.030

M2 (v/v)

0.296

0.145

0.427

0.121

0.249

0.140

0.356

0.077

0.312

0.102

0.458

0.174

T1 (ºC)

4.724

3.974

5.868

1.006

6.221

3.002

7.922

2.789

3.411

3.163

3.742

0.299

T2 (ºC)

13.136

11.144

14.398

1.746

15.358

12.549

18.562

3.026

13.387

11.525

15.491

1.994

M1 (v/v)

0.055

0.002

0.093

0.047

0.038

0.001

0.075

0.037

0.090

0.086

0.092

0.003

M2 (v/v)

0.153

0.104

0.227

0.065

0.351

0.263

0.408

0.077

0.121

0.102

0.132

0.017

T1 (ºC)

8.005

-

-

-

4.861

-

-

-

8.666

-

-

-

T2 (ºC)

11.506

-

-

-

11.292

-

-

-

17.264

-

-

-

M1 (v/v)

0.014

-

-

-

0.077

-

-

-

0.020

-

-

-

M2 (v/v)

0.216

-

-

-

0.127

-

-

-

0.454

-

-

-

T1 (ºC)

7.320

-

-

-

6.104

-

-

-

7.323

-

-

-

T2 (ºC)

12.394

-

-

-

13.017

-

-

-

13.123

-

-

-

M1 (v/v)

0.007

-

-

-

0.010

-

-

-

0.032

-

-

-

M2 (v/v)

0.136

-

-

-

0.205

-

-

-

0.100

-

-

-

6

Table S5. Climatic characteristics of sampled sites. Climatic characteristics (mean ± SD) and trends (τ, Pearson correlation coefficient) obtained for the studied Abies species considering annual and seasonal mean temperatures (T) and total precipitation (P). Data were obtained using the homogenized and quality-checked CRU TS 3.23 meteorological dataset gridded at a 0.5º spatial resolution for the period 1950–2010 (60).

Mean Annual

Tree species

A. alba

A. borisii-regis

A. cephalonica

A. cilicica

A. pinsapo A. maroccana A. tazaotana

Mean winter

Mean spring

Site code

T (ºC)

Trend (T)

P (mm)

Trend (P)

T (ºC)

Trend (T)

P (mm)

Trend (P)

T (ºC)

Trend (T)

P (mm)

Trend (P)

ROSE PESC IANN CUGN PANT GARA GUAR PORH PORL SJPE BOUM MONT ABRH ABRM ABRL KARP1 PERT AINO KARP2 PARN leb2 syr2 tu24 DEAC IBAC APSB APSG APSN ABMC ABTZ

15.1 ± 0.5 13.3 ± 0.5 13.9 ± 0.5 13.8 ± 0.5 13.8 ± 0.5 17.0 ± 0.5 10.9 ± 0.7 6.6 ± 0.7 6.6 ± 0.7 6.6 ± 0.7 9.4 ± 0.6 13.4 ± 0.5 10.0 ± 0.5 10.0 ± 0.5 10.0 ± 0.5 12.1 ± 0.5 10.7 ± 0.5 15.9 ± 0.5 12.1 ± 0.5 16.2 ± 0.6 11.8 ± 0.7 17.2 ± 0.6 13.0 ± 0.6 11.5 ± 0.6 11.5 ± 0.6 17.7 ± 0.4 15.7 ± 0.5 15.7 ± 0.5 15.1 ± 0.4 15.1 ± 0.4

0.48 0.51 0.51 0.51 0.51 0.57 0.68 0.70 0.70 0.70 0.65 0.53 0.30 0.30 0.30 0.34 0.30 -0.12 0.34 0.40 0.34 0.25 0.14 0.06 0.06 0.43 0.56 0.56 0.37 0.37

719 ± 144 754 ± 145 665 ± 140 665 ± 140 665 ± 140 763 ± 170 722 ± 159 1029 ± 191 1029 ± 191 1029 ± 191 832 ± 171 692 ± 170 1061 ± 182 1061 ± 182 1061 ± 182 770 ± 141 871 ± 153 688 ± 125 770 ± 141 529 ± 118 647 ± 171 712 ± 149 681 ± 145 643 ± 139 643 ± 139 628 ± 231 612 ± 202 612 ± 202 638 ± 204 638 ± 204

-0.11 -0.10 -0.01 -0.01 -0.01 0.24 -0.18 -0.19 -0.19 -0.19 -0.22 -0.19 -0.31 -0.31 -0.31 -0.19 -0.30 -0.10 -0.19 -0.05 -0.20 -0.11 -0.11 0.02 0.02 0.03 -0.04 -0.04 -0.01 -0.01

7.4 ± 0.8 5.7 ± 0.8 6.6 ± 0.8 6.7 ± 0.7 6.7 ± 0.7 10.3 ± 0.7 3.2 ± 1.1 -0.7 ± 1.1 -0.7 ± 1.1 -0.7 ± 1.1 2.3 ± 1.0 7.1 ± 0.9 1.7 ± 0.8 1.7 ± 0.8 1.7 ± 0.8 4.1 ± 0.7 2.1 ± 0.9 8.9 ± 0.8 4.1 ± 0.7 8.5 ± 0.9 2.9 ± 1.2 8.1 ± 1.2 4.0 ± 1.0 2.1 ± 1.2 2.1 ± 1.2 12.6 ± 0.7 9.5 ± 0.7 9.5 ± 0.7 9.5 ± 0.7 9.5 ± 0.7

0.10 0.15 0.10 0.10 0.10 0.16 0.37 0.37 0.37 0.37 0.32 0.24 -0.02 -0.02 -0.02 -0.07 -0.05 -0.44 -0.07 -0.03 0.02 -0.02 -0.12 -0.19 -0.19 0.39 0.41 0.41 0.23 0.23

211 ± 76 212 ± 77 207 ± 66 207 ± 66 207 ± 66 282 ± 101 160 ± 77 247 ± 102 247 ± 102 247 ± 102 166 ± 89 133 ± 76 379 ± 126 379 ± 126 379 ± 126 292 ± 96 299 ± 102 283 ± 84 292 ± 96 224 ± 86 367 ± 123 369 ± 114 363 ± 126 314 ± 111 314 ± 111 283 ± 163 257 ± 149 257 ± 149 276 ± 142 276 ± 142

-0.18 -0.17 -0.03 -0.03 -0.03 0.11 -0.06 -0.10 -0.10 -0.10 -0.08 0.06 -0.32 -0.32 -0.32 -0.17 -0.22 -0.06 -0.17 0.01 -0.25 -0.18 -0.09 -0.12 -0.12 0.04 -0.12 -0.12 -0.06 -0.06

13.4 ± 0.8 11.6 ± 0.7 11.9 ± 0.8 11.9 ± 0.7 11.9 ± 0.7 14.8 ± 0.7 9.5 ± 0.9 4.9 ± 0.9 4.9 ± 0.9 4.9 ± 0.9 7.8 ± 0.8 11.4 ± 0.8 8.7 ± 0.7 8.7 ± 0.7 8.7 ± 0.7 10.4 ± 0.8 9.4 ± 0.8 13.8 ± 0.7 10.4 ± 0.8 14.0 ± 0.9 10.3 ± 0.9 15.7 ± 0.8 11.4 ± 0.8 10.1 ± 0.9 10.1 ± 0.9 15.9 ± 0.7 13.9 ± 0.7 13.9 ± 0.7 13.1 ± 0.7 13.1 ± 0.7

0.36 0.37 0.41 0.41 0.41 0.49 0.41 0.46 0.46 0.46 0.40 0.28 0.33 0.33 0.33 0.34 0.34 0.01 0.34 0.36 0.22 0.22 0.13 0.06 0.06 0.22 0.34 0.34 0.21 0.21

174 ± 64 173 ± 65 153 ± 54 153 ± 54 153 ± 54 157 ± 64 210 ± 84 292 ± 103 292 ± 103 292 ± 103 233 ± 69 189 ± 81 242 ± 72 242 ± 72 242 ± 72 175 ± 53 207 ± 59 143 ± 43 175 ± 53 113 ± 47 160 ± 65 178 ± 72 155 ± 61 163 ± 58 163 ± 58 157 ± 75 165 ± 72 165 ± 72 170 ± 76 170 ± 76

0.12 0.11 0.05 0.05 0.05 0.24 -0.04 0.01 0.01 0.01 -0.01 -0.08 -0.14 -0.14 -0.14 -0.09 -0.09 0.02 -0.09 -0.07 -0.16 -0.10 -0.20 0.01 0.01 -0.17 -0.29 -0.29 -0.18 -0.18

7

Table S5. Continued

Mean summer

Tree species

A. alba

A. borisii-regis

A. cephalonica

A. cilicica

A. pinsapo A. maroccana A. tazaotana

Site code ROSE PESC IANN CUGN PANT GARA GUAR PORH PORL SJPE BOUM MONT ABRH ABRM ABRL KARP1 PERT AINO KARP2 PARN leb2 syr2 tu24 DEAC IBAC APSB APSG APSN ABMC ABTZ

Mean autumn

T (ºC)

Trend (T)

P (mm)

Trend (P)

T (ºC)

Trend (T)

P (mm)

Trend (P)

23.2 ± 0.9 21.4 ± 0.9 21.8 ± 0.9 21.8 ± 0.9 21.8 ± 0.9 24.5 ± 0.8 19.4 ± 1.0 14.7 ± 1.1 14.7 ± 1.1 14.7 ± 1.1 17.3 ± 1.0 20.7 ± 0.9 18.7 ± 0.9 18.7 ± 0.9 18.7 ± 0.9 20.8 ± 0.9 19.7 ± 0.9 23.8 ± 0.8 20.8 ± 0.9 24.9 ± 1.0 19.8 ± 0.8 25.5 ± 0.7 22.2 ± 0.7 21.0 ± 0.7 21.0 ± 0.7 23.2 ± 0.6 22.7 ± 0.7 22.7 ± 0.7 21.5 ± 0.6 21.5 ± 0.6

0.43 0.45 0.48 0.48 0.48 0.51 0.66 0.67 0.67 0.67 0.57 0.45 0.30 0.30 0.30 0.43 0.33 0.35 0.43 0.51 0.46 0.38 0.37 0.37 0.37 0.40 0.53 0.53 0.35 0.35

113 ± 55 126 ± 61 82 ± 42 82 ± 42 82 ± 42 50 ± 34 159 ± 67 227 ± 86 227 ± 86 227 ± 86 209 ± 78 161 ± 66 119 ± 72 119 ± 72 119 ± 72 66 ± 39 106 ± 56 36 ± 26 66 ± 39 30 ± 24 4±8 18 ± 22 29 ± 24 41 ± 27 41 ± 27 17 ± 16 26 ± 22 26 ± 22 27 ± 23 27 ± 23

-0.01 -0.03 0.04 0.04 0.04 0.05 -0.28 -0.32 -0.32 -0.32 -0.34 -0.23 0.14 0.14 0.14 0.03 0.05 0.20 0.03 -0.03 0.02 0.06 0.15 0.03 0.03 -0.10 -0.19 -0.19 -0.04 -0.04

16.2 ± 0.8 14.5 ± 0.8 15.1 ± 0.7 15.7 ± 0.8 15.7 ± 0.8 18.4 ± 0.7 11.8 ± 0.8 7.7 ± 0.8 7.7 ± 0.8 7.7 ± 0.8 10.4 ± 0.8 14.5 ± 0.7 11.0 ± 0.8 11.0 ± 0.8 11.0 ± 0.8 13.2 ± 0.7 11.6 ± 0.8 17.3 ± 0.7 13.2 ± 0.7 17.3 ± 0.9 14.1 ± 0.8 19.4 ± 0.9 14.6 ± 0.9 12.9 ± 0.9 12.9 ± 0.9 18.9 ± 0.7 16.9 ± 0.8 16.9 ± 0.8 16.4 ± 0.7 16.4 ± 0.7

0.22 0.25 0.23 0.23 0.23 0.30 0.47 0.50 0.50 0.50 0.46 0.38 0.09 0.09 0.09 0.14 0.10 -0.26 0.14 0.20 0.31 0.19 0.08 0.02 0.02 0.15 0.24 0.24 0.11 0.11

221 ± 73 242 ± 86 224 ± 83 224 ± 83 224 ± 83 274 ± 121 192 ± 91 263 ± 108 263 ± 108 263 ± 108 224 ± 80 209 ± 86 320 ± 96 320 ± 96 320 ± 96 238 ± 81 258 ± 83 226 ± 77 238 ± 81 172 ± 78 117 ± 56 147 ± 57 133 ± 68 127 ± 63 127 ± 63 169 ± 95 161 ± 83 161 ± 83 165 ± 82 165 ± 82

-0.10 -0.05 -0.01 -0.01 -0.01 0.13 0.01 0.03 0.03 0.03 0.01 -0.17 -0.14 -0.14 -0.14 -0.05 -0.22 -0.14 -0.05 0.01 0.09 0.14 0.01 0.18 0.18 0.10 0.09 0.09 0.23 0.23

8

Table S6. Summary of the Generalized Least Squares regressions model. We used it to forecast site chronologies of tree-ring width indices (TRWip). Abbreviations: wi, winter; au, autumn; sp, spring; su, summer. Numbers after climate variables indicate months, whereas the subscript “p” indicates the previous year. Tree species

A. alba

A. borisii-regis

A. cephalonica

A. cilicica

A. pinsapo A. maroccana A. tazaotana

Site code

R2 (adj)

Equation

ROSE PESC IANN CUGN PANT GARA GUAR PORH PORL SJPE BOUM MONT ABRH ABRM ABRL KARP1 PERT AINO KARP2 PARN leb2 syr2 tu24 DEAC IBAC APSB APSG APSN ABMC ABTZ

0.74 0.84 0.70 0.73 0.75 0.77 0.88 0.69 0.55 0.56 0.73 0.64 0.48 0.32 0.57 0.33 0.76 0.81 0.52 0.77 0.85 0.57 0.62 0.65 0.60 0.54 0.89 0.77 0.58 0.57

1.29-0.05t1+0.03t3-0.03t6-0.04t10p+0.09pwi+0.04psp+0.15e9+0.24e27+0.05e21p 1.26-0.07t1-0.02t10p-0.15pwi-0.08psp+0.03psup+0.04paup-0.14e1-0.09e5-0.13e21+0.19e27+0.06e16p+0.04e24p 0.99-0.05t10+0.07tsp-0.04t10p+0.07p5+0.05psu-0.04e17+0.06e25+0.09e24p 0.95+0.05t4-0.05t8-0.05t9p-0.10p1+0.05p5-0.07e5-0.22e8+0.23e22+0.06e24+0.36e9p+0.21e24p 0.97-0.04t10-0.06t10p-0.02p2+0.05psup-0.08e3p+0.07e5p-0.05e15p+0.10e24p 1.00+0.03t3-0.03p1-0.04p10-0.03e3-0.02e16-0.04e17+0.10e26-0.03e18p+0.06e20p 0.94-0.04t1+0.07t2+0.03t5+0.03t11p-0.03p9+0.06pau+0.10e8-0.12e16+0.04e18p+0.05e25p 0.92-0.07t1+0.06t2+0.02t3-0.08t10+0.08t11p-0.05taup+0.06p4+0.03psu-0.11e8p+0.01e18p 0.95 -0.04t1+0.06t2+0.03t7-0.02t8-0.02t9p-0.06taup+0.04p6-0.09e5+0.04e18+0.01e21p 0.66+0.05t2+0.03t5+0.04t11p-0.06taup+0.05psu-0.30e2+0.22e6-0.85e14-0.13e23-0.31e3p-0.09e8p 1.09+0.03t6-0.03t10+0.08tsp-0.04t9p+0.02p8+0.01paup-0.92e2+1.16e13-0.11e23-0.29e8p+0.06e17p 1.36+0.04t4-0.02tsp-0.01t9p+0.03t12p+0.06p2-0.03p10+0.03psup-0.07e2+0.97e13+0.14e15-0.34e6p 1.09+0.03t4-0.02t9-0.04tsu-0.03t9p+0.02p7+0.17e7p-0.05e18p+0.03e21p 0.91+0.09psu+0.02p11p+0.03psup-0.13e7 0.97-0.02t10+0.06t12p-0.07taup+0.04p5+0.04p7-0.02e16-0.05e23+0.17e8p-0.18e15p-0.05e18p+0.07e24p 0.98-0.03p4+0.05p6-0.03e19+0.01e25p 0.89+0.03t3-0.04t10p+0.07p5+0.05psu+0.03p9p-0.43e9-0.04e16-0.52e2p-0.08e5p-0.26e10p-0.08e23p 0.98-0.04twi-0.06p3+0.04p6-0.02psup-0.05e3+0.03e5-0.14e17+0.18e26-0.07e16p-0.04e17p+0.09e23p+0.11e27p 0.92 +0.02t9-0.02p4+0.06p6-0.16e9+0.04e18-0.08e19-0.08e11p-0.02e15p+0.04e24p 1.06-0.05t6-0.02t11p+0.02p2+0.02p6-0.02p12p-0.47e13-0.09e20+0.05e24+0.08e26+0.42e2p-0.14e5p+0.08e24p 0.92-0.02t3+0.04t7-0.06twi+0.03psp-0.13p11p+0.10paup-0.12e16+0.04e20+0.07e24-0.05e15p 0.92-0.07p10-0.02p12p+0.18e5-0.55e6+0.01e21+0.04e25+0.10e27-0.27e15p 1.10+0.03t4-0.04taup+0.05p8+0.01pau-0.06p11p-0.18e1-0.07e3-0.07e19-0.09e24+0.13e26+0.17e7p 1.18-0.04t9p+0.08t12p+0.08p1+0.06p5-0.08e16-0.07e18+0.05e19-0.07e5p-0.47e14p-0.11e23p 1.03+0.05t4-0.05t9p+0.06p5+0.04p7+0.08e26+0.16e27-0.24e22p+0.12e24p 1.23-0.04t5-0.06t6+0.06p7-0.05p8+0.06p10p-0.22e23+0.73e14p+0.19e18p-0.11e19p-0.13e23p 1.02-0.03t5-0.03t8-0.02t11p+0.03p4+0.04paup-0.43e2-0.04e3+0.06e5-0.18e10+0.05e17-0.06e19+0.11e24-0.11e15p-0.06e16p 1.06-0.04t3-0.03t10p-0.03p2+0.09psp+0.07p9p-0.61e2+0.17e6-0.23e7-0.10e20-0.08e23+0.13e24+0.48e9p-0.13e16p 0.97-0.05taup+0.01e24+0.03e25-0.05e1p+0.05e19p-0.08e21p-0.09e23p 0.93+0.04t2+0.06t4+0.04t7+0.03p4+0.04psu+0.05p10p+0.03psup-0.06e17+0.10e18-0.08e23p

9

Table S7. Statistics of the Bayesian estimation of the site-by-site tuned VS-Lite growth response parameters (T1, T2, M1, and M2) considering two 21st century periods; and the two IPCC AR5 emission scenarios (RCP 4.5 and RCP 8.5). Different letters indicate significant (P < 0.05) differences between periods based on Tukey HSD post-hoc tests. IPCC emission scenarios RCP 4.5 Tree speciescountry

Abies alba-Italy

Abies alba-Spain

A. borisii-regis

A. cephalonica

Parameter (unit)

RCP 8.5

2010-2049

2050-2100

2010-2049

2050-2100

Mean

Min

Max

SD

Mean

Min

Max

SD

Mean

Min

Max

SD

Mean

Min

Max

SD

T1 (ºC)

4.501b

3.261

6.892

1.340

5.896a

3.687

8.581

1.657

4.067b

3.072

5.525

0.873

5.358a

4.632

6.584

0.719

T2 (ºC)

10.869b

10.249

12.008

0.625

15.022a

11.294

21.247

4.121

14.957a

10.635

22.833

4.665

13.188b

10.763

21.719

4.374

M1 (v/v)

0.058a

0.011

0.086

0.031

0.037b

0.0.2

0.084

0.027

0.052a

0.022

0.085

0.023

0.028b

0.000

0.061

0.026

M2 (v/v)

0.214

0.124

0.264

0.047

0.244

0.126

0.371

0.098

0.233a

0.148

0.322

0.076

0.149b

0.103

0.173

0.030

T1 (ºC)

3.154b

2.318

4.246

0.176

4.707a

2.881

6.156

1.308

3.407

1.908

6.671

1.940

3.491

2.133

6.868

1.845

T2 (ºC)

12.850b

10.388

23.426

5.192

13.703a

10.538

22.481

4.564

13.905b

10.917

22.995

4.569

14.598a

10.209

23.608

4.996

M1 (v/v)

0.043b

0.002

0.078

0.033

0.056a

0.017

0.078

0.029

0.057b

0.002

0.082

0.030

0.070a

0.007

0.095

0.032

M2 (v/v)

0.279

0.100

0.366

0.126

0.280

0.161

0.357

0.094

0.280

0.107

0.378

0.113

0.230

0.110

0.305

0.091

T1 (ºC)

6.662

3.148

8.581

2.150

6.843

5.258

8.238

1.257

5.359

2.152

8.626

2.989

5.329

1.628

7.277

2.362

T2 (ºC)

14.852a

10.237

23.108

6.067

13.656b

10.236

22.484

4.997

13.704b

10.420

23.419

5.531

12.416a

10.322

18.197

3.336

M1 (v/v)

0.042b

0.013

0.075

0.026

0.055a

0.018

0.097

0.039

0.077

0.015

0.094

0.035

0.065

0.017

0.097

0.033

M2 (v/v)

0.354a

0.313

0.381

0.030

0.277b

0.190

0.355

0.081

0.331a

0.307

0.363

0.027

0.287b

0.155

0.480

0.138

T1 (ºC)

4.623b

3.353

6.738

1.844

6.031a

3.428

7.403

2.255

5.617a

4.306

7.240

1.492

5.310b

3.511

7.541

2.050

T2 (ºC)

11.412

10.965

12.084

0.593

11.261

10.465

12.555

1.131

11.549b

10.440

13.269

1.511

13.177a

10.465

18.022

4.205

M1 (v/v)

0.051a

0.028

0.080

0.027

0.030b

0.007

0.074

0.038

0.077a

0.064

0.091

0.013

0.034b

0.001

0.064

0.032

M2 (v/v)

0.392a

0.310

0.440

0.072

0.265b

0.214

0.298

0.045

0.225a

0.183

0.305

0.069

0.171b

0.103

0.215

0.060

10

Table S7. continued. IPCC emission scenarios RCP 4.5 Tree species

A. cilicica

A. pinsapo

A. maroccana

A. tazaotana

Parameter (unit)

RCP 8.5

2010-2049

2050-2100

2010-2049

2050-2100

Mean

Min

Max

SD

Mean

Min

Max

SD

Mean

Min

Max

SD

Mean

Min

Max

SD

T1 (ºC)

5.737b

4.489

6.088

0.699

6.769a

4.633

8.675

1.841

6.058b

4.437

7.277

1.089

6.806a

3.699

8.163

1.770

T2 (ºC)

16.134a

11.065

22.237

4.361

15.092b

10.993

19.858

3.315

16.198b

13.159

20.680

2.797

19.652a

12.635

23.109

4.219

M1 (v/v)

0.053a

0.005

0.089

0.037

0.030b

0.001

0.077

0.034

0.048a

0.008

0.096

0.038

0.027b

0.001

0.067

0.032

M2 (v/v)

0.106b

0.101

0.111

0.004

0.157a

0.102

0.348

0.107

0.119

0.101

0.164

0.026

0.125

0.100

0.175

0.034

T1 (ºC)

6.690

6.223

6.973

0.407

6.720

5.851

7.172

0.753

6.243b

5.531

7.295

0.930

6.965a

6.183

7.885

0.859

T2 (ºC)

13.312b

13.279

21.221

4.291

15.878a

12.686

21.894

5.213

18.134b

13.000

20.739

4.446

18.510a

14.007

22.639

4.328

M1 (v/v)

0.032

0.003

0.059

0.028

0.038

0.001

0.093

0.049

0.030b

0.004

0.082

0.045

0.058a

0.010

0.093

0.043

M2 (v/v)

0.216a

0.104

0.339

0.118

0.110b

0.100

0.121

0.011

0.145a

0.101

0.168

0.038

0.117b

0.106

0.130

0.012

T1 (ºC)

7.001a

−

−

−

4.135b

−

−

−

5.303b

−

−

−

5.928a

−

−

−

T2 (ºC)

12.865

−

−

−

12.518

−

−

−

11.855b

−

−

−

21.294a

−

−

−

M1 (v/v)

0.069a

−

−

−

0.048b

−

−

−

0.095

−

−

−

0.090

−

−

−

M2 (v/v)

0.199a

−

−

−

0.100b

−

−

−

0.131

−

−

−

0.111

−

−

−

T1 (ºC)

3.626b

−

−

−

7.099a

−

−

−

7.603a

−

−

−

4.648b

−

−

−

T2 (ºC)

11.732b

−

−

−

17.522a

−

−

−

15.472a

−

−

−

12.685b

−

−

−

M1 (v/v)

0.004b

−

−

−

0.035a

−

−

−

0.012b

−

−

−

0.090a

−

−

−

M2 (v/v)

0.170a

−

−

−

0.128b

−

−

−

0.352a

−

−

−

0.101b

−

−

−

11

Table S8. Extreme Climate Indices (for details see http://etccdi.pacificclimate.org/list_27_indices.shtml). Code 1 2 3 4 5

ID FD SU ID TN

Indicator Name

Indicator Definitions

Units

Frost days

Annual count when daily minimum temperature <0 °C

days

Summer days

Annual count when daily max temperature >25 °C

days

Ice days

Annual count when daily maximum temperature <0 °C

days

Tropical nights

Annual count when daily minimum temperature > 20 ºC

days

GSL

Growing season length

Annual (1st Jan to 31st Dec in NH, 1st July to 30th June in SH) count between first span of at least 6 days with TG >5 °C and first span after 1st July (1st January in SH) of 6 days with TG < 5 °C (where TG is daily mean temperature)

days

6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

TXx TNx TXn TNn TN10p TX10p TN90p TX90p WSDI CSDI DTR Rx1day Rx5day SDII R10mm R20mm R1mm CDD CWD R95p R99p PRCPTOT

Hottest day

Monthly maximum value of daily max temperature

°C

Warmest night

Monthly maximum value of daily min temperature

°C

Coldest day

Monthly minimum value of daily max temperature

°C

Coldest night

Monthly minimum value of daily min temperature

°C

Cool nights

Percentage of time when daily min temperature < 10th percentile

%

Cool days

Percentage of time when daily max temperature < 10th percentile

%

th

Warm nights

Percentage of time when daily min temperature > 90 percentile

%

Warm days

Percentage of time when daily max temperature > 90th percentile

%

Warm spell duration index

Annual count when at least six consecutive days of max temperature >90th percentile th

days

Cold spell duration index

Annual count when at least six consecutive days of min temperature < 10 percentile

days

Diurnal temperature range

Monthly mean difference between daily max and min temperature

°C

Max 1  day precipitation amount

Monthly maximum 1 day precipitation

mm

Max 5  day precipitation amount

Monthly maximum consecutive 5 days precipitation

mm

Simple daily intensity index

Ratio of annual total precipitation to the number of wet days (≥ 1 mm)

mm/day

Number of heavy precipitation days

Annual count when precipitation ≥10 mm

days

Number of very heavy precipitation days

Annual count when precipitation ≥ 20 mm

days

Number of precipitation days

Annual count when precipitation ≥ 1 mm

days

Consecutive dry days

Maximum number of consecutive days when precipitation < 1 mm

days

Consecutive wet days

Maximum number of consecutive days when precipitation ≥1mm

days

Very wet days

Annual total precipitation from days > 95th percentile

mm

Extremely wet days

Annual total precipitation from days > 99th percentile

mm

Annual total wet day precipitation

Annual total precipitation from days ≥ 1 mm

mm

12

Table S9. Spearman (rs) correlation coefficients (see color legend) calculated between indexed ring-width chronologies (TRWi) and current (A) and previous (B) annual Extreme Climate Indices (see SI Appendix, Table S8 for indices codes) considering the studied Mediterranean fir species (ABAL, A. alba; ABBR, A. borisii-regis; ABCE, A. cephalonica; ABCI, A. cilicica; ABPI, A. pinsapo; ABMC, A. maroccana; ABTZ, A. tazaotana). Crosses indicate significant correlations (P < 0.05). (A)

ABAL

ABBR

ABCE

ABCI

ABPI ABMC ABTZ

Extreme Climate Indices* Site code ROSE PESC IANN CUGN PANT GARA GUAR PORH PORL SJPE BOUM MONT ABRH ABRM ABRL KARP1 PERT AINO KARP2 PARN leb2 syr2 tu24 DEAC IBAC APSB APSG APSN ABMC ABTZ

e1

e2

e3

e4

X

X

e5

e6

e7

e8

e9

e10

e11

e12

e13

e14

e15

e16

e17

e18

e19

e20

e21

e22

e23

e24

e25

X X

X

X

e26

e27

X

X

X

X

X X

X

X X

X

X

X

X

X

X

X

X

X

X

X

X

X

X

X

X

X

X

X X

X X

X

rs > 0.4 0.3 < rs < 0.4 0.2 < rs < 0.3 0.1 < rs < 0.2 0 < rs < 0.1 0 0 < rs < -0.1 -0.1 < rs < -0.2 -0.2 < rs < -0.3 -0.3 < rs < -0.4 rs < -0.4

X

X

Legend

X X

X

X X X X

X

X

X

X

X

X

X

X X X

X X

X

X

X

X

X

X X

X

X

X

X

X

X

X

X X

X

X

X

X

X X

X X

X

X

X

X

X

X

X

X X

X

X

X

X

X

X

X

X

X

X

X

X

13

(B)

ABAL

ABBR

ABCE

ABCI

ABPI ABMC ABTZ

Extreme Climate Indices* Site code ROSE PESC IANN CUGN PANT GARA GUAR PORH PORL SJPE BOUM MONT ABRH ABRM ABRL KARP1 PERT AINO KARP2 PARN leb2 syr2 tu24 DEAC IBAC APSB APSG APSN ABMC ABTZ

e1

e2

e3

e4

e5

e6

e7

e8

e9

e10

e11

e12

e13

e14

e15

e16

e17

e18

e19

e20

e21

e22

e23

e24

X

e25

e26

X

X

e27

X X

X

X

X

X

X

X X

X

X

X

X X X X X

X X

X

X X

X

X

X

X X

X

X

X

X

X

X

X

X

X

X

X

X

X

Legend rs > 0.4 0.3 < rs < 0.20.4 < rs < 0.10.3 < rs < 0 <0.2 rs < 0.1 0 0 < rs < -0.10.1 < rs < -0.2-0.2 < rs < -0.3-0.3 < rs < rs -0.4 < -0.4

X X X

X

X

X X

X X X

X

X

X X

X

X

X

X

X

X

X

X

X

X

X

X

X

X

X X

X X

X

X

X X

X

X

X

X

X

X

X

X X

X

X

-X

X

X

X

X X

X

X

X

X

X

X

X

X

X

X

X

X

-X

X

X

*e1: FD; e2: SU; e3:ID; e4: TR; e5: GLS;e6: TXx; e7: TNx, e8: TXn; e9: TNn; e10: TN10p; e11: TX10p; e12: TN90p; e13: TX90p; e14: WSDI; e15: CSDI; e16: DTR; e17: Rx1day; e18: Rx5day;e19: SDII; e20: R10mm; e21: R20mm; e22: R1mm; e23: CDD; e24: CWD; e25: R95pTOT; e26: R99pTOT; e27: PRCPTOT

14

Supplementary Figures

Fig. S1. Distribution of the sampled Mediterranean fir forests and climatic diagram for each study site. Total monthly precipitation (P) and mean monthly temperature (T) corresponding to the 1950-2010 period (see site codes in SI Appendix, Table S1).

15

Fig. S2. Pearson correlation coefficients calculated between TRWi and climate variables. Pearson correlation coefficients calculated between site chronologies of tree-ring width indices (TRWi) and mean temperature (cross symbols) and total precipitation (square symbols) for the studied Abies species in the common period 1950-2010. The analyzed temporal window spans from previous (capital letter) July to current October and seasonal periods. The values located outside the dashed lines indicate significant (P < 0.05) correlation coefficients. 16

Fig. S3. Tree-ring width indices (TRWi) calculated for each fir species and extreme growth years (Cropper index, right y axis). Means are shown as dark green lines and grey areas represent the maximum and minimum TRWi values considering all studied sites per species. The lower plots show negative pointer years (bars) expressed as Cropper indices. The red bars and numbers indicate selected climatically extreme years causing a sharp growth reduction. Horizontal lines frame extreme Cropper values. 17

Fig. S4. Percentage of Growth Change (%) to the observed extreme negative years. Percentage of Growth Change (%) to the extreme negative years and percentages of trees showing drought- or heat-induced pointer years during the period 1950-2010 for each of the studied Abies species. 18

Fig. S5. Growth projections according to RCP 4.5 emission scenarios (2). Projected tree-ring width indices (TRWip) for each Abies species site

considering the emission scenario RCP 4.5. Species means are shown as black lines and grey polygons represent the TRWip maximum and minimum values considering all sites per species. The graph below shows the projected change in TRWi reduction (GR) computed considering the 5-year moving average values from 2011 to 2100 period. The red bars and numbers indicate selected climatically extreme years causing reduced growth in each species during the 21st century. Bold lines frame extreme values (GR > 25%).

19

Fig. S6. Maps displaying Pearson correlation coefficients calculated at site level between observed (TRWi) and VS-Lite model projected (TRWiVSL) tree-ring indices for Abies species. Plots show the entire period of analyses (1950–2010, green symbols) and two sub-periods (1950–1980, orange symbols; and 1980–2010, red symbols). Correlation values higher than 0.25 are significant at P < 0.05. 20

Fig. S7. Simulated monthly growth response curves (gT, gM) for the period 1950-2010. Simulated (VS-Lite, period 1950-2010) monthly growth response curves considering temperature (gT, red lines-mean and areas-SD) and soil moisture limitations (gM, blue lines-mean and areas-SD) for each species. Selected extreme events are indicated in different colors (see Figs. 2 and 3). See SI Appendix, Table S1 for site codes. 21

Fig. S8. Projected monthly growth response curves (gT, gM) for the period 2050-2100 for RCP 4.5 and RCP 8.5 emission scenarios. Projected monthly growth response curves corresponding to the RCP4.5 (upper map) and RCP 8.5 emission scenarios (below map) with VS-Lite for the period 2050-2100 considering temperature (gT, red lines-mean and areas-SD) and soil moisture limitations (gM, blue lines-mean and areas-SD) for each selected site. Selected extreme events are indicated in different colors (see Figs. 4 and 5). Map displaying projected growth trends (Pearson correlation coefficients of site mean tree-ring width indices series) considering the 2011–2100 period following growth projections based on respective emission scenarios. Values higher than 0.25 and lower than 0.25 are significant at P < 0.05. See SI Appendix Table S1 for site codes. 22

Fig. S8. Continued.

23

Fig. S9. Simulated winter (gT) and summer (gM) growth responses for the observed and forecasted periods. Simulated (VS-Lite) mean growth response representing winter (DJF) temperature (gT, red circles-mean and SD) and late summer (JAS) soil moisture limitations (gM, blue triangles-mean and SD) for each Abies species. Values were fitted for one observed 1950-2010 period and two 21st century (2011-2049 and 20502100) periods based on RCP 8.5 emission scenario. Selected observed growth responses during extreme events (1950-2010 period, see Figs. 2 and 5) with extreme low growth rates are indicated in different colors on the left. These values defined the thresholds of vulnerability for each climate refugia under forecasted growth response for RCP 8.5 climate scenario. Different letters indicate significant (P < 0.05) differences between observed and projected periods based on Tukey’s HSD post-hoc tests (lowercase and uppercase letters correspond to the 2011–2049 and 2050–2100 periods, respectively). 24