Mutually inclusive mechanisms of drought-induced tree mortality

An extreme summer drought caused unprecedented tree dieback across Central Europe in 2018, highlighting the need for a better mechanistic understanding of drought-induced tree mortality. While numerous physiological risk factors have been identified, the principal mechanisms, hydraulic failure and carbon starvation, are still debated. We studied 9,435 trees from 12 temperate species planted in a diversity experiment in 2013 to assess how hydraulic traits, carbon dynamics, pest infestation, tree height and neighbourhood competition influence individual mortality risk. We observed a reduced mortality risk for trees with wider hydraulic safety margins, while a rising sugar fraction of the non-structural carbohydrate pool and bark beetle infestation were associated with higher risk. Taller trees had a lower mortality risk. The sign and magnitude of neighbourhood effects on mortality risk depended on the species-identity of the involved trees, with most species having beneficial and some having deleterious effects on their neighbours. While severe tissue dehydration causing hydraulic failure precedes drought-induced tree mortality, we show that the probability of this event depends on a series of mutually inclusive processes including pest infestation and starch depletion for osmotic adjustment, and is modulated by the size and species identity of a tree and its neighbours.


Introduction
Worldwide, forests are exposed to rises in temperature, atmospheric vapour pressure deficit and intensity and frequency of severe drought events (Dai 2013;Yuan et al. 2019;Zhou et al. 2019).
As a consequence, large-scale tree mortality events are documented for all forest biomes (Allen et al. 2015;Hartmann et al. 2018;Senf et al. 2018;Schuldt et al. 2020). Challenges in forecasting the impacts of drought-related tree mortality based on the structure and composition of forests highlight the need for a robust mechanistic understanding of the processes involved (Choat et al. 2018;Brodribb et al. 2020).
Two interrelated physiological mechanisms are central for drought-induced tree mortality: hydraulic failure, which is the partial or complete loss of xylem functionality due to embolism formation, and carbon starvation, the depletion of non-structural carbohydrates (NSC) due to negative drought impacts on photosynthesis (McDowell et al. 2008). While there is clear evidence that hydraulic failure is a universal component of the processes preceding tree death under drought (Rowland et al. 2015;Anderegg et al. 2016;Adams et al. 2017;Correia et al. 2019;Powers et al. 2020;Li et al. 2020), the role of NSC reserves in these processes is still debated (Adams et al. 2017;Blackman et al. 2019;Kannenberg & Phillips 2020;O'Brien et al. 2020). Both processes interact with antagonistic biotic agent demographics (McDowell et al. 2008(McDowell et al. , 2011, as drought effects on the water and carbon balance may additionally predispose trees to pathogen or insect attacks (Raffa et al. 2008;Sturrock et al. 2011;Hart et al. 2014;Huang et al. 2020), exacerbating climate-driven forest mortality (McDowell et al. 2013(McDowell et al. , 2020Seidl et al. 2017).
Many risk factors have been associated with drought-induced mortality (O'Brien et al. 2017).
Specifically, increasing tree size is often found to increase drought vulnerability (Peng et al. 2011;Bennett et al. 2015;O'Brien et al. 2017;Stovall et al. 2019). Further, water shortage may intensify interspecific tree competition depending on species-specific water use behaviour and above-and below-ground niche partitioning (Ammer 2019; Grossiord 2019), thus affecting individual mortality risk via variation in local neighbourhood composition (Vitali et al. 2018;Fichtner et al. 2017Fichtner et al. , 2020. In the recent past, several distinct mechanisms of drought-induced mortality have been proposed for a variety of tree species in different biomes (Adams et al. 2017;O'Brien et al. 2017;Choat et al. 2018;Brodribb et al. 2020). Mechanistic insight advancing future processbased models is most likely to result from research that accounts for the multitude of simultaneous processes at play during tree death under drought (McDowell et al. 2019).
Consequently, field-based studies that jointly analyse multiple potential causes of droughtinduced tree mortality across a range of functionally different tree species are urgently needed.
In 2018, Central Europe experienced an extreme summer drought, followed by a non-typical dry winter and spring (Hari et al. 2020). This global-change type drought event resulted in the highest average growing season temperature and vapour pressure deficit ever recorded (3.3°C and 3.2 hPa above the long-term average for the reference period from 1961-1990, respectively) and caused unprecedented drought-induced tree mortality in many species (Schuldt et al. 2020).
We used observations from 9,435 trees belonging to six angiosperm and six coniferous species from Europe and North America grown in a young tree diversity and nutrient enrichment experiment established in 2013 in Southern Germany (Wein et al. 2016) to try to untangle the complexity of possible mechanisms explaining drought-induced mortality in trees. Based on a Bayesian hierarchical modelling approach, we predicted the individual mortality risk following the drought using data on hydraulic traits, carbon dynamics, pest infestation, tree height and interactions with neighbouring trees. We ask whether (i) a tree species' susceptibility to drought-induced mortality depends on both hydraulic traits and NSC dynamics, and whether (ii) the probability of dying on the individual level is modulated by tree height and the influence of biotic (pest infestation and interaction with neighbour trees) and abiotic (fertilisation and environmental) factors.

Cross-species patterns of drought-induced tree mortality rates
In response to the anomalous severe drought of the year 2018 and the prolonged dry conditions that followed throughout the winter and spring of 2019, 34% of the 9,435 trees present in our experimental plantation died. Tree mortality increased dramatically compared to a total loss of 9% of trees between 2013 and early 2018 (see Methods). Mortality rates varied strongly across the 12 studied species, ranging from 0.6% for A. saccharum to 79.5% for L. laricina (Table 1, Fig. 1). The average probability of dying predicted by our model closely reflected the specieslevel mortality rates observed in the field (Fig. 1). Our model reached full convergence (Table   S2.1) and explained 51.8% of the variance in the observed individual mortality (Fig. 2).
Mortality rates were credibly elevated over the average mortality rate (i.e. their 95% highest posterior density interval (HDI) excluded the overall average) for Larix spp., Betula spp. and P. abies, while they were much lower specifically for Acer spp. and Q. rubra (credibly below 10% in all three cases).

Hydraulic safety margins and non-structural carbohydrates
Both hydraulic traits and carbohydrate dynamics differed considerably among species, with observed hydraulic safety margins (HSM) ranging from 0.63 to 3.65 MPa, and changes in the fraction of soluble sugars (glucose, fructose and sucrose) to total leaf NSC (Δsugar) over the growing season ranging from -33.2 to +17.0% (Table 1). An increase in relative leaf sugar content is indicative of starch depletion and higher need of soluble sugars for osmotic adjustment (Hsiao 1976;Thalmann & Santelia 2017). We detected credible effects (i.e. 95% HDI excluding zero) of both HSM and Δsugar on the mortality risk, with species with a wider HSM being at lower risk, and species with a more pronounced increase in relative leaf sugar content dying at higher rates (Fig. 3, Fig. 4, Table S2.1). This suggests that drought-induced mortality in trees is jointly influenced by both hydraulic failure and carbon starvation  Δsugar in these genera (Table 1) indicates that the reduced carbon gain caused by this strategy may have contributed to their high mortality rates.
Both access to carbohydrate reserves and their utilization rate have been reported to be controlled by water availability (Sevanto et al. 2014), providing evidence that carbohydrate use is controlled by the functionality of the water transport system (Atkin & Macherel 2009) rather

Tree size effects
In addition to the species-level effects of HSM and NSC dynamics, our model identified a positive within-species effect of tree height on the survival of all species except for Q. rubra, though the magnitude of that height effect varied strongly among species ( Fig. 4; Fig. S2.2).
This finding is in contrast to mounting evidence that larger trees tend to be more susceptible to drought-induced tree mortality across biomes (Phillips et al. 2010;Lindenmayer et al. 2012;Bennett et al. 2015;Grote et al. 2016;Olson et al. 2018;Stovall et al. 2019Stovall et al. , 2020, although some authors report opposite patterns (Van Mantgem et al. 2009;Peng et al. 2011). This apparent contradiction might result because studies reporting higher mortality for the largest trees usually focus on adult trees in mature forest stands, where crown exposure (Stovall et al. 2019) and path length-related constraints on water transport (Ryan & Yoder 1997;West et al. 1999) are more relevant. Our analysis, in turn, was based on densely spaced and young trees (planted in 2013 at the age of 1-3 years) with limited variation in tree height (Table 1)

Bark beetle effects on conifers
Individuals diagnosed with bark beetle infestation in 2018 had a much higher probability of dying (Fig. 4), consistent with a drought-driven predisposition for biotic attacks (McDowell et al. 2011). The prevalence of bark beetles differed greatly between species, with strong outbreaks only for P. abies and L. decidua (Table 1). Most likely, a decline in resin exudation and thus defence capability caused by the reduction of relative tissue water content made P.
abies vulnerable to insect infestation (Netherer et al. 2015). Despite the low overall prevalence, the infested individuals among the native conifers P. abies, L. decidua and P. sylvestris died with high probability. As bark beetles cause physical damage to the phloem and xylem that is likely to affect both carbon metabolism and xylem water transport, their role for plant death cannot be understood without acknowledging its link to other processes driving droughtinduced mortality (McDowell et al. 2008(McDowell et al. , 2011. Since both healthy and infested trees died during the drought event while a number of infested trees were able to survive, biotic attack was unlikely to be the main driver of mortality in the affected species. Therefore, pest infestations should rather be interpreted as an additional stressor amplifying the risk of death by other causes. The presence of bark beetles in neighbouring trees had no overall effect on mortality (Fig. 4), and on species level only affected P. abies, the native species that suffered the strongest beetle infestation ( Fig. S2.2). While neighbourhood effects can profoundly impact the predisposition to insect herbivory (Castagneyrol et al. 2018), the low influence in our study may be related to the timing of the assessment of pest infestation in late 2018, when the drought and associated insect outbreaks hat likely already run its course and further spread was unlikely.

Interactions with neighbour trees
Neighbourhood tree species diversity might help to mitigate drought impacts due to complementarity effects, i.e. temporal or spatial resource partitioning that results in higher   Table 1. (Fig. 5a). The number of neighbours of a particular species surrounding each focal tree varied depending on the species mixture in the plot and pre-drought background mortality. The presence of more and larger neighbours of a given species resulted in a reduced mortality risk in the majority of possible combinations. A. platanoides was the only species whose presence on average increased the mortality risk of its neighbours (Fig. 5a). Across pairs of tree species, beneficial interactions were much more common than detrimental ones (Fig. 5b). The correlation between the directed neighbourhood effects in a pair of neighbours was not credibly different from zero ( = -0.25, 95% HDI: -0.76 -0.22), indicating that an advantage for one species in a pair did not always translate into a disadvantage for its neighbour. Instead, there was a tendency of pairs of neighbour species to be either one-sided or mutually beneficial, while not a single combination resulted in an increased mortality risk for both species (Fig. S2.4). In general, species that suffered more from drought increased the survival probability of their neighbours and vice versa. While A. platanoides, the species that most strongly increased the mortality of its neighbours, was among the species least affected by drought, four of the five species with positive effects on the survival of their neighbours suffered losses of over 50% ( Fig. 1, Fig. 5a, Table 1). Besides an increased mortality of P. pungens in the neighbourhood of P. abies, the only species that reduced the survival of their neighbours were A. platanoides and Q. rubra. The presence of these two dominant overstory species increased the mortality especially for suppressed, subordinate species (Fig. 5). The increased survival of trees with neighbour species that suffered extreme losses may be a result of competitive release, especially in the case of Larix spp. and Betula spp., which did not compete for water after their early leaf shedding in 2018.
An important implication of our model is that drought-related interactions between neighbouring trees are directed and species-specific, as previously observed for adult trees preserve the information about the identity of interacting species allow for a more nuanced understanding of the underlying processes, and provide valuable additional information that e.g. permits identification of mutually beneficial species combinations. We hope that our approach to addressing neighbourhood effects contributes to a more thorough focus on the nature of directed neighbourhood interactions in plant drought responses, as in their sum they ultimately determine the resilience of a community against drought events.

Nutrient and environmental effects
Conditional on the effect of the other predictor variables, none of the nutrient-enhanced treatments (N, P or NP addition) differed in mortality compared to the control treatment for any of the analysed species (Fig. 3, Fig. S2. While no direct fertilization effects were observed, a considerable amount of the species-wise plot effects were credibly different from zero ( Fig. S2.5), indicating the presence of small-scale environmental influences on tree survival. Though the precise origin of these differences is unknown, our results indicate they had a non-negligible contribution to the total variance in mortality risk, which may be associated with plot-specific differences, e.g. in soil texture, rooting depth or wind exposition.

Conclusions
The extreme natural drought event in 2018 that killed over one third of the trees at the IDENT tree diversity experiment in Freiburg provided a unique opportunity to compare the mutually inclusive mechanisms underlying drought-induced mortality in one setting. While we support the view that severe tissue dehydration causing hydraulic failure marks the final stage in the process of tree death under drought, our results illustrate that hydraulic traits alone are not sufficient to predict individual mortality risk. Instead, the individual probability of reaching lethal desiccation thresholds depends on a series of mutually inclusive processes, notably starch depletion, known to be important for osmotic adjustment. Given the projected increase in drought exposure in most of the world's forested biomes, improving the mechanistic understanding of processes involved in drought-induced tree mortality is essential. To predict future changes in the structure and composition of forests, it is crucial to adopt a holistic view on the drivers of tree death under drought. Specifically, there is a need for an improved understanding of the processes that link carbohydrate metabolism and water relations, and the ways mortality risk is modulated by pests, environmental conditions, and neighbourhood interactions among trees varying in size and species composition in a non-zero-sum list of effects.

Tree mortality and height
We surveyed the survival and health status of all 20,335 planted trees at pre-drought (2017/18 inventory), directly after the drought in autumn 2018, and post-drought in early summer of 2019. Trees were classified as healthy, damaged (leaf browning or leaf loss), severely damaged (partial canopy dieback) or dead. Our analysis was based on the subset of the inner 25 trees of each plot that were alive during the pre-drought inventory, resulting in a total of 9,435 trees.
These trees represent 90.9% of the initially planted individuals. Since boreholes of the bark beetle Pityogenes chalcographus were observed on several trees in early July 2018, we additionally recorded the beetle infestation for all conifers.
Pre-drought tree height (stem base to the apical branch) was measured during December 2017 and January 2018 for the 25 central trees of each plot as part of a yearly inventory. Additionally, a remote sensing campaign was carried out in July 2018 using a drone (Octo8XL, Mikrokopter GmbH, Saldenburg, Germany) equipped with an optical RGB camera with a 16 mm lens (Sony A5000, Stuttgart, Germany). The heights of the edge trees (for computation of neighbourhood effects) were then imputed based on these remote-sensing derived aggregate and rescaled by the average height of the corresponding species in each plot to account for species differences (see Supplementary Material S1 for details).

Hydraulic traits
For measuring xylem embolism resistance, 100-150 cm long branch samples were collected from 12-19 live trees per species with the exception of the genus Quercus (3-5 trees × 4 plots × 10 species = 171 samples in total) in late summer 2019, immediately wrapped in moist towels and stored in black plastic bags to prevent dehydration during transport. Samples were re-cut under water directly before xylem vulnerability curves were constructed according to standard protocols for short-vesseled coniferous and diffuse-porous species using the flow-centrifuge technique (Cochard et al. 2005;Delzon et al. 2010;Schuldt et al. 2016). Vulnerability curves were fitted in R v4.0.0 (R Core Team 2020) with nonlinear least squares using a sigmoidal model based on the raw conductance measurements (Pammenter & Vander Willigen 1998;Ogle et al. 2009 Material S1 for details).

Non-structural carbohydrates
Leaf samples for the analysis of non-structural carbohydrate contents were taken during predrought (21 May to 18 June 2018) after leaf enfolding was completed in monocultures (48 plots), and post-drought (1-4 October 2018) on the full design, yielding ~1,000 sampled trees in total. Leaves were sampled from a random subset of the inner 5 × 5 trees representing the most common tree size whenever possible. Buffer trees were sampled if less than three trees per species were available for sampling in the core area. Each three to five leaves or 50-100 needles from a vital branch of the upper and middle part of the sun canopy were pooled, immediately stored at 4°C in the field and, within the same day, oven-dried at 60°C for 48 h to stop enzymatic activity (Landhäusser et al. 2018).
The low-molecular-weight sugars (glucose, fructose and sucrose) and starch were analyzed following the protocol of Wong (1990),

Neighbourhood effects
The influence of up to eight immediate neighbour trees on each of i in I central trees was expressed in form of a × neighbourhood matrix N containing the sum of the heights hj of the trees of species j of J species divided by the height hi of the central tree, weighted by their relative distance dij (i.e. with a weight of 1 for first and 1/ √2 for second order neighbours; see Supplementary Material S1).

= ∑ ℎ ℎ
In this formulation, the computed value can be interpreted as a partial Hegyi competition index for a species (Hegyi, 1974), thus effectively decomposing the contribution of different competitor species to the total neighbourhood effect.

Statistical modelling
In a hierarchical Bayesian modelling framework based on the Stan probabilistic programming language (Stan v. 2.21.0; Carpenter et al. 2017), we described the observed tree mortality status in 2019 as a function of the individual mortality risk.
We assumed that the observed state of a tree ( [ ] ; 0: tree alive, 1: tree dead) for each tree i belonging to species j in plot k for a total of I trees, J species and K plots is subject to observation errors with a species-specific probability of erroneously classifying a living tree as dead (e.g. because of leaf shedding). To account for this misclassification, we modelled the observed mortality state of individual trees as a one-inflated Bernoulli process with a treespecific probability of dying and a species-specific misclassification probability . Material S1).
We expressed the probability of dying of each tree as a logit-linear function of the true species-level Δsugar and HSM values, a term specifying species specific intercept and effects of tree height, bark beetle infestation and nutrient treatment ( × predictor matrix X), a term specifying neighbourhood effects ( × neighbourhood matrix N, elevated to a power of c), and species-specific plot effects .

logit( ) = [ ] ⋅ + Δ [ ] ⋅ Δ + + +
The × species-specific regression coefficients β were described by a multivariate normal distribution with covariance matrix Σβ: ∼ MVN( 0 , ) To account for a possible correlation between the directed neighbourhood effects within a pair of species, the × neighbourhood effects matrix γ for focal species m and neighbour species n was parameterized as follows: where is a 2 × 2 correlation matrix allowing for a correlation of the directed neighbourhood effects in pairs of non-identical species (the off-diagonal elements of γ).
The design effects δ for each plot and species were modelled by a normal distribution centred around zero with standard deviation , i.e. as a varying intercept for each combination of plot and species: ∼ Normal(0, ) Modelling was performed via R using rstan v. 2.19.3 (Stan Development Team, 2020). The model reached a ̂ statistic of under 1.002 and effective sample sizes above 1,500 for all estimated parameters, indicating full model convergence (Table S2.1). The estimated speciesspecific misclassification rates were low (on average 0.9%; Supplementary Material S1), though higher rates of up to 4.2% of trees erroneously classified as dead were found for some species (Table S2.1, Fig. S2.1), resulting in a visible mismatch between observed and predicted mortality for these species in Fig. 1

S3:
Digital supplement: A repository containing the code for data processing and fitting the model described in S1 and the full dataset is available upon request and will be made public after publication of the article under the following address: https://github.com/r-link/mutually_inclusive_mechanisms