The role of nutritional impairment in carbon-water balance of silver fir drought-induced dieback.

Rear-edge populations at the xeric distribution limit of tree species are particularly vulnerable to forest dieback triggered by drought. This is the case of silver fir (Abies alba) forests located in Southwestern Europe. While silver fir drought-induced dieback patterns have been previously explored, information on the role played by nutritional impairment is lacking despite its potential interactions with tree carbon-water balances. We performed a comparative analysis of radial growth, intrinsic water-use efficiency (iWUE), oxygen isotopes (δ18 O) and nutrient concentrations in leaves of declining (DD) and non-declining (ND) trees in silver fir in four forests in the Spanish Pyrenees. We also evaluated the relationships among dieback predisposition, intraspecific trait variation (wood density and leaf traits) and rhizosphere soil physical-chemical properties. The onset of growth decline in DD trees occurred more than two decades ago, and they subsequently showed low growth resilience against droughts. The DD trees presented consistently lower foliar concentrations of nutrients such as P, K, Cu and Ni than ND trees. The strong effects of foliar nutrient status on growth resilience indices support the key role played by mineral nutrition in tree functioning and growth before, during and after drought. In contrast, variability in wood density and leaf morphological traits, as well as soil properties, showed weak relationships with tree nutritional status and drought performance. At the low elevation, warmer sites, DD trees showed stronger climate-growth relationships and lower δ18 O than ND trees. The uncoupling between iWUE and δ18 O, together with the positive correlations between P and K leaf concentrations and δ18 O, point to deeper soil/bedrock water sources and vertical decoupling between nutrient and water uptake in DD trees. This study provides novel insights into the mechanisms driving silver fir dieback and highlights the need to incorporate tree nutrition into forest dieback studies.

INTRODUCTION

Climate change‐induced forest disturbances have increased in the past decades (Seidl et al., 2017). Global air temperatures have increased along with a rising frequency and severity of extreme climatic events such as droughts and heat waves (Christidis et al., 2015; Dai, 2013). Model‐based climate projections for the 21st century anticipate further rises in temperature and vapour‐pressure deficit accompanied by increased precipitation variability and/or reduced rainfall in some regions (IPCC, 2021; Zhou, Zhang, et al., 2019). Consequences for forest ecosystems comprise declining tree vigor and growth (Camarero, Gazol, & Sánchez‐Salguero, 2021), reduction of primary productivity (Ciais et al., 2005) and ultimately forest dieback and mortality (Allen et al., 2010, 2015; Senf et al., 2020). Moreover, projected warming rise is expected to lead to further changes in forests dynamics and structure, tending toward younger stands with faster turnover, while old‐growth forests with stable dynamics are dwindling (McDowell et al., 2020). Therefore, a better understanding about the mechanisms that allow trees to survive and adapt to hotter drought is needed to improve our projections of impending and future forest responses to ongoing climate change (Sala et al., 2010).

Tree survival under drought depends on its vulnerability to hydraulic failure, its capacity to avoid carbon starvation or the combination of both (McDowell et al., 2011). Hydraulic failure is the disruption of water transport due to massive xylem embolism that leads to plant tissues desiccation, whereas carbon starvation occurs when prolonged closure of stomata to avoid water loss causes depletion of internal carbon pools thus preventing the maintenance of metabolism (McDowell et al., 2008). Although both mechanisms are non‐mutually exclusive, hydraulic failure has received stronger support across taxa, particularly in gymnosperms (Adams et al., 2017). However, drought‐induced tree dieback and mortality are complex phenomena, and the contribution of many other factors such as biotic agents or nutritional imbalance may accelerate decline and death of already climatically stressed trees (Anderegg et al., 2015; Gessler et al., 2017; Hevia et al., 2019; Trugman et al., 2021).

The role of tree nutritional status in drought‐induced forest dieback has received little research attention despite its crucial feedbacks with tree responses to water deficit (but see González de Andrés et al., 2021; Hevia et al., 2019). Adequate levels of essential nutrients are needed for key physiological processes such as photosynthesis, hydraulic function or stomatal regulation that determine how trees cope with drought (Güsewell, 2004; Sardans et al., 2013). Furthermore, droughts have large impacts on soil nutrient availability and uptake by trees (Rozas & Sampedro, 2013). Nutrient mineralization is decreased under low soil moisture conditions due to reduction in soil microbial activity and ion mobility (Kreuzwieser & Gessler, 2010). The reduction of the active transpiration flux along the soil–plant–air continuum due to decreases in stomatal conductance would also diminish the acquisition of nutrients (Schlesinger et al., 2016). Furthermore, plant nutrient uptake from nutrient‐rich topsoil layers can be hampered by the functional impairment or death of fine roots and their symbiotic mycorrhizal fungi as a result of prolonged topsoil desiccation (León‐Sánchez et al., 2018; Querejeta et al., 2021). Either being the cause or consequence, the role played by the nutritional status in drought‐induced forest dieback requires further attention if we want to understand how trees respond to drought.

Several mechanisms have been proposed to operate complex interplays between tree nutrition and drought. For instance, Salazar‐Tortosa et al. (2018) reported the existence of a detrimental feedback loop between reduced transpiration and acquisition of nutrients in a field experiment with European pine saplings. Besides, intraspecific and seasonal differences have been detected in several tree species that are able to shift their water uptake depth as a response to drying of the shallow soil layers (e.g., Brinkmann et al., 2019; Ripullone et al., 2020; Volkmann et al., 2016; Voltas et al., 2015; Zhou, Zhao, et al., 2019). Indeed, experimental studies and recent modelling approaches have revealed the widespread usage of deep soil/bedrock moisture for water uptake under high evaporative demand (McCormick et al., 2021; Querejeta et al., 2007; Stocker et al., 2021). However, water availability in trees relies mainly on deep soil layers or weathered bedrock during prolonged drought, which may induce nutrient deficiencies, as many essential nutrients decrease as soil depth increases (Jobbágy & Jackson, 2001; Querejeta et al., 2021). Subsequent nutritional impairment can exacerbate and amplify the negative effects of drought on tree physiological functioning and productivity and increase the risk of tree mortality (Gessler et al., 2017; González de Andrés et al., 2021; Hevia et al., 2019).

Considering tree growth as a sensitive indicator of the trees carbon‐water balance, as well as their nutritional status, the retrospective analysis of tree growth trends, assessed by the study of tree rings, provides reliable way to investigate the above‐mentioned drought‐induced unbalances. Common patterns, although more common in gymnosperms, have been identified as early warning signals of tree dieback and mortality including low long‐term growth rates, high inter‐annual growth variability or low growth resilience to drought events (Cailleret et al., 2017; Camarero et al., 2015, 2018; DeSoto et al., 2020; Serra‐Maluquer et al., 2021a). However, the specific individual traits and environmental factors that predispose trees to growth decline and death still remain unclear (Kannenberg et al., 2020). Further, the analysis of carbon and oxygen isotope ratios in tree tissues can help to disentangle the factors driving forest decline and mortality (Billings et al., 2016; Cherubini et al., 2021; Gessler et al., 2018). Stable carbon isotope ratios (δ13C) of plant tissues are determined by the relative importance of diffusive limitation and stomatal conductance versus biochemical limitations of photosynthesis, which in turn allow for the estimation of time‐integrated intrinsic water‐use efficiency (iWUE) (Farquhar et al., 1982). The analysis of oxygen isotopes (δ18O) can further help discern the effects of environmental changes on stomatal conductance, since δ18O is not substantially affected by photosynthetic rates, but it is influenced by leaf‐level evaporative effects and post‐photosynthetic isotope fractionation processes (Gessler et al., 2014; Scheidegger et al., 2000). However, caution is needed when interpreting the so‐called dual‐isotope approach since plant δ18O may not only reflect evaporative processes at the leaf level but also variations in the δ18O of the source water used by the tree related to inter‐annual variability and changing water uptake depths (Barbour, 2007; Roden & Siegwolf, 2013). In fact, Treydte et al. (2014) found stronger signals of source water isotopic composition than leaf‐level evaporative enrichment on the δ18O signature of Larix decidua tree rings, which is in line with the predominant effect of soil water δ18O on the isotopic composition of tree rings reported by previous studies (e.g., Barbeta & Peñuelas, 2017; González de Andrés et al., 2021; Sarris et al., 2013).

Last but not least, the dual isotope approach exclusively focuses on carbon and water relations in trees, albeit overlooking the potentially relevant contribution of nutrient availability and uptake to dieback processes (González de Andrés et al., 2021; Hevia et al., 2019; Houle et al., 2007). Growth responses to drought and mineral nutrition of trees can co‐vary according to functional trait variability, microsite and soil conditions and competitive pressure. Wood density and leaf morphological traits (e.g., leaf area and leaf mass per area) have been used to described strategies in resource acquisition and mortality risk along environmental gradients (Chave et al., 2009; Greenwood et al., 2017; Wright et al., 2004), although inconclusive results have been reported at intraspecific levels (Anderegg et al., 2018; Fajardo, 2016; Serra‐Maluquer et al., 2021b). Belowground processes including soil nutrient dynamics can both impact and be impacted in complex ways during forest dieback (Gazol et al., 2018a; Rozas & Sampedro, 2013; Štursová et al., 2014). Therefore, the combination of all these different complementary approaches could provide valuable insights into the multiple links between nutrient‐related processes and drought‐induced forest dieback.

Here, we compare the radial growth patterns and foliar isotopic and nutrient composition of coexisting declining and non‐declining trees of silver fir (Abies alba Mill.) forests of the Spanish Pyrenees. We also compare those patterns with functional traits (wood density, leaf traits) and the soil nutritional status of sampled trees. Silver fir is a keystone montane and subalpine conifer widely distributed across Europe that reaches its southern and xeric distribution limit in the Spanish Pyrenees (Caudullo et al., 2017). Given that the tree growth sensitivity to climate is expected to increase under climate change scenarios at these rear‐edge populations (Gazol et al., 2020; Sánchez‐Salguero et al., 2017), advancing the understanding of drought‐induced dieback of silver fir populations near its rear edge is fundamental for understanding how it will respond to further warming and drying climate trends in other European regions. In this study, we aim at (i) analyzing the long‐ and short‐term growth trends of declining and non‐declining silver fir trees showing contrasting vigor (estimated as different crown defoliation); (ii) characterizing their foliar carbon and oxygen isotope ratios and nutrient composition as a time‐integrated proxy of water and carbon balance and nutritional status (iii) and relating growth, isotopic and elemental composition data to each other and to intraspecific trait variability and soil properties in order to outline mechanistic linkages among them. We hypothesized that nutritional impairment in declining trees amplifies the negative impacts of recurrent hotter droughts on tree performance and thus plays a key but overlooked role in drought‐induced forest dieback processes in rear‐edge silver fir forests.

MATERIALS AND METHODS

Study sites and climatic data

The study sites were montane silver fir (Abies alba Mill.) stands located in north‐western Aragón, in the central‐western Spanish Pyrenees (Figure S1, Table 1). The populations inhabiting these marginal forests represent the southwestern limit of the distribution of the species in Europe. Here, silver fir usually grows in mesic sites on north‐facing slopes, forming pure or mixed forests with European beech (Fagus sylvatica L.) or Scots pine (Pinus sylvestris L.) and the understory vegetation is commonly dominated by several shrubs such as European box (Buxus sempervirens L.).

TABLE 1

Characteristics (mean ± standard error) of sampled sites and silver fir trees for non‐declining (ND) and declining (DD) trees

	Cotatuero (CO)		Paco Ezpela (PE)		Paco Mayor (PM)		Selva de Oza (SO)
	ND	DD	ND	DD	ND	DD	ND
Latitude N (º)	42º39′10″		42°44′29″		42º42′36″		42º49′45″
Longitude W (º)	00º02′41″		0°49′37″		0º 38′51″		0º42′15″
Elevation (m a.s.l.)	1530		1075		1320		1285
Aspect	S – SW		N – NE		N – NE
Temperature (ºC)	7.19 ± 1.06		9.39 ± 0.80		9.15 ± 0.93		8.25 ± 0.92
Precipitation (mm)	1229 ± 230		995 ± 234		1150 ± 253		1411 ± 272
DBH (cm)	38.79 ± 1.22	45.63 ± 2.92	32.31 ± 1.84	31.73 ± 1.67	30.65 ± 1.39	33.54 ± 1.76	49.57 ± 4.59
Tree age (y)	118.3 ± 4.4	120.4 ± 4.2	98.9 ± 6.8	114.8 ± 4.9	89.3 ± 5.3	94.5 ± 6.6	103.5 ± 12.2
DCI	1.01 ± 0.08	0.85 ± 0.10	0.94 ± 0.08	1.22 ± 0.09	0.84 ± 0.09	1.00 ± 0.10	0.77 ± 0.07
No. sampled trees (No. radii)	15 (28)	16 (29)	20 (37)	19 (33)	19 (34)	19 (40)	20 (40)
EPS > 0.85 since	1892	1900	1905	1914	1912	1919	1926
TRW*	1.18 ± 0.07	1.43 ± 0.16	1.66 ± 0.12	1.47 ± 0.07	1.96 ± 0.21	1.86 ± 0.18	2.78 ± 0.17
Rbar *	0.51 ± 0.02	0.48 ± 0.05	0.49 ± 0.03	0.48 ± 0.05	0.58 ± 0.03	0.55 ± 0.04	0.53 ± 0.03
AC*	0.81 ± 0.02	0.82 ± 0.02	0.72 ± 0.03	0.76 ± 0.03	0.75 ± 0.04	0.77 ± 0.03	0.74 ± 0.04
MS*	0.17	0.21	0.24	0.25	0.25	0.24	0.19

Variables’ abbreviations: diameter at breast height (DBH), distance‐dependent competition index (DCI), year from which expressed population signal (EPS) is higher than 0.85, tree‐ring width (TRW), mean inter‐series correlation (Rbar), first‐order autocorrelation (AC) and mean sensitivity (MS).

Calculated for the period 1950–2019 on raw (TRW, AC) or standardized (Rbar, MS) ring‐width values.

We selected four different silver fir stands, three of them characterized by abundant dominant trees with high levels of crown defoliation (Cotatuero, CO; Paco Ezpela, PE; and Paco Mayor, PM), and the other site inhabited only by healthy fir trees (Selva de Oza, SO). Historical management also differed among study forests. PE, PM and SO had signs of intense logging activity in the past with abundant stumps and wood trails (Cabrera, 2001; Camarero et al., 2011). Meanwhile, CO had not been logged for at least 50 years since it is located in the “Ordesa and Monte Perdido” National Park, where strict conservation policies have been implemented since the establishment of the park in 1918.

Study stands are located on areas with marls and limestones lithologies and the climate in the study area is continental with relatively cool and wet summers. Long‐term records of climatic data of each study site were obtained from the 1‐km2 gridded E‐OBS v. 22.0e database (Cornes et al., 2018). The remaining elevation distance between grids’ average and study sites was accounted for by assuming an altitudinal lapse rate of air temperature of −0.0055°C m−1. The mean annual temperature across the study sites is 8.2°C ranging from 7.2°C in the coldest site (CO, 1530 m a.s.l.) to 9.4°C in the warmest site (PE, 1075 m a.s.l.). Annual precipitation ranges from 995 mm in PE to 1411 mm in SO, the wettest site. To characterize drought severity, 1.1‐km2 resolution series of the standardized precipitation‐evapotranspiration index (SPEI) for the period 1961–2015 were used (Vicente‐Serrano et al., 2017). This is a standardized multi‐scalar drought index based on the accumulated water deficit, in which negative values indicate a negative cumulative water balance (Vicente‐Serrano et al., 2010).

Field sampling and dendrochronological methods

In total, 128 dominant and codominant, mature silver fir trees were selected for the present study. At sites with canopy dieback (CO, PE and PM), we sampled pairs of non‐declining (ND) and declining (DD) trees (15–20 pairs per site totaling 108 trees). At sites without canopy dieback (SO) only non‐declining trees were sampled (20 trees). Tree vigor classes were established based on defoliation by visual assessment of crown transparency (Dobbertin, 2005). We selected the trees based on their recent defoliation levels by considering neighboring individuals showing high (>60%, DD trees) or low (<40%, ND trees) canopy defoliation, respectively. This defoliation threshold was selected based on previous research on this species (Camarero et al., 2011, 2015; Hevia et al., 2019). At the healthy site (SO), mature trees were randomly selected and sampled. For each sampled tree, diameter at breast height (DBH) and height of the top of the tree were measured using tapes and a laser rangefinder (Nikon Forestry Pro II), respectively. In addition, a sun‐exposed branch was cut using a telescoping pole and fully developed and undamaged leaves from the last cohort were placed in sealed plastic bags and transported to the laboratory. Field sampling was conducted in autumn 2019 and summer 2020.

We characterized the neighborhood of every focal tree by measuring distance and DBH and annotating species identity of the nearest neighbor in each cardinal point from the focal tree. The degree of competition around each sampled tree was assessed using a distance‐dependent competition index (DCI; Hegyi, 1974). The index was calculated as follows: where DBHi is the DBH of the focal tree, DBHj is the DBH of the neighbor tree, and distij is the distance between both trees. This index is derived from the hypothesis that the competitive effect of a neighbor tree increases with increasing size and proximity (Tomé & Burkhart, 1989).

Radial growth of sampled trees was evaluated using dendrochronology. Two cores at 1.3 m height were extracted from each tree, perpendicular to the maximum slope, using 5 mm Pressler increment borers (Haglöf, Sweden). The wood samples were air‐dried, glued onto wooden mounts and polished until the xylem cellular structure was visible (Fritts, 1976). All samples were visually cross‐dated, and the tree‐ring width (RW) was measured to a precision of 0.01 mm using a LINTAB measuring device (Frank Rinn). Cross‐dating was further validated using the COFECHA software, which calculates moving correlations among individual tree series (Holmes, 1983). Tree‐ring width series were transformed to basal area increment (BAI) series because it is a two‐ dimensional measure of stem increment in area that is known to better reflect growth of the whole tree than the one‐dimensional ring width (Biondi & Qeadan, 2008). BAI series were calculated using the following equation and assuming concentric rings: where R 2 and R 2 −1 are the radii corresponding to the years t and t − 1, respectively.

To detrend each individual tree‐ring width series, we applied a cubic regression spline with a 50% frequency response cutoff at 30 years to the raw ring‐width series (RW). Afterwards, an autoregressive model was applied to each detrended series to remove the first‐order autocorrelation therefore building residual, pre‐whitened ring‐width index chronologies. To obtain residual site chronologies (RWIres), individual series were averaged year‐by‐year using a biweight robust mean. All the processes of tree‐ring series detrending and chronology computation were performed using the package dplR (Bunn et al., 2020) in R software (version 4.1.1., R core Team, 2021).

Drought event selection and growth performance

We used bootstrapped correlations to identify the time scale at which drought (SPEI) impacted most on radial growth. Correlations were performed among site RWIres chronologies and SPEI at the 1‐, 3‐, 6‐, 9‐, 12‐ and 24 time scale for the months from previous September to September of the year of tree‐ring formation. At every study site, the strongest growth response was found to August SPEI calculated at 3‐month time (SPEI.3. Aug) (Figure S2), which was consistent with previous analyses (Pasho et al., 2011). Relationships with climate were assessed using the R‐package treeclim (Zang & Biondi, 2015). Then, values of SPEI.3Aug below the 90th percentile during the 21st century were identified as drought years because they were the most severe and recent dry spells affecting the study forests. We selected drought events only during the last two decades to minimize chances of morphological and chemical changes in the sampled trees. Selected drought events were 2005 and 2012 (Figure S3).

Growth response to extreme droughts was measured as resistance (Rt), recovery (Rc) and resilience (Rs) indices proposed by Lloret et al. (2011). The indices based on the ratios of pre‐drought, drought and post‐drought growth BAI values were calculated as follows: where BAID is the growth during the drought event and BAIpreD and BAIostD are the mean BAI of the three years preceding and following, respectively, the drought when no drought conditions occurred. We chose consistent preD and postD periods of three years, as previous studies have shown that they represent a good compromise between drought intensity and the short‐term growth response in silver fir (Gazol et al., 2017). Values of Rt, Rc and Rs of 2005 and 2012 were averaged to a single value per tree, thus reflecting the tree performance across the two droughts.

Nutrient concentrations and isotopic composition of leaves

We used leaf C and O isotopic composition and nutrient concentrations as time‐integrated proxy measures of carbon and water relations and nutrient status of trees, respectively. Leaves were oven‐dried 72 h at 70°C, milled and homogenized to a fine powder using a ball mill (Retsch ZM1, Haan, Germany). Leaf phosphorus (P), potassium (K), aluminum (Al), calcium (Ca), copper (Cu), magnesium (Mg), manganese (Mn), nickel (Ni), silicon (Si) and strontium (Sr) concentrations were measured by inductively coupled plasma optical emission spectrometry (ICP−OES; Thermo Elemental Iris Intrepid II XDL) after a microwave‐assisted digestion with HNO2:H2O2 (4: 1, v/v). In addition, we calculated two elemental ratios (Mn:Ca and Mn:Al) on a mass basis that have been previously proposed as early warning signals of tree dieback (Hevia et al., 2019; Houle et al., 2007).

The carbon isotopic composition (δ13C) and C and N concentrations of leaves were determined by continuous flow dual isotope analysis using a CHNOS elemental analyzer interfaced to an IsoPrime100 mass spectrometer. Leaf δ13C isotope composition is expressed in delta notation (‰) relative to the Vienna Pee Dee Belemnite standard (V‐PDB). The leaf oxygen isotope composition (δ18O) was determined in continuous flow using an Elementar PYRO Cube interfaced to a Thermo Delta V mass spectrometer. Leaf δ18O isotope composition is expressed in delta notation (‰) relative to VSMOW for δ18O. Long‐term external precision values for leaf δ13C and δ18O determinations are ±0.10‰. and ±0.20‰, respectively. All stable isotope analyses were conducted at the Center for Stable Isotope Biogeochemistry, University of California‐Berkeley (USA).

The foliar δ13C was used to estimated intrinsic water‐use efficiency (iWUE) as the ratio between the photosynthetic rate (A) and the stomatal conductance rate (g), following Farquhar et al. (1982): where Ca and Ci are CO2 concentrations in the atmosphere and the intercellular space, respectively, and 0.625 is the relation among conductance of H2O and CO2. To determine Ci, we used the next equation: where δ13Ctree and δ13Catm are the tree and atmospheric C isotope compositions, respectively, a is the diffusion fractionation across the boundary layer and the stomata (+4.4‰) and b is the Rubisco enzymatic biologic fractionation (+27.0‰).

Functional trait measurements

Trait measurements were conducted following standard protocols (Pérez‐Harguindeguy et al., 2013). The leaf area (LA) of 20 needles per tree including both ND and DD trees were scanned and measured with ImageJ software (Schneider et al., 2012) while being still fresh. Afterwards, sampled leaves were oven‐dried during 72 h at 70 ºC prior obtaining dry weights. Then, leaf mass per area (LMA), which is a measure of biomass investment in leaves per unit light interception and gas exchange area (Poorter et al., 2009), was calculated by dividing leaf dry weight by LA. In addition, we measured the stem wood density (measured as wood specific gravity, WSG). WSG describes the carbon investment in woody tissues and has important implications in mechanical support and water transport (Chave et al., 2009). Stem cores (10 mm diameter) were extracted from each tree and taken to the laboratory, where fresh volume of the first 5 cm of sapwood was estimated by the volume replacement method after bark removal. Cores were then oven‐dried during 72 h at 70°C and dry weight was measured. WSG was estimated as the ratio between dry weight and fresh volume.

Soil sampling and physico‐chemical analyses

We collected three soil samples below the canopy projection of each tree using a 5‐cm diameter soil borer from the uppermost 15 cm where most fine roots are found (Janssens et al., 2002). Litter was removed before soil sample collection. Soil subsamples were mixed in sealed plastic bags to create a unique composed sample per tree. Due to field work conditions, we did not collect soil samples beneath ND trees at PM site. Samples were air‐dried on a glasshouse and sieved with a 2‐mm mesh size. Soil C (organic and total), N and P were determined with an elemental analyzer (Element Analyzer VarioMAX N/CM). Soil texture was determined with a laser diffraction method in a particle analyzer (Coulter Mastersizer 2000), and clay content was corrected following Taubner et al. (2009). To integrate the different components of soil texture into one single variable, the exponent of the Saxton equation (Saxton et al., 1986) was calculated as follows: where less negative values of b indicate sandy soils with lower soil water retention capacity.

Statistical analyses

We used the Kruskal Wallis test to evaluate differences between ND and DD trees (overall and at site level) regarding foliar C and O isotopic composition and nutrient concentrations, functional traits (height, LA, LMA and WSG) and soil characteristics. The Wilcoxon rank‐sum test was employed to check if the changes through time of BAI differed between tree‐vigor classes at each study site with canopy dieback.

We used generalized additive mixed models (GAMM; Wood, 2017) to describe silver fir long‐term growth trends and individual responses to climate. For each site, we adjusted individual BAI as a function of a smooth function of calendar year using thin‐plate regression splines with a maximum of five degrees of freedom (Wood, 2003). Different growth trends between ND and DD trees were assessed by including an interaction between the smooth function of calendar year and vigor class. Models were compared with and without the interaction terms. To account for growth responses to climate, we included as linear predictors the SPEI.3Aug and its interaction with vigor class. We also included tree identity as a random effect, since BAI series are repeated measures over the same individual, and a first‐order autocorrelation structure. The GAMMs were fitted using the R‐package mgcv (Wood, 2011).

To obtain a multidimensional overview of the nutritional status of trees at site and regional (i.e., considering all sites together) scales and to compare the two vigor classes (ND vs. DD trees), we performed non‐metric multidimensional scaling (NMDS) with Euclidean dissimilarity for nutrient concentrations (N, Ca, Cu, K, Mn, Ni, P, Si, Sr) (Legendre & Legendre, 2012). Significant differences between ND and DD trees were tested with permutational multivariate analysis of variance (PERMANOVA; Anderson, 2001). The NMDS axes were correlated with iWUE and δ18O and projected into the ordination space. Multivariate analysis was performed with the R‐package vegan (Oksanen et al., 2019). Pairwise Spearman rank correlations were used to test associations between isotopic composition (iWUE and δ18O) and individual nutrient concentrations of leaves. Covariations between tree nutrient status (tree scores along NMDS1 and NMDS2), intraspecific trait variation, competition pressure (DCI) and soil conditions were also evaluated using Spearman correlations. For the sake of clarity, correlations were calculated using tree scores along NMDS axes of the analysis performed at regional scale, i.e., considering all sites together.

We used linear mixed‐effects models (LMMs; Pinheiro & Bates, 2000) to determine which factors influenced short‐term growth responses to drought. One model was adjusted for each resilience index (Rt, Rc and Rs). Predictor variables were nutritional status (scores of the NMDS1 and NMDS2 of ordination of foliar nutrient concentrations), foliar isotopic composition (iWUE and δ18O), leaf and stem traits (tree height, LA, LMA and WSG), soil physico‐chemical properties (Saxton b coefficient, pH, N and P assimilable concentrations and C:N ratio) and their interaction with tree vigor class. Predictors were standardized by subtracting the mean and dividing by the standard deviation prior to the analyses to obtain comparable parameters across covariates. Growth response indices were log‐transformed to achieve normality assumptions. LMMs also included study site as random effect (Zuur et al., 2009). All combinations of predictors in the saturated model were fitted, and the model with the lowest corrected Akaike information criterion (AICc) was selected (Burnham & Anderson, 2002). The goodness of fit of models was evaluated with the coefficient of determination for GLMMs (R_GLMM 2) proposed by Nakagawa et al. (2017). Marginal R 2 (R_GLMM(m) 2) accounts for the proportion of variance explained by the fixed effects, and conditional R2 (R_GLMM(c) 2) is interpreted as the variance explained by the entire model, that is, fixed plus random effects. Estimated marginal means of linear trends were computed for significant interactions, and pairwise comparisons were calculated to test for significant differences among vigor classes. The R‐packages nlme (Pinheiro et al., 2020) and MuMIn (Barton, 2019) were used to fit models and calculate goodness of fit, respectively, whereas emmeans (Lenth, 2020) was used to test for significant differences among vigor classes.

RESULTS

Long‐term radial growth and drought performance

BAI of DD trees was significantly lower than that of ND trees since 2011 in PE, 2015 in CO and 2016 in PM. Meanwhile, DD trees in CO grew more than ND trees during the period 1926–1974 (Figure 1). According to GAMMs, the effect of year differed between vigor classes (Table 2; Table S1). The negative effect of year for about two decades in DD trees (Figure 2a–d) pointed to the growth decline of this vigor class since late 20th–early 21st centuries. Although no canopy dieback was observed at the SO site, tree growth also showed a negative trend in the last few years.

FIGURE 1

Interannual variation of basal area increment (BAI) of non‐declining (solid green lines) and declining (dashed red lines) silver fir trees of the four study sites. Solid lines represent the means and shaded areas around them the standard error of the mean. The grey filled areas indicate the periods when BAI of tree vigor classes significantly (p < .05) differed according to Wilcoxon rank‐sum tests. Dash vertical lines represent drought years for which growth response was evaluated using resilience indices

TABLE 2

Generalized additive mixed effect models (GAMMs) selected to explain growth trends (basal area increment, BAI) in silver fir at each study site

	Cotatuero		Paco Ezpela		Paco Mayor		Selva de Oza
	β ± SE	t	β ± SE	t	β ± SE	t	β ± SE	t
Class (ND)	−3.09 ± 2.07	−1.49	1.84 ± 1.30	1.42	−1.15 ± 1.28	−0.90	−	−
SPEI	0.19 ± 0.07	2.94*	0.38 ± 0.06	6.19*	0.39 ± 0.07	5.53*	1.09 ± 0.13	8.15*
Class (ND) × SPEI	−0.03 ± 0.08	−0.34	−0.61 ± 0.08	−7.11*	−0.38 ± 0.11	3.62*	−	−
	Edf	F	Edf	F	Edf	F	Edf	F
ND × Year	3.63	40.26*	3.60	7.14*	3.85	35.93*	2.55	2.03*
DD × Year	3.23	56.38*	3.65	22.53*	3.73	43.46*
R 2	0.194		0.378		0.444		0.086

For the variables included as linear terms in the models, the regression coefficient (β ± SE) and the t statistics are shown. For the variables included as spline functions in the model the degrees of freedom (Edf) and the F statistics are shown. Class refers to non‐declining (ND) and declining (DD) trees.

Significant terms (p < .05).

FIGURE 2

Long‐ and short‐term growth trends of non‐declining (green solid lines and empty bars) and declining (red dashed lines and striped bars) silver fir trees. (a–d) Effect of calendar year in the growth trends according to GAMMs. Solid lines represent the means and shaded areas around them the standard error of the mean. (e–h) Resistance (Rt), recovery (Rc) and resilience (Rs) against the driest spells of the 21st century (2005 and 2012). Significant differences (p < .05) between non‐declining and declining trees according to Kruskal Wallis test are indicated with different letters

Our results indicated sensitivity of the radial growth of silver fir to summer drought since SPEI.3Aug showed the strongest correlations with site‐level RWIres (Figure S2). We found different growth responses to water deficit between vigor classes in PE and PM but not in CO (Table 2). In PE and PM, DD individuals were more sensitive to water deficit than ND trees. Tree growth responses to the most severe droughts of the 21st century differed between sites and vigor classes (Figure 2e–h). ND trees showed better drought performance than DD trees in terms of Rc (CO and PM) and Rs (CO, PE and PM). No site‐level differences between vigor classes were found regarding Rt.

Patterns in foliar composition, functional traits and soil properties

Regional comparison revealed higher foliar δ13C and iWUE in DD trees compared to ND trees, pointing to higher drought stress of the former vigor class, although this was strongly driven by the low iWUE of ND trees growing at the moist SO site. However, this result does not agree with the within‐site patterns where significant differences were only found in PE in the opposite direction, that is, higher foliar δ13C and iWUE in ND than in DD trees (Figure 3). In the case of the foliar δ18O composition, we found a general pattern of lower δ18O values in DD trees across sites that was consistent with within‐site results as similar significant differences were found in PE and PM. Both foliar isotopic signatures (δ13C/iWUE and δ18O) were not significantly correlated with each other (CO, PM and SO) or were negatively correlated (PE).

FIGURE 3

Relationships between foliar oxygen isotopic ratio (δ18O) and intrinsic water‐use efficiency (iWUE) of silver fir trees in each study site. Shaded areas represent the centroid of the values for each tree‐vigor class. The results of the Kruskal–Wallis test and associated significance (*p < .05) evaluating differences between non‐declining (green dots) and declining trees (red triangles) regarding iWUE and δ18O are indicted in the top right of each graph. Significant associations between isotopic signatures are represented with solid black lines showing Spearman ρ statistic and associated p‐value. Note that x and y axes have the same scale for allowing comparison between graphs

Regarding nutrient concentrations in leaves, differences between vigor classes were consistent across scales (Table 3; Figure S4). Leaf concentration of P, K, Cu and Ni were higher in ND trees than in DD trees, whereas the opposite situation occurred with concentration of Ca, Mg, Si and Sr, which were higher in DD trees. Although, Mn concentration was not significantly different when all sites were considered together, it showed higher values in DD trees as compared to ND trees in PE. Ordination analysis helped us to characterize the nutritional status of trees using the foliar elemental composition as a time‐integrated proxy. The first axis (NMDS1) was negatively related to P, K and Cu concentrations and positively correlated to Ca, Mn, Si and Sr concentrations. The second axis (NMDS2) encompassed most of the variability of foliar Al and Ni concentrations, which showed a negative association with N and Mg concentration in leaves (Figure 4). Site‐level ordinations produced similar results (Figure S5). Vigor classes were significantly different in their overall nutritional status at both regional and site scales based on PERMANOVA test, and these results agree with the individual nutrient comparisons. Global differences in the Mn:Al ratio, which was higher in DD trees than in ND trees, remained significant at the site‐level. Meanwhile, no significant differences between vigor classes were found regarding the Mn:Ca ratio.

TABLE 3

Mean (± standard errors) isotopic composition and nutrient concentrations in leaves, and functional traits of silver fir trees considering non‐declining and declining trees for all study sites

		Non‐declining trees	Declining trees
Isotopic composition	δ13C (‰)	−30.29 ± 0.18a	−29.74 ± 0.15b
	δ18O (‰)	32.68 ± 0.15a	32.08 ± 0.21b
	iWUE (µmol mol−1)	52.24 ± 2.18a	58.81 ± 1.77b
Nutrient concentrations	N (%)	1.27 ± 0.02a	1.25 ± 0.03a
	P (mg g−1)	1.48 ± 0.06a	0.96 ± 0.04b
	K (mg g−1)	8.61 ± 0.26a	5.82 ± 0.30b
	Al (µg g−1)	81.98 ± 7.49a	57.49 ± 4.63a
	Ca (mg g−1)	4.74 ± 0.17a	7.28 ± 0.40b
	Cu (µg g−1)	4.13 ± 0.21a	3.07 ± 0.21b
	Mg (µg g−1)	66.53 ± 1.70a	72.73 ± 2.69b
	Mn (µg g−1)	261.82 ± 21.01a	347.68 ± 41.82a
	Ni (µg g−1)	2.26 ± 0.23a	0.85 ± 0.09b
	Si (µg g−1)	88.09 ± 7.19a	149.36 ± 10.52b
	Sr (µg g−1)	7.47 ± 0.58a	11.71 ± 1.07b
	Mn:Ca	0.56 ± 0.04a	0.45 ± 0.04a
	Mn:Al	4.02 ± 0.24a	5.98 ± 0.50b
Functional traits	Height (m)	20.57 ± 0.58a	18.68 ± 0.45b
	Leaf area (cm2)	598.74 ± 16.51a	429.26 ± 13.84b
	LMA (mg mm−2)	0.12 ± 0.01a	0.15 ± 0.01b
	WSG (g cm−3)	4.3 10−3 ± 1.6 10−4a	6.2 10−3 ± 3.1 10−4b
Soil properties	Saxton b	−4.89 ± 0.06a	−5.02 ± 0.08a
	pH	6.17 ± 0.09a	6.76 ± 0.15b
	N (%)	0.41 ± 0.03a	0.49 ± 0.05a
	C:N	17.10 ± 0.32a	17.20 ± 0.28a
	P assimilable (ppm)	25.41 ± 2.82a	29.98 ± 4.41b

Different letters indicate significant differences (p < .05) between tree vigor classes according to the Kruskal Wallis test.

FIGURE 4

Non‐metric multidimensional scaling (NMDS) biplot of nutrient concentrations in leaves considering all study sites together of non‐declining (green dots) and declining trees (red triangles). Shaded areas represent the centroid of the values for each tree‐vigor class (green area, non‐declining trees; red area, declining trees). Black arrows indicate loadings of nutrient concentrations. The correlation between intrinsic water‐use efficiency (iWUE) and oxygen isotope composition (δ18O) and NMDS axes was projected in the ordination diagram and represented with blue arrows. F statistic and associated p‐value of the PERMANOVA test comparing tree vigor classes is shown

Both δ18O and iWUE were negatively correlated with tree scores along NMDS1 axis considering all sites together (Figure 4). These relationships remained significant at site scale for δ18O in PE, PM and SO, whereas iWUE showed significant correlation with NMDS1 only in PE (Figure S5). There were no significant correlations between isotopic signals in leaves and tree scores along NMDS axes at the CO site. Pairwise Spearman rank correlations between foliar isotope composition and individual nutrient concentrations further supported these results (Table S3).

At the regional scale, ND trees showed consistently higher LA and lower WSG values than DD trees (Table 3), and such global differences remained significant at the site level (Figure 5b,d). ND trees were taller than DD trees when all sites were considered together (Table 3). However, this pattern was not noticeable at the site‐level; therefore, it was likely mediated by the higher height of trees in SO (Table S2; Figure 5a). A similar situation occurs regarding LMA, for which the general discrepancy between vigor classes was strongly driven by the large difference encountered in PM (Figure 5c).

FIGURE 5

Variation in functional traits between sites and vigor: non‐declining (ND; green empty bars) and declining trees (DD; red striped bars) of studied silver fir forests including: (a) height, (b) leaf area (LA), (c) leaf mass per area (LMA), and (d) wood specific gravity (WSG). Error bars represent standard errors. Different letters on the top of each bar indicate significant differences (p < .05) between vigor classes of the same study site based on the Kruskal–Wallis test

Comparisons of soil physico‐chemical properties beneath the canopy of different vigor classes in CO and PE showed significant differences only for C:N ratio (CO and PE) and assimilable P (PE). Additionally, the high concentrations of N and P in the rhizosphere soil that were found at CO in comparison with the other study sites are noteworthy (Figure S6).

Tree nutritional status was coordinated with the intraspecific trait variation revealed by the correlations that tree scores along NMDS1 axis showed with WSG (negative correlation at all study sites), LMA (negative correlation at CO and PM) and LA (positive correlation at CO, PE and PM) (Figure S7). That is, trees with higher wood density and smaller and more xerophyll leaves also exhibit lower P, K and Cu and higher Ca, Si and Sr concentrations in leaves, indicating that lower wood hydraulic conductivity and smaller and more xerophyll leaves are coupled to impaired leaf nutrient status in silver fir. Indeed, LA and WSG were also strongly negatively correlated with each other at all the study sites. Besides, a higher tree‐to‐tree competition pressure (DCI) was also associated with poorer leaf nutritional status at three of the four study sites (PE, PM and SO). Regarding rhizosphere soil properties, no consistent relationships with the tree nutritional status were found except for a positive association between tree scores along NMDS1 axis and the soil concentration of assimilable P at the sites showing dieback (Figure S7), which paradoxically indicates higher soil P availability in the rhizosphere of P‐deficient declining trees.

Impacts of tree characteristics and soil properties on drought growth performance

The tree nutritional status had a strong impact on drought performance of tree growth as we found significant effects of the scores of the first and second ordination axes of the foliar nutrient concentrations on Rt, Rc and Rs (Table 4). Trees with lower NMDS1 values (i.e., higher foliar P, K, Cu and lower of Ca, Mn, Sr and Si concentrations) showed an overall better response against droughts (Figure 6). The tree scores along NMDS2 also had a significant impact on Rt and Rs (only in ND trees in the latter), with positive effects of increasing Al and Ni concentrations on drought performance (Table S4). Tree height had a negative effect on every drought performance index, whereas bigger leaves enhanced Rt and Rs in ND trees. The most consistent effect of soil properties was exerted by pH, which showed a positive relationship with growth performance under drought. In addition, DD trees showed poorer post‐drought recovery capacity at microsites with higher C:N ratios in soil. Finally, we found no significant effect of foliar isotopic composition or competitive pressure on short‐term responses to recent droughts in the study silver fir forests.

TABLE 4

Selected linear mixed‐effects models characterizing the resistance, recovery and resilience growth indices against droughts

	Resistance	Recovery	Resilience
Vigor class	4.129*	3.058+	9.540**
NMDS1	16.174***	4.154*	8.583**
NMDS2	7.272**		24.117***
Height	3.718+	10.521**	5.089*
LA	0.337	0.063	1.008
pH	4.853*	7.786**	4.798*
C:N		5.536*
Vigor class × NMDS2			4.813*
Vigor class × LA	6.199*		4.346*
Vigor class × C:N		9.488**
R_GLMM(m) 2	0.263	0.235	0.337
R_GLMM(c) 2	0.279	0.251	0.349

Abbreviations: LA, leaf area. pH and C:N ratio refer to soil variables; NMDS1 and NMDS2, scores of first and second axis of ordination of foliar nutrient concentrations.

For each variable, the F statistic and the associated probability (+ p < .1; *p < .05; **p < .01; ***p < .001) are shown. R_GLMM(m)2 (marginal coefficient of determination, i.e. proportion of variance explained by the fixed factors) and R_GLMM(c)2 (conditional coefficient of determination, i.e. proportion of variance explained by the entire model, including both fixed and random effects; Nakagawa et al., 2017) are presented for each model.

FIGURE 6

Relationship between tree nutritional status (tree scores along the first axis of non‐metric multidimensional scaling of foliar nutrient concentrations; NMDS1) and tree growth drought performance (a, Rt (resistance index); b, Rc (recovery index); c, Rs (resilience index)) of non‐declining (green dots) and declining (red triangles) silver fir trees. Solid lines represent predicted relationships by generalized least square models (see Table 4) and dashed lines are 95% confidence intervals

DISCUSSION

Linkages between growth patterns and isotopic signals

The onset of growth decline in DD trees several decades ago (Figures 1 and 2) is aligned with previous studies reporting canopy dieback in sub‐Mediterranean silver fir populations in the Spanish Pyrenees (Büntgen et al., 2014; Camarero et al., 2011, 2015; Gazol et al., 2018b, 2020; Hevia et al., 2019). The occurrence of successive severe droughts and increasingly warmer temperatures during the 1980s–1990s may have impaired growth and resilience and led to very low growth rates that preceded dieback (Cailleret et al., 2017; Camarero et al., 2015, 2018; DeSoto et al., 2020). Indeed, the extreme 1985–1986 drought has been identified as a tipping point for the growth decline of several silver fir populations in the southwestern Spanish Pyrenees (Camarero et al., 2011), being aggravated by successive droughts in 2005 and 2012 (Camarero et al., 2018). This impact is more pronounced at low‐elevation, warm sites where drought is more intense than at high‐elevation, mesic sites (Camarero & Gazol, 2021; Gazol et al., 2015). Furthermore, management history and its interaction with climate may also have played an important role in determining dieback of disturbed silver fir forests in the Pyrenees (Sangüesa‐Barreda et al., 2015). Intense logging during the 1950s shaped the current populations as severe and selective thinning removed most of the large trees and promoted the persistence of the smallest and slowest growing trees, which seem not able to cope with water deficit and competition (Camarero et al., 2011). By contrast, the absence of forest management in CO due to the protection of the national park would have prevented a logging‐induced negative selection of trees.

Our results are also in line with the well‐established growth sensitivity of silver fir to water deficit and the associated low soil moisture and elevated evaporative demand during summer (Camarero et al., 2015; Gazol et al., 2019; Pasho et al., 2011; Sánchez‐Salguero et al., 2017; Vicente‐Serrano et al., 2015), as we found the strongest correlation of radial growth with August SPEI at 3‐month scale. At the warmer and drier, low‐elevation, managed sites, DD trees showed stronger climate–growth relationships and lower resilience against drought than ND trees, in agreement with previous results of forest stands undergoing dieback (Colangelo et al., 2017; González de Andrés & Camarero, 2020; Linares & Camarero, 2012b; Sangüesa‐Barreda et al., 2015). Therefore, different susceptibility to decline among individuals appear to be driven by intraspecific differences in growth sensitivity to drought stress among neighboring individuals. This is partly supported by the lower iWUE displayed by DD trees in PE (Figure 3b). Less efficient use of water by more defoliated trees has been previously found in silver fir, associated with a loss of responsiveness to rising Ca (Linares & Camarero, 2012a; Pellizzari et al., 2016) and other tree species (González de Andrés & Camarero, 2020; Hentschel et al., 2014; Salmon et al., 2015), although the opposite pattern is also found in the literature (Camarero, Colangelo, Gazol, et al., 2021; González de Andrés et al., 2021; Pellizzari et al., 2016). Interestingly, Puchi et al. (2021) have recently found different water‐use strategies between vigor classes of Araucaria araucana depending on precipitation and site conditions. This pattern suggests loss of hydraulic performance of DD trees as indicated by the smaller lumen areas of tracheids of this vigor class previously reported in silver fir (Pellizzari et al., 2016).

Whether the lower iWUE of DD trees is driven by reductions of photosynthetic rates and/or less tight stomatal control could be resolved using leaf δ18O information (Gessler et al., 2018). The absent or even negative iWUE–δ18O associations (Figure 3) revealed an uncoupling between leaf‐level processes and oxygen isotope signature (Barbour et al., 2000), suggesting that carbon source‐sink dynamics are largely driven by photosynthetic capacity rather than by stomatal responses to soil moisture and evaporative demand (Barbour et al., 2002; Billings et al., 2016). Actually, there is growing evidence supporting the predominant effect of the variability in source‐water δ18O composition on leaf and tree‐ring δ18O isotopic signals in many ecosystems (Ding et al., 2021; González de Andrés et al., 2021; Sarris et al., 2013; Shestakova et al., 2014; Treydte et al., 2014; Voltas et al., 2015). The steep vertical gradient in evaporative isotopic enrichment of soil water that develops in drying soils during prolonged rainless periods causes strong isotopic enrichment near the surface that exponentially decreases with depth (Allison et al., 1983). So, the approximate depth of soil/bedrock water uptake by trees can often be roughly inferred from the δ18O signals imprinted on foliar cellulose by the isotopic composition of their water sources (not measured in the present study). The lower foliar δ18O signature in DD trees at warmer, drier, formerly managed sites hinted at their utilization of deeper water sources than their co‐occurring healthy neighbors (Figure 3), since deep water pools usually show more depleted δ18O values than topsoil water, due to their lower exposure to evaporative isotopic enrichment in deep layers (Sarris et al., 2013). Similar findings were reported for a Nothofagus dombeyi population undergoing dieback (González de Andrés et al., 2021). However, other studies on oak species reported the preferential use of shallow water resources by DD trees, which would be less capable to reach deeper groundwater than less vulnerable ND trees (Ripullone et al., 2020).

Nutritional impairment of declining trees

An outstanding result of this study is the consistently lower foliar nutrient concentrations of DD trees compared to ND trees across study sites (Table 2; Figure 4 and Figure S3) as previously Hevia et al. (2019) showed in tree‐ring nutrient series. Two non‐mutually exclusive hypotheses may explain the observed patterns in tree nutritional status. On one hand, at the low‐elevation, warm sites, the greater utilization of water reserves from deeper soil/bedrock layers by drought‐stressed DD trees may imply disadvantages in terms of nutrient availability and uptake, as most soil nutrients are located in the topsoil horizons (Jobbágy & Jackson, 2001). This interpretation is strongly supported by the frequent positive correlations between foliar δ18O values and the leaf concentration of P, K and essential micronutrients such as Cu and Ni at these sites (Table S3; Figure 4 and Figure S4), suggesting enhanced nutrient uptake in trees using a greater proportion of shallow, isotopically enriched water from fertile topsoil horizons. This interpretation is also consistent with the conceptual framework proposed by Querejeta et al. (2021), who found that heavier exploitation of subsoil/bedrock water induced by top‐soil desiccation leads to the deterioration of plant mineral nutrition due to the vertical decoupling between water and nutrient availability and uptake in soil/bedrock profiles. Differences in water uptake depth and nutrient status between vigor classes may originally arise from microtopographic features or root habits (e.g., root dept and density), e.g., by DD trees growing on drier, less fertile or more exposed microsites or in sites with shallower soils with lower water‐holding and storage capacity. Individual genetic predisposition to drought‐induced dieback could also drive divergences between ND and DD trees, and interact with microenvironmental conditions, given that drought resilience of silver fir has been associated to genes linked to photosynthesis and drought stress (Heer et al., 2018); although, it does not seem a plausible explanation in this study since we compared trees within the same populations. Eventually, severe defoliation of DD trees should lead to sharp declines of overall canopy conductance and transpiration, thereby reducing the transpiration‐driven mass flow of soil nutrients to roots and subsequently decreasing nutrient uptake by trees (Salazar‐Tortosa et al., 2018). Poor nutrient status and selection of conservative traits in DD trees may also be a consequence of management legacies at the low‐elevation, managed sites (Camarero et al., 2011). We found a consistently negative correlation between leaf nutritional status and the competitive pressure at the declining sites (Fig. S7). Tree‐to‐tree competition has been regarded as a modulating factor of the adaptation capacity of tree growth to climatic stress in Abies species (Lebourgeois et al., 2014; Linares et al., 2010). Our study hints at increasing competition for water due to drought‐induced topsoil dessication that in turn impairs nutrient uptake and thus tree physiological functioning under drought.

Whatever the initial triggering mechanism, the deteriorated mineral nutrition of DD trees likely exacerbated and amplified the deleterious effects of droughts on physiological functioning and growth (Gessler et al., 2017). Foliar nutrient concentrations at the non‐declining site (SO) are considerably higher than those found at the declining sites, particularly for P, K; Cu and Ni (Fig. S4), so reinforcing the idea of a tight linkage between tree water balance and nutrition under water deficit. Moreover, the positive correlations between iWUE and foliar concentration of P and K (Figure 4) suggest nutrient limitation of photosynthesis in silver fir considering the relevant role of these essential elements in carbon assimilation and primary productivity (Güsewell, 2004; Sardans et al., 2013). Therefore, diminished photosynthetic capacity of highly defoliated trees showing poorly functional leaves with low nutrient contents could foster the growth decline patterns found in DD trees. Moreover, P and K also play key roles in mechanisms of avoidance of water stress, stomatal function, regulation of hydraulic conductivity and internal plant osmoregulation (Sardans & Peñuelas, 2007; Trifilò et al., 2011). Thus, the reduced concentration of P and K in leaves of DD trees could also increase the likelihood of hydraulic failure under severe droughts.

The significant concentration of Ca and certain micronutrients such as Mn, Sr and Si in leaves of DD trees (Table 2; Figure S4) may have several explanations. First, Ca and K are highly antagonistic and interfere with each other during absorption by the roots (Chapin, 1980; Diem & Godbold, 1993). Due to the limestone lithology of the study sites, the soil solution may be highly Ca saturated, thereby enhancing K deficit in trees. Deeper soil water uptake by drought‐stressed DD trees at low‐elevation, managed sites induced by topsoil desiccation could intensify this pattern since K and Ca availability decreases and increases with soil depth, respectively. The negative correlations between foliar δ18O and Ca concentrations (Table S3; Figure 4 and Figure S4) further suggest greater Ca uptake and accumulation with greater water uptake from deeper subsoil/bedrock layers in limestone lithologies. Second, the roots of phosphorus‐deficient DD trees may exudate organic acids and carboxylates as a mechanism for improving mineralization and mobilization of inorganic and organic P in their rhizospheres (Lambers et al., 2013). This strategy also enhances the availability and uptake of other micronutrients such as Mn, thereby increasing their concentration in plant tissues (Lambers et al., 2015), coupled with an accumulation of Ca (Mizuno et al., 2013). However, the elevated carbon cost of releasing metabolically expensive organic acids and carboxylates to rhizosphere soil for P scavenging may have contributed to carbon depletion and loss of vitality in DD trees. We also found consistent higher values of the Mn:Al ratio in DD trees across scales (Figure S4). Nutrient imbalances related to Mn have been considered as predisposing or contributing factors in forest dieback (González de Andrés et al., 2021; Hevia et al., 2019; Houle et al., 2007; Kogelmann & Sharpe, 2006). Accumulation of Mn in leaves has been associated to impaired photosynthesis and growth due to metabolic interferences with other nutrients (St. Clair et al., 2005). This study provides evidence in support the usefulness of Mn and related ratios as early warning signals of forest dieback in response to drought.

Coordination of leaf and stem traits in relation to drought responses

Growth resistance, recovery and resilience to drought were greater in trees with higher foliar concentrations of essential nutrients (Figure 6). The strong effects of foliar nutrient concentrations on drought performance (Table 4) support the key role played by mineral nutrition in tree functioning and growth before, during and after drought (Kreuzwieser & Gessler, 2010). Reduced foliar concentrations of essential nutrients such as P and K can intensify drought effects on photosynthesis and stomatal regulation, and thus have potential impacts and negative feedbacks on hydraulic failure and carbon starvation (McDowell et al., 2011). After drought, tree nutritional status also modulates re‐establishment of physiological functions and the repair of tissues damaged during drought (Gessler et al., 2017).

Differences between vigor classes also extended to leaf and wood traits, which showed tight coordination with tree mineral nutrition (Fig. S7). However, wood density and leaf morphological traits appeared to be poor predictors of within species variability in drought sensitivity, which concurs with previous studies on silver fir (George et al., 2015; Mihai et al., 2021; Serra‐Maluquer et al., 2021b). In general, our results suggest that intraspecific variability towards trait constellations denoting more conservative resource use strategies prompts growth decline and nutritional impairment of silver fir trees. As predicted by the leaf economic spectrum (Wright et al., 2004), slow‐growing DD trees showed smaller leaves with lower nutrient concentrations (Figure 5). These results are consistent with patterns reported within and across species (e.g., Greenwood et al., 2017; Valladares & Sánchez‐Gómez, 2006). Conversely, studies analyzing the relationship between wood density and drought response have provided unclear, scale‐dependent results. Greenwood et al. (2017) found that species with denser wood are less susceptible to drought in a global study. In contrast, positive correlations between wood density and growth rate and drought‐induced mortality have been also found (Hoffmann et al., 2011). Furthermore, translating the interspecific patterns to the intraspecific level remains challenging (Fajardo, 2016).

The negative relationship between foliar P, K and Cu leaf concentrations and soil assimilable P (Fig. S7) agrees with the work of Sardans and Peñuelas (2007), who found increases in soil P and K concentrations as a result of reduced tree uptake owing to drought. The complexity of the interactions between tree functioning, stand structure and soil would explain the absence of a direct effect of soil properties on short‐term growth responses to drought (Serra‐Maluquer et al., 2021a). Only pH exerted a consistent, positive impact on drought resilience (Table 4). The negative impact of soil acidity on tree growth and vitality have been previously reported on silver fir (Pinto et al., 2008; Rozas & Sampedro, 2013; Tallieu, 2020), and it could be explained by increased leaching losses of base cations, impairment of fine roots development and/or hampered mineralization under low pH conditions (Högberg et al., 2006; Viet et al., 2013).

Finally, trees at the high‐elevation, unmanaged site (CO) did not follow the same patterns in climate–growth relationships and isotopic composition, as no significant differences in leaf or wood traits were found between vigor classes (Table 2; Figure 3a). These results, together with the higher growth of DD trees during most of the 20th century (Figure 1), suggest that canopy dieback at site CO may be induced by a combination of different factors. The detection of the fungal pathogen Heterobasidion annosum in the study area (personal observation) could have contributed to the dieback at this site. Although typically associated to logging activities (Oliva & Colinas, 2010), H. annosum can also disperse from fallen, uprooted and dead trees in unmanaged forests with abundant gaps such as CO (Sangüesa‐Barreda et al., 2015). Severe dry spells usually reduce carbon reserves of trees so reducing investment in defenses, what facilitates the access of pathogens (McDowell, 2011; Oliva et al., 2014). H. annosum may have impaired water and carbon balance at CO by forcing trees to allocate photoassimilates to defense responses at the expense of radial growth, which results in lower hydraulic conductivity and sapwood storage (Oliva et al., 2012). The bigger size and higher growth rates of DD trees (Table 1; Figure 1) is consistent with the preferential attack of this primary pathogen on previously dominant or fast‐growing trees (Cherubini et al., 2002).

CONCLUSIONS

Impaired nutritional status of declining trees of Pyrenean silver forests was supported by their lower concentrations of essential nutrients such as P, K and Cu in leaves, a time‐integrated proxy of tree nutrition. This may become crucial due to feedbacks between nutrients and tree physiological functioning involved in coping with drought. Declining trees also showed stronger climate sensitivity and lower growth resilience against drought than non‐declining trees. Site‐dependent patterns in isotopic composition prevent us to outline consistent conclusions on water and carbon balances in silver fir dieback process at regional scale. However, the absence or negative relationships between carbon and oxygen signatures point to the uncoupling of the latter from leaf‐level processes. Frequent positive correlations between δ18O and essential nutrients in two of the three declining forests suggest the exploitation of subsoil/bedrock water induced by topsoil dessication in declining trees considering the vertical decoupling between water and nutrient availability in soil/bedrock layers.

Our results suggest that forthcoming studies may benefit from incorporating nutritional aspects for disentangling mechanisms underlying forest dieback and mortality. Indeed, we have found that tree‐level differences in nutritional status provide similar, on even greater, insights into the drought performance patterns than intraspecific trait variation and soil physical‐chemical properties.

CONFLICTS OF INTEREST

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript or in the decision to publish the results.

AUTHOR CONTRIBUTIONS

AG, RS‐S, JCL, JIQ and JJC planned and designed the research; EGA, AG, MC, RS‐S, JCL and JJC collected the data; EGA, JIQ and JMI analyzed the data; EGA, AG, JIQ and JJC led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Supplementary Material

Click here for additional data file.