Climate change-induced background tree mortality is exacerbated towards the warm limits of the species ranges

Abstract

Key message

An influence of the recent changes in temperature or rainfall was demonstrated, increasing background tree mortality rates for 2/3 of the 12 studied tree species. Climate change-induced tree mortality was exacerbated towards the warm or dry limits of the species ranges, suggesting in these areas a progressive replacement by more xeric species.

Context

Despite the identification of climate change effects on tree mortality in various biomes, the characterization of species-specific areas of vulnerability remains poorly understood.

Aims

We sought to assess if the effects of temperature and rainfall changes on background tree mortality rates, which did not result from abrupt disturbances, were linked to climate change intensity only, or if they also depended on the tree’s location along climatic gradients.

Methods

We modelled background mortality for 12 of the most common European tree species using 265,056 trees including 4384 dead trees from the French national forest inventory. To explain mortality, we considered variables linked to tree characteristics, stand attributes, logging intensity and site environmental characteristics, and climate change effects.

Results

We found an influence of temperature and rainfall changes on 9 species out of 12. For 8 of them, climate change-induced tree mortality was exacerbated towards the warm or dry limits of the species ranges.

Conclusion

These results highlight that tree mortality varies according to the climate change intensity and the tree location along temperature and rainfall gradients. They strengthen the poleward and upward shifts of trees forecasted from climate envelope models for a large number of European tree species.

1 Introduction

The extent to which climate change has affected the geographical distribution of species has been widely studied, with reports of upward and poleward shifts of many plant species worldwide (Parmesan and Yohe 2003). Climate change-induced range shifts in tree species were detected in different areas. For example, in Western Europe; tree seedlings have shifted in elevation during the last half of the twentieth century (Lenoir, Gégout et al. 2009). In North-East Spain, Fagus sylvatica stands have moved upward in elevation and are being progressively replaced at lower elevations by Quercus ilex, a Mediterranean xeric species (Peñuelas, Ogaya et al. 2007) while in the Swiss Alps, Quercus pubescens stands are progressively replacing Pinus sylvestris stands on dry sites (Rigling, Bigler et al. 2013).

Forecasts from climate envelope models with climate change scenarios at the European scale indicate important changes in species composition and habitat shifts by the end of the twenty-first century, with either winning (with a net gain in potential range) or losing (with a net loss in potential range) species (Dyderski, Paź et al. 2018). These potential changes of European forest species composition could result in important ecological and economic loss as compared to present, due to widespread shifts from high value productive timber species to low productivity Mediterranean species (Hanewinkel, Cullmann et al. 2012).

Tree death is a difficult phenomenon to predict because it is rarely the consequence of a unique factor but rather of a variety of linked factors influencing death at different time scales and intensities (Waring 1987). Catastrophic mortality, which is observed in a stand after severe disturbances like wildfires, extreme climatic events, storms, or pest outbreaks, can be distinguished from background mortality, which is the mortality observed in a stand in the absence of abrupt disturbances (Das, Stephenson et al. 2016). Drought stress can lead to tree mortality by hydraulic failure, occurring at different levels of loss of xylem hydraulic conductance depending on species that have contrasting tolerance levels, or carbon starvation (Choat, Brodribb et al. 2018). Heat stress can also affect the ecophysiological functioning of trees by reducing photosynthetic activity and damaging tissues, with different resistance levels depending on species and tree ontogeny (Niinemets 2010; Teskey, Wertin et al. 2015).

Several studies reported an effect of climate change on mortality associated with changes in temperature and rainfall over the past decades (Carnicer, Coll et al. 2011; Luo and Chen 2015; Senf, Pflugmacher et al. 2018). They highlighted climate change-related excess mortality rates on a large number of species based on correlations between climate change patterns and spatial variability of background mortality (Taccoen, Piedallu et al. 2019). While these studies provided insights into the link between long-term temperature or rainfall changes and tree mortality, they did not analyze if the impact of warming and drying on background tree mortality varied according to the tree location along its climatic gradient.

What happens to trees toward the warm edges of the species ranges is not widely understood. Bussotti, Pollastrini et al. (2015) listed 4 possible outcomes for trees located in range contractions areas in an extensive review: (1) persistence and acclimatization thanks to phenotypic plasticity, (2) evolution through genotype changes, (3) migration, and (4) extinction processes. The centre–periphery hypothesis assumes that populations are at higher risk of extinction near the edge of their geographic range (Brown 1984). Several studies reported decreased growth of European tree species located at the warm edges of their distribution range (Herrero, Rigling et al. 2013; Charru, Seynave et al. 2017). However, observations of increased tree mortality in these areas are rare and only concern a limited number of species and trees located on the driest sites. In the inner Alps, excess mortality rates of Pinus sylvestris were evidenced on the driest sites of a valley in Switzerland, while Quercus pubescens abundance increased on the same sites between 1983 and 2003 (Rigling, Bigler et al. 2013). However, species distribution limit is not always equivalent to species vulnerability to climate change. Genetic, demographical and trait features interact with climatic effects and may allow edge populations to persist, while tree populations living at the distribution core may experience mortality under severe climatic extremes (Vilà-Cabrera, Martínez-Vilalta et al. 2011).

This article complements a previous study that evaluated the impact of climate change on a large set of European tree species, and demonstrated increased tree mortality with temperature increase or rainfall decrease for 18 of the 43 studied species (Taccoen, Piedallu et al. 2019). Here, we aimed to evaluate if mortality of 12 important European angiosperm and gymnosperm tree species was driven by temperature and rainfall changes only or if it varied according to tree location along climatic gradients. We modelled tree mortality using 24 candidate predictors characterizing tree characteristics, stand attributes, logging intensity, and local environmental conditions. Using 12 additional variables describing changes in temperatures and rainfall over the past decades, we evaluated for each species if a climate change effect remained significant on tree mortality when all these other potential factors of mortality were taken into account. When an impact of climate change was detected, we assessed if it was exclusively linked to warming or drying intensity, (hypothesis H1, Fig. 1); or if the climate change effect was aggravated towards the warm or dry limits of a species’ distribution range (hypothesis H2, Fig. 1). Finally, we analyzed the spatial variability of climate change effects on background tree mortality and provided maps of climate-change-related excess probabilities of mortality.

Fig. 1

Schematic representation of potential climate change effects on tree mortality for spatially homogeneous or heterogeneous warming or drying for a given tree species

2 Materials and methods

2.1 Species selection

We focused on the major European forest tree species in terms of forest cover (Brus, Hengeveld et al. 2012). We discarded among these species for which massive and unpredictable health issues are documented in France: Castanea sativa, Fraxinus excelsior and Fraxinus angustifolia (Bakys, Vasaitis et al. 2009). We obtained a list of 12 species including 6 gymnosperms: Abies alba Mill. (Abal), Larix decidua Mill. (Lade), Picea abies (L.) Karst (Piab), Pinus pinaster Ait. (Pipi), Pinus sylvestris L. (Pisy), and Pseudotsuga menziesii Mirb. (Psme), and 6 angiosperms species: Betula pendula Roth (Bepe), Carpinus betulus L. (Cabe), Fagus sylvatica L. (Fasy), Quercus robur L. (Quro), Quercus petraea (Mattus.) (Qupe), and Robinia pseudoacacia L. (Rops).

2.2 Forest area and selection of plots

We used data collected over the afforested areas of France (551 500 km²), a territory that covers 10 out of the 14 climate types found in western Europe according to the Koppen-Geiger classification (Peel, Finlayson et al. 2007). We used non-permanent forest inventory plots inventoried by the National Geographic Agency (IGN) over the period 2009-2015. Approximately 6500 plots were sampled yearly at different locations (Robert, Vidal et al. 2009). They were uniformly distributed across French forests in a 1-km²-based grid (Hervé, Wurpillot et al. 2014). Trees with a diameter at breast height (dbh) > 7.5 cm were considered countable. Trees were surveyed in concentric subplots of 6, 9, and 15 m radius, for trees of 7.5 < dbh < 22.5 cm; 22.5 ≤ dbh < 37.5 cm, and 37.5 cm ≤ dbh, respectively. A dead tree was recorded when it had a dbh > 7.5 cm, had no living parts over 1.3 m height, and was presumed to have died within the last 5 years prior to inventory (the death date was visually assessed). In order to avoid confusion between mortality factors related to abrupt disturbances not necessarily related with climate change and those related to gradual changes of temperature and rainfall regimes, all plots with evidence of perturbations were removed from the analysis. Because we focused on the effects of temperature and rainfall, most of the dead trees uprooted by wind, felled, or broken, representing 1.8% of the trees (versus 3.2% for standing dead trees), were also discarded. Additionally, as the proximity of forest edges induces different growth conditions from closed forest (Harper, Macdonald et al. 2005), all plots with edges and small groves were removed from the dataset (corresponding to approximatively 10% of the plots). Plots composed of less than 10% of living trees with a diameter > 7.5 cm were removed due to the impossibility to calculate reliable synthetic indices such as basal area (BA). The final number of trees was 265,056 trees, including 4384 dead trees (Table 1, Appendix Text 1 and Fig. 7).

Table 1 List of the species. Code corresponds to the abbreviations of the genus and species names. Plots with dead trees contained at least one dead tree with a dbh > 7.5 cm

2.3 Candidate variables tested in the models

We used 24 variables in order to take well-known effects on tree mortality into account (Taccoen, Piedallu et al. 2019). They concerned tree characteristics, stand attributes, logging intensity, soil properties, and reference period climate (Table 2). In addition, we compared the performance of 6 variables describing climate change intensity only (Hypothesis H1) and 6 variables combining historical climate and climate change intensity (Hypothesis H2, Table 2).

Table 2 Description of the 36 explanatory variables used in the models Abb = abbreviation. The Source column indicates the origin of the data: collected by the NFI (Field), calculated using the NFI information (Calc) or extracted from models available with geographic Information System (Mod.). See Appendix Table 5 for the values.

2.4 Tree characteristics and stand attributes

Tree characteristics were described by tree circumference at 1.3 m height (Circ), measured on the field, and tree relative circumference (RelCirc). RelCirc is the ratio of the circumference of each individual tree and of the mean of the circumference of all of the trees within the plot. To describe stand density and stage development, we calculated 4 variables using field measured tree attributes. BA is the weighted sum of all trees basal area within the plot (m²/ha), calculated from the measured Circ. The number of trees per hectare (NB) is the sum of the number of trees inventoried on the plot (nb/ha). Canopy cover (CC) is the proportion of the floor covered by the vertical projection of all measured tree crowns on the plot (%). The plot quadratic mean diameter (QMD), representing the weighted mean quadratic diameter of the trees present in the plot, was also computed. To evaluate forest composition, we calculated the total number of tree species inventoried on the plot (Nb_sp) and the proportion of BA occupied by the species within each plot (PropBA, %). To evaluate stand structure heterogeneity, we calculated the Gini index of inequality of tree diameters (Gini, unitless) on the plot (Latham, Zuuring et al. 1998) that varies between 0 and 1 with increasing diameter unevenness. Because trees were not inventoried on the same surfaces, weights were attributed to each tree according to the size of the inventory subplot and were used to calculate all these indices except Circ and Nb_sp.

2.5 Logging intensity

Logging intensity can considerably influence observed mortality (Kahl and Bauhus 2014). Fewer dead trees can be observed in the most accessible and managed forests as a consequence of silvicultural practices leading to the removal of dead and declining trees, and it can reduce competition and also sometimes increase wind throw risk. Skid trails (trails) is an indicator of the presence of already existing skid trails. The skidding distance (skid) is an indicator of the distance from the plot center to the nearest existing skid trail. At last, recent cut (cut) is the intensity, if existing, of recent cuts on the plot (less than 5 years) evaluated with the number of stumps on the plot.

2.6 Soil properties

Unfavorable soil conditions are predisposing factors to tree decline and mortality (Manion 1981). Six variables were used, describing water availability, waterlogging, and soil nutrition. Available Water Content (AWC) represents the maximum amount of water a soil can store and was calculated from the soil surveys made by the NFI on each plot. The textural description, the depth and the proportion of stone content of the soil were used with pedotransfer functions to estimate AWC for each plot (Piedallu, Gégout et al. 2011). pH, carbon to nitrogen ratio (CN), temporary and permanent waterlogging indices (TW and PW) were calculated using indicator values (IV) of understory plant species present in the plot (ter Braak et Looman, 1986) defined using the EcoPlant database (Gégout, Hervé et al. 2003). The floristic surveys made on each plot to estimate pH, CN, TW, and PW by averaging the indicator values of each species recorded on the plot. Finally, Topo aims to estimate the lateral water flow qualitatively using the site topography measured on field in three factors: water departure, water arrival, and equilibrium.

2.7 Climate-related variables

As temperature acts on plant growth and survival during all year and precipitation act as a water-stress driver during the growing season (Niinemets 2010), climate variables were assessed using digital maps at 1km resolution, with variables describing seasonal temperature (Tm) for all seasons of the year and rainfall (RF) for spring and summer seasons. Seasonal values were calculated by averaging monthly values for Tm and by summing average monthly rainfall for RF, with: December to February; March to May; June to August; September to November for winter, spring, summer and autumn, respectively. First, we considered average climate conditions on a historic period. Secondly, we considered the changes of temperature and rainfall between this period and contemporary period (hypothesis H1). Finally, we considered the interactions between average climate conditions and climate evolution (hypothesis H2).

2.8 Reference period climate

Climate in France is characterized by several decades of stable climate conditions before a strong increase of temperature since 1988 (Appendix Fig. 8). We averaged monthly temperature and rainfall values on a reference period spanning the period 1961–1987. Using the methods developed by Ninyerola, Pons et al. (2000) based on regression analysis between the climatic variables and different environmental predictors available with geographical information systems, we modelled and mapped the climatic conditions using 237 homogenized meteorological stations (Gibelin, Dubuisson et al. 2014) for temperature (calculated as the mean between maximum and minimum temperature), and 435 stations for precipitation (Piedallu, Cheret et al. 2019). The monthly maps were calculated using a bootstrap linear modelling approach (Bertrand, Lenoir et al. 2011). The R² values for the cross-validation of mean annual temperature and cumulated annual rainfall observations were found to be 0.93 and 0.77 respectively. From this dataset, we calculated 6 variables for different seasons: Tm_winRef, Tm_sprRef, Tm_sumRef, Tm_autRef, RF_sprRef, and RF_sumRef.

2.9 Climate change intensity

We used variables characterizing climate evolution as compared to the reference period to assess the effects of climate change on mortality independently of the tree location inside the species range (Hypothesis 1). The average temperature in France over the 1988–2015 period was 1.1 °C higher than the average temperature over the 1961–1987 period (Appendix Fig. 8). To avoid combining model uncertainties calculated for different periods, we did not calculate climate evolution based on the models calculated for two different time slices, but by interpolating the changes calculated for all the weather stations (200 stations for temperature and 1119 stations for rainfall), using MétéoFrance data (Gibelin, Dubuisson et al. 2014) (Appendix Figure 9). Only homogenized meteorological stations with complete records over the 1961-2015 period were used. For each meteorological station, we calculated the difference between the 1961 and 1987 reference period and periods contemporary to the NFI sample plots recorded each year between 2009 and 2015. Because a tree can die a long time after climatic disturbances (Das, Stephenson et al. 2013), we considered for each plot 15 years shifting periods of time before the date of survey to fully capture the effects of delayed mortality. Changes were calculated monthly for each meteorological station and interpolated with ordinary kriging using Arcgis 10.5., showing changes according to the seasons (Fig. 2a, b). We obtained 6 variables characterizing climate change evolution (Tm_winEvo, Tm_sprEvo, Tm_sumEvo, Tm_autEvo, RF_sprEvo, RF_sumEvo).

Fig. 2

Spatial representations of climate change anomalies per season between the reference period and the 1988–2015 period for temperature (°C, a) and rainfall (mm, b). c The interactions between reference period temperature and temperature evolution (°C²). d The relative rainfall evolution (%) as compared to the reference period average rainfall

2.10 Interaction between reference period climate and climate change intensity

Finally, we calculated interactions between the reference period climate and climate change intensity to assess if climate change effects were more important toward the warm or dry edges of the species ranges (Hypothesis 2). For temperatures, they were calculated for each season as the product between reference temperatures and their changes (Fig. 2c). For rainfall, they were calculated as the quotient between rainfall changes and reference rainfall for spring and summer (Fig. 2d). We obtained 6 variables characterizing these interactions (Tm_winRef*Evo, Tm_sprRef*Evo, Tm_sumRef*Evo, Tm_autRef*Evo, RF_sprevo/RF_sprRef, RF_sumevo/RF_sumRef).

2.11 Statistical model

We used logistic regression to model the mortality status of each tree (0: alive, 1: dead) for the 12 species using the 36 candidate variables. Variable selection was made for each variable one after the other, with a forward stepwise procedure based on the maximum decrease in residual deviance (McCullagh and Nelder 1989). If the relation between a variable and the presence-absence of mortality was significant at the 0.05 level, it was kept in the model. Else, the following variable in terms of residual deviance decrease was tested. We continued the variable selection process until no variable added a significant deviance reduction. To assess for the commonly observed U-shaped and bell-shaped response curves to tree size and competition intensity (Holzwarth, Kahl et al. 2013), we tested tree stand attributes variables with their quadratic forms (Ter Braak and Looman 1986). Model quality was evaluated with the area (AUC) using the receiver operating characteristics (ROC) curve (Fielding and Bell 1997) and True Skill Statistic (Allouche, Tsoar et al. 2006). We cross validated the models using 80% of the trees for calibration and 20% of randomly selected trees for validation. To evaluate the relative importance of each predictor in the models, we calculated the drop contribution of each variable i (Marmion, Luoto et al. 2009). It is the difference in explained deviance when the variable is dropped from the multivariate model:

$$ {\mathrm{DropContribution}}_i={\mathrm{ResidualDeviance}}_{\mathrm{wit}\mathrm{hout}\_i}-{\mathrm{ResidualDeviance}}_{\mathrm{wit}{\mathrm{h}}_i} $$

For each variable, the ratio of its drop contribution and the sum of the drop contributions of all the variables in the multivariate model was calculated, representing the relative importance of the variable in the model (RI):

$$ \operatorname{Re} lative\operatorname{Im} por\tan {ce}_i=\frac{DropContribution_i}{\sum_{i=1}^n DropContribution} $$

With

n = number of variables in the model

i = Variable in the model

We also calculated and mapped the excess probability of mortality as compared to a climate change-free context for each tree of each species, as the difference between mortality predictions from the model including all the variables and mortality predictions from the same model where effects related to climate change were set to 0. All non-climate change-related variables were set to their mean value for all the trees of the considered species (Muller and MacLehose 2014). We obtained, for each plot where the species was present, a single value in %, corresponding to the changes of probability to find a tree dead since less than 5 years attributed to climate change. We mapped and interpolated these values with ordinary kriging using ArGIS 10.5 in order to produce maps of excess probability of mortality per species.

3 Results

3.1 Model quality and predictability of background mortality

The AUC for the 12 mortality models varied from 0.75 to 0.91 (mean ± s.d of 0.82 ± 0.05) and TSS from 0.42 to 0.69 (0.53 ± 0.10) in the validation dataset (Table 3). Mean AUC of gymnosperms models (0.85) was slightly but not significantly higher than that of angiosperms models (0.83, T test, P = 0.46).

Table 3 Performance of the mortality models for the 12 species in terms of AUC and TSS. Calibration was performed on 80% of the trees and validation on 20%

3.2 Variables selected in mortality models

Mortality models for the 12 species had on average 7.3 variables selected (Table 4). The average relative importance (RI) of the variables highlight that mortality is primarily driven by tree characteristics and stand attributes variables. For all of the 12 species, at least one tree characteristics variable and one stand attribute variable were selected during model building (Fig. 3, Table 4). Soil properties variables and reference period climate variables were each selected for 4 species. Finally, climate change variables were selected for 9 species among 12 (Fig. 3, Table 4).

Table 4 Responses of the 36 variables evaluated from models fitted using logistic regression, for the 12 studied species. The ∪ symbol indicates a significant U-shaped response curve, the ∩ symbol a bell-shaped response curve, the ↘ and ↗ symbols indicates respectively a decreasing and an increasing response curve. For qualitative variables (Dist, Trails, and Cut), a symbol is for coefficient values indicating in average a significant lower mortality with increasing logging intensity. A + symbol is for higher mortality with increasing logging intensity. When there is no symbol the variable was not significant (P > 0.05). RI = Relative importance. Nb= total number of species for which the variable was significant. See Table 1 for the corresponding species name and Table 2 for the corresponding variables name. See Taccoen et al. 2019 for the coefficients of the model (Taccoen, Piedallu et al. 2019).

Fig. 3

Relative importance (RI) of the variables per category for the 12 species mortality models. Species are sorted by decreasing RI of climate change effects inside each group. Here, climate change variables do not distinguish between climate change intensity variables (H1) and interactions between climate and climate change (H2)

3.3 Effects of variables that are not related to climate change

Tree characteristics variables were significant for all species (Table 4). RelCirc (mean RI = 40%) response curves indicated a pattern of greater mortality for the most dominated trees within a stand and slight resumption of mortality for the most dominant trees (Appendix Table 5). Observed trends for Circ (RI = 9%) indicated that the typical dead tree, regardless of the species, was a small and dominated individual. PropBA (RI = 16%) showed different responses while BA (RI = 6%, Table 4) indicated globally higher mortality rates in stands with higher biomass and likely higher competition intensity. The Gini index (RI = 5%) indicated a lower mortality for stands with very irregular and very regular stand structures. At least one of the three Logging intensity variables was significant for 8 of the species (Table 4), with a homogeneous response corresponding to decreasing observed mortality with increasing logging intensity. Significant mortality responses to soil properties variables were scarce and of low relative importance. Temperature variables were significant for 4 species, and indicated higher mortality when summer or autumn temperatures were high (Table 4).

3.4 Effects of climate change

Climate change variables were significant for 9 species out of 12 with an average RI of 5.3% (Table 4). Eleven of the 12 climate change variables selected in the models were interactions between the reference period and climate change, while only one of them concerned climate change anomalies (Table 4, Fig. 4a). This indicates that for a given warming or changes in rainfall regimes, for these species, probabilities of mortality are higher in the warmest or driest areas of the species range. With 6 species concerned, the interactions between temperature increase and mean temperature were the preponderant climate-change-related effects in our models, with more mortality in the warmest and warming areas in spring and summer (Quro, Lade, Rops, Piab, and Psme), autumn (Pipi), and winter (Quro). Bepe expressed a different response to temperature than all the other species, with a decreasing mortality with increasing winter temperature. Effects of changes in relative rainfall amount concerned three species (Fig. 4), with increasing mortality with decreasing relative rainfall in summer for Fasy, in spring for Pisy, and a slight decrease in mortality with decreasing relative spring rainfall for Pipi. Effects of changes in temperature or rainfall without interaction with climate variables only concerned Bepe (Fig. 4).

Fig. 4

Number of occurrences (a) and cumulated relative importance (b, %) of climate change-related variables corresponding to the hypothesis H1 and H2 for the 12 species. In (b), values above the 0 value indicate an increase in mortality in warmer or drier ecological conditions Values below 0 indicate a decrease in mortality (see Table 1 for the corresponding species name)

3.5 Comparing the vulnerability of species to climate change

The excess probability of mortality as compared to a climate change free context was calculated for each species (Fig. 5, Appendix Table 7). It is the highest for Rops (mean = + 1.80%/5 years), Piab (+ 0.98%/5 years), and Quro (+ 0.69%/5 years). Pipi, Psme, and Lade had comparable excess probability of mortality values (with means of + 0.37%/5 years, + 0.30%/5 years, and + 0.29%/5 years, respectively). Pisy and Fasy both had close to zero mean excess probability of mortality values (+ 0.14%/5 years and − 0.05%/5 years, respectively). For Bepe, the mean value reached − 0.64%/5 years.

Fig. 5

Excess probability of mortality (%) as compared to a climate change free context for the 12 species. The climate change-related excess probability of mortality is calculated as the difference between predicted values for the full model and a partial model where climate change-related effects have been set to 0. In order to assess the importance of the sole climate-change effects on tree mortality, all non-climate change-related variables were set to their mean value for all the trees of the considered species, while qualitative variables were set to their most frequent value. The upper and lower boundaries of the boxplots are the 25th and 75th percentiles and the central line is the median (see Table 1 for the corresponding species name)

Finally, we mapped the change of probability of mortality induced by climate change for the 9 species responding to climate change variables (Fig. 6). For four species, Fasy, Pisy, Quro, and Rops, the highest excess probability of mortality rates was reached in southeastern part of their range where both temperature and temperature changes were high. For the three mountainous species, Piab, Lade, and Psme, it was observed in lowland areas, where average temperature is higher than at high elevations. The other two species have different predicted patterns of mortality, with increasing mortality in the colder part of the range for Pipi, and decreasing mortality rates for Bepe.

Fig. 6

Climate change-related excess probability of mortality for the 9 species sensitive to its effects. When changes in temperature and rainfall lead to more mortality, the excess mortality rate is positive (represented by warm colors) and negative when they lead to less mortality (represented by cold colors). Spatial representations were limited to the current distribution area of the species, obtained with indicator kriging using the presence/absence of the species extracted from the National Geographic Institute database

4 Discussion

4.1 Links between climate change and background tree mortality

Climate variables were selected in all the models, with an influence of recent climate change for 2/3 of the studied tree species. Among the 12 species we studied, 6 were specifically monitored at the European scale between 1992 and 2017 as part of the ICP Forests (level 1) damage survey (Michel and Seidling 2017). Important defoliation rates increases were recorded in this period, up to + 30% for Pipi, Quro and Qupe, + 20% for Fasy, and + 10% for Piab, more pronounced in the Iberian Peninsula (Carnicer, Coll et al. 2011). At the European scale, different studies have characterized potential range shifts of the major tree species in the next decades. Upward and poleward shifts were forecasted for most of European species, with potential range shifts from low to high elevations and latitudes (Dyderski, Paź et al. 2018). However, forecasts from climate envelope models are characterized by important variability and uncertainties between the different models (Cheaib, Badeau et al. 2012). Here, we highlighted a link between tree mortality and climate change, corroborating observations of increasing mortality rates and modeling studies demonstrating species shifts toward higher latitudes or elevations. We failed to distinguish between mortality linked to abiotic and biotic factors. Some of the species were more sensitive to pathogens (e.g., Piab, Abal) than others (e.g., Cabe, Fasy). The comparable model performance for these different species suggests that stands with vulnerable ecological conditions could be preferentially threatened by biotic agents; the effect of pathogens was probably indirectly taken into account in our models.

4.2 An increased vulnerability of species toward the warm edges of their distribution

From the 9 species identified to be sensitive to climate change variables, the effects of changes in temperature and rainfall regimes on tree mortality were exacerbated towards the warm edges of the species distribution for 8 of them, corroborating our hypothesis H2. Since the early twenty-first century, several dieback events and trends of tree growth decline were reported in hot and dry areas as a consequence of important droughts and heat waves (Breda, Huc et al. 2006). In South-West France, where Quro grows close to the warm edge of its distribution, its abundance decreased with decreasing water availability, while for Quercus ilex, a xeric species for which this region represents the cold range of distribution, the abundance remained constant (Urli, Lamy et al. 2015). From a physiological point of view, Quro appeared to be much more vulnerable to xylem cavitation than Quercus ilex, suggesting important demographic changes and a future replacement of Quro by Quercus ilex with increasing warming and drought frequencies. Mortality rates of Pisy were also found to be abnormally high in a Swiss valley, where it evolves close to the dry limits of its distribution (Rigling, Bigler et al. 2013). Identically, in Spain, Fasy stands located in warm and dry areas are being progressively replaced by Quercus ilex, for which recruitment rates up to 3 times higher than that of Fasy were measured (Peñuelas, Ogaya et al. 2007). The increasing mortality in the colder part of the range we observed only for Pipi could be explained by important differences in management practices between the South-West of France, where this species is widely used for timber production in mono-specific even-aged plantations, and other parts of its distribution, where this species can be found in mixed and uneven aged forests.

Along altitudinal gradients, different responses of mountain tree species to climate change were expected with potentially improved ecological suitability at high elevation, where temperature is the growth limiting factor, and increased vulnerability at low elevations where water availability is the limiting factor (Niinemets 2010). Such trends were previously evidenced in terms of growth decline at low elevations for several European coniferous species (Charru, Seynave et al. 2017). However, these observations of increased climate-change-related mortality towards the warm edges of the species ranges either occurred after extreme drought events (Breda, Huc et al. 2006) or concerned specific species and sites located at the very warm end of the species range (Vilà-Cabrera, Martínez-Vilalta et al. 2011; Rigling, Bigler et al. 2013). Here, we found a maximum climate change-related mortality in the warmest and driest part of the species distribution ranges for most of the studied species. However, the interaction between the reference period and climate change variables used in our models does not allow us to distinguish between the effects of historical climate and the intensity of changes. For example, significant warming when temperatures were low during the reference period or little warming when temperatures were high during the reference period led to the same values.

4.3 A possible decreasing mortality in response to climate change toward the cold edge of the species distribution

For several species, we highlighted mortality decreases in response to climate change on the totality or their distribution ranges (Bepe) or in specific areas (Fasy, Pisy). For Bepe, we detected a decreasing mortality with increasing winter temperatures. This trend is consistent with growth increases observed in Western Siberia between the late nineteenth century and the early twenty-first century in response to longer vegetation periods and mean temperature increases (Kharuk, Kuzmichev et al. 2014). Identically, forecasts for 2100 in Finland indicate growth increases for Bepe in parallel with decreases for Piab and Pisy (Kellomäki, Strandman et al. 2018). The distribution of mortality decreases for these species is limited to cold and moist sites where available water content is high (Rameau, Mansion et al. 1989). In this context, increasing temperatures in the absence of water deficit could result in lower mortality rates.

4.4 Species apparently insensitive to changes of temperature and rainfall

We did not highlight climate change-related effects on the mortality of Cabe, Abal, or Qupe. Cabe is considered a drought-sensitive species from an ecological point of view (Zapater, Breda et al. 2013). However, the current distribution of this species in France is corresponding to the core of its climate envelope, while the warm edges of its distribution are located at lower latitudes (Sikkema, Caudullo et al. 2016). For Abal, important growth decreases were highlighted in France at low elevations and increases at high elevations, suggesting an increased vulnerability at the species warm range limit (Lebourgeois, Rathgeber et al. 2010). Massive mortality events due to bark beetle outbreaks were also recorded in recent years (Marini, Økland et al. 2017), but they occurred after 2016 and the data used for this study was collected between the years 2009 and 2015, which corresponds to a period with fewer bark beetle attacks in France. Although we did not highlight climate change effects on these three species, they responded significantly to average temperatures in models, indicating that probabilities of mortality were higher in the warmest areas of their distribution ranges. For some species, climate variables describing the reference period could be selected instead of climate change variables, due to their greater spatial variability compared to climate change variables, mainly in mountainous areas. Therefore, future temperature increases might also trigger a climate change response of these species.

4.5 Insights on the link between tree and stand characteristics and background tree mortality

We found results in line with previous background mortality studies, ensured that the main drivers related to stand characteristics were taken into account, and limited the risk of inadequate selection of the climate change variables. The effects of competition and the tree development stage, approximated by the circumference, were the predominant drivers of mortality (Ruiz-Benito, Lines et al. 2013; Taccoen, Piedallu et al. 2019). These effects reflect the self-thinning-related mortality observed in all stands that can induce up to 95% decrease of the total number of trees on a species diameter range, as evidenced in pure even-aged stands in France (Charru, Seynave et al. 2012). Additionally, we highlighted that tree likelihood of mortality was influenced by tree sizes distribution within the stand and species composition (Fridman and Stahl 2001; Eid and Oyen 2003). The use of these indicators should be strongly considered in future mortality studies.

Large-scale tree mortality studies based on forest inventories also face the problem of forest harvesting practices. In western and central Europe forests have been intensively managed for centuries (Hermy and Verheyen 2007). The removal of dead and declining trees leads to a lower number of observed dead trees in situations where logging intensity is high (Nyland 2016). We found that logging intensity was a critical parameter to consider in mortality studies to remove potential bias factors in models. However, due to the difficulty to accurately account for logging intensity in inventory data, the proportion of dead trees is probably underestimated both in inventory databases and in our calculations.

4.6 Mapping tree species vulnerability to climate change

The identification of the main drivers of tree mortality and knowledge about the spatial distribution of recent climate change allowed us to map tree species vulnerability at the national scale. Forest adaptation to climate change needs to address the different risk factors according to the region and to prioritize actions according to local conditions (Lee, Maggini et al. 2015). Vulnerability is most of the time assessed using expert knowledge, species distribution models calculated for different periods of time (Gavrutenko, Gerstner et al. 2021), or climate change vulnerability indices (Williams, Shoo et al. 2008). The present study demonstrates that mortality models can be used for spatial risk assessment. Such an approach has the advantage of reflecting observed tree health in forests over a specific period of time. Species vulnerability maps can be useful tools for policy-makers and forest managers in a context of severe observed diebacks (Allen, Breshears et al. 2015). We only mapped the mortality risk linked to climate change variables, but we could also include spatial data describing reference period climate and soil properties, or even stand characteristics when available. In that case, modifying the stand characteristics in the models can make it possible to simulate the effects of management on the species vulnerability level.

5 Conclusion

We provided an overview of how recent climate change has differentially impacted the mortality rates of the main European forest tree species and how this effect varied along the heat and drought gradients. We also demonstrated the high spatial variability of the responses, with mostly negative effects of climate warming and drying, but also sometimes positive effects linked to winter warming or rainfall increase. Forecasting from climate envelope models showed potential range shifts from low to high elevations and latitudes for European species (Cheaib, Badeau et al. 2012; Dyderski, Paź et al. 2018). Here, through a large-scale statistical study, we showed that a large number of European tree species is concerned by increased climate change-related mortality toward the warm edges of the species distribution.

With model projections by 2100 of global average temperature increases ranging between + 2.0 °C and + 4.9 °C as compared to 1960-2010 period (Raftery, Zimmer et al. 2017), the impact of climate-change on tree mortality is expected to rise, suggesting important changes in future forest species composition toward the warm limits of the species ranges, with a progressive replacement by more xeric Mediterranean species. The identification of areas where climate change-related excess mortality has already occurred opens promising prospects for the monitoring of climate change effects and the adaptation of silvicultural practices in the near future.

Appendices

Appendices

Text 1 Details about the survey periods used for modelling. The NFI survey with the current protocol has started in 2005, but different changes in the protocol occurred. In 2005 and 2006, trees were considered dead if they were unlikely to survive more than 1 year. From 2007 to 2015, each tree without visible living parts above 1.30 m height was considered dead. Moreover, the sampling protocol with concentric circles was applied to all living trees from 2005 to 2015 and to dead trees from 2009 to 2015: from 2005 to 2008, dead trees were measured across the full 15-m circle regardless of their diameter, resulting in a significant drop in the proportion of dead trees in the 2009–2015 period as compared to the 2005–2008 period. Due to these significant changes in the sampling protocol, we excluded the years 2005–2008 and only kept data from the 2009–2015 period for our study.

Fig. 7

Location of the forest plots with living (grey) and dead trees (black) for the 12 studied species. N = total number of trees of the considered species. n = number of dead trees

Fig. 8

Yearly temperature deviation between 1961 and 2017 compared to the 1961–1987 average (°C) in France. Temperature values were extracted on each NFI plot and yearly averaged. We calculated the difference between yearly temperature values and the average 1961–1987 temperature value. The average value for the 1988–2015 period is of + 1.1 °C as compared to the average 1961–1987 average temperature value

Fig. 9

Location of the homogenized climate series used to model climate evolution with respectively 200 and 1119 series for mean temperature (TM) and rainfall (RF). Data used is provided by Météo-France, the national reference institute for climatic data

Table 5 Median, 1st, and 99th centile value of the potential explanatory variables for the 12 species. See (in the main text) Table 1 for the corresponding species name and Table 2 for the corresponding variables name

Table 6 Critical values (i.e., values at which the tangent to the response curve is horizontal) for the variables with U-shaped and bell-shaped response curves. When the response is strictly decreasing or increasing, no values are provided because no inflexion point exists. Critical values have the same units as the variables (see Table 2 in the main text)

Table 7 Importance of climate change effects on the predicted probabilities of mortality of the 12 species. The base mortality rate corresponds to an average predicted value where all non-climate change-related variables were set to their average value on all the trees of the considered species. The 6 columns located in the right part of the table indicate how predicted probabilities of tree mortality in %.5 years⁻¹ vary depending on the intensity of climate change. One value per plot where the species is present is calculated.