Introduction 

Human activities have caused approximately 1.0 °C of global warming in 2020 above values recorded in pre-industrial times, with important changes in the intensity and the frequency of some climate and weather extremes ([38]). Warming caused by anthropogenic greenhouse gases emissions (GHGs) will persist for centuries to millennia ([46], [68]) and will continue to cause further long-term changes either in the climate system or in the biosphere. In particular, atmospheric CO₂ concentration has been constantly rising in the last centuries because of human activities, reaching in the 2000s the highest levels over the last 160,000 years ([37]). In forest ecosystems, climate change is altering productivity, species distribution, tree physiology and pest disease severity and frequency, with widespread tree mortality phenomena ([2], [3], [41], [65], [64]). However, increase in CO₂ atmospheric concentration can stimulate tree growth, can modify plant physiology and finally affect forest dynamics ([44], [16], [39], [57]). Indeed, elevated CO₂ concentrations can stimulate tree growth due to enhanced photosynthesis ([1]). Moreover, higher atmospheric CO₂ levels (C_a) usually lead to stomata closure, thus decreasing transpiration rates and increasing assimilation rates ([66], [80]). This leads to an improvement of water use efficiency (WUE) defined as the ratio between photosynthesis (A) and transpiration. Multiplying WUE by vapour pressure deficit (VPD) yields to the intrinsic water use efficiency (iWUE = A/g_s, where g_s is the stomatal conductance), which is not sensitive to increased transpiration driven by abiotic changes ([19], [61]). Higher iWUE values, coupled with faster tree growth measured by tree-ring width and converted into basal area increment (BAI), have been previously reported for several forests and linked to the increase in atmospheric CO₂ ([22], [42], [63]). However, this relationship has been not always confirmed ([45], [76], [73]) as an increase in iWUE alone may not directly translate into higher BAI as other factors (e.g., high temperature, recurrent drought, nutrient limitation) may negatively influence tree growth ([72]). More recently, Guerrieri et al. ([33]) reported an increase in iWUE in the last thirty years in 8 North America mature temperate forests and linked it to an overall increase in photosynthesis rather than to a reduction in stomatal conductance at most sites.

Tree rings offer insight into lifetime growth patterns, allowing climate impacts on trees to be evaluated ([27]). On the other hand, stable carbon (C) isotopes represent a useful tool to better understand changes in tree growth and productivity related to climate conditions through time. The variability of the C isotopic signature (δ¹³C) in wood gives information about the plant’s sensitivity to different climatic and environmental conditions ([48]). C isotope discrimination (Δ¹³C) can be related to CO₂ flux through stomata and to changes in C_a to intercellular CO₂ concentration (C_i) gradient and, by consequence, to iWUE ([25], [21]). However, the extent to which rising C_a has affected long-term iWUE, and whether climate could explain deviations from expected C_a-induced growth enhancement, is still poorly understood ([72]). On the other hand, environmental variation of δ¹⁸O in tree rings reflects the wide range of variation in meteoric water, soil water and transpiration ([29]). Barbour & Farquhar ([5]) found a strong relationship between leaf water evaporative enrichment (Δ¹⁸O_L) and g_s. Moreover, the same authors proposed an equation describing the relationship between Δ¹⁸O_L and ¹⁸O enrichment of cellulose (Δ¹⁸O_c), taking into account the exchange of carbonyl oxygen with unenriched xylem water. Based on these findings, the derived Δ¹⁸O_L was used to constrain the contribution of stomatal conductance in driving observed changes in iWUE ([33]).

Old-growth stands have developed for a long period of time without relevant human impact and/or stand replacing or major natural disturbances, and have three main structural characteristics: old and large trees, abundant coarse woody debris in different decay stages and a multilayered vertical structure ([69]). Old-growth forests have an important role in maintaining animal and plant diversity, are important C reservoirs and are privileged sites to study long-term plant physiological responses, long-term forest dynamics and climate change impact on forest ecosystems ([49]). Old-growth forests have been often considered to be insignificant as carbon sinks because of an equilibrium between assimilation and respiration ([54]). This hypothesis is based on studies showing a decline in net primary productivity at stand level with stand age ([81], [31], [59]) and in photosynthesis at tree level ([36]) and the general idea that ecosystem respiration increases with stand age ([54]). Some recent works have suggested that such a decline is not always observed at stand level ([43], [82]) and that old-growth forests are still accumulating C, as stand structure rather than age determines the overall forest capacity to absorb C from the atmosphere ([84], [43], [56]). However, as far as single trees are concerned, it has been proved that taller trees differ physiologically from shorter, younger trees ([58]): several studies have highlighted how old-living trees undergo age-related declines with hydraulic limitations (i.e., reduction in stomatal conductance and in leaf-specific hydraulic conductance) and reduction in photosynthesis ([81], [60], [58]). In this context, we focused our attention on old-living silver fir (Abies alba Mill.) trees in three mixed old-growth forests in the Balkans with the aim to better understand the relationship between growth dynamics, atmospheric CO₂ concentration increase and climate conditions through time. By measuring cellulose δ¹³C in tree-rings, the present study hypothesised an increase in the iWUE of old-living silver trees during the last two centuries. Moreover, using cellulose δ¹⁸O data and BAI we tested the hypothesis that the increase in iWUE was mostly related to enhanced photosynthesis rather than reduced stomatal conductance.

  Materials and methods 

Study area

Three study sites were selected in three different old-growth forests, along a 100 km transect from North-west to South-east, in the Balkans, in South-eastern Europe (Fig. S1 in Supplementary mayterial). All of the three sites are in the montane belt and are mixed with silver fir (Abies alba Mill.), Norway spruce (Picea abies Karst), and European beech (Fagus sylvatica L.).

The first site (LOM) is located in the Lom forest reserve (44° 27′ N, 16° 28′ E; 1250-1500 m a.s.l.) in the Dinaric Alps in Bosnia and Herzegovina. Climate is continental with maritime airstream influences. Mean annual precipitation is about 1600 mm and mean annual temperature is 7.6 °C ([10]).

The second site (PER) is located in the Perućica forest reserve (43° 19′ N, 18° 40′ E; 600-2800 m a.s.l.), inside the Sutjeska National Park in the southern Dinaric Mountains, Bosnia and Herzegovina. Climate is a mix of Mediterranean and continental, with mean annual precipitation of 1400 mm and mean annual temperature of 11.3 °C ([51]).

The third site (BIO) is located in the National Park of Biogradska Gora (42° 53′ N, 19° 36′ E; 830-2100 m a.s.l.) in the north-eastern part of Bjelasica mountain range, Montenegro ([50]). The mean annual precipitation at Biogradsko lake (1093 m a.s.l.) is 1962 mm, with a mean annual temperature of 5 °C.

Tree sampling and dendrochronological analysis

To characterize the forest’s structure, a regular 120-m grid was superimposed to the 1:10.000 raster map of each forest and 30 to 40 sampling points were randomly selected at each site (in 2014, 2016 and 2017 in PER, LOM and BIO sites, respectively). At each point, a 615.5 m² circular plot was set up for the measurements of the living trees (dbh ≥ 7.5 cm), tree species were assessed and tree mean diameter at 1.30 m (dbh) was measured ([50]). Moreover, one permanent plot (1 ha) was set up for long-term monitoring purposes at each forest.

For the purposes of this study, we decided to consider only silver fir as it was the dominant species in term of basal area at all sites (Tab. 1), because of its longevity and because it mostly occupied the upper canopy layer (PER is reported in Fig. S2 as example). Thus, one tree core was taken at the bottom of the stem from 7-43 individuals of silver fir within the top canopy layer (i.e., the tallest trees in the stand) depending on the site (Tab. 1, Fig. S2 in Supplementary material), using an increment borer. Cores were stored in cardboard and, once in laboratory, samples were glued on woody supports and sanded with paper with 60, 120, 240, 320 and 600 grit to better highlight tree rings.

Tree-ring widths were measured for each core using a LINTABTM (Frank Rinn, Heidelberg, Germany) and then were visually cross-dated using TSAPWin® software ver. 4.81 (Rinntech-Metriwerk GmbH & Co. KG, Heidelberg, Germany). COFECHA software ([32], [34]) was used for the statistical cross-date and to check chronologies quality ([13]).

Isotopes analysis

From all the samples at each site, the six longest chronologies were selected. Then, tree rings for each core were grouped using 10-year intervals, cut and separately milled using an ultra-centrifugal mill with a sieve of 0.5 mm (ZM 100®, Retsch Tecnology, Haan, Germany). In order to have enough replications at each plot, we considered only five consecutive time-intervals along each time-series for the further analysis (1800-1810; 1850-1860; 1900-1910; 1950-1960; 2000-2010). Wood powder was treated with a two-step digestion to extract cellulose according to Green ([30]). Briefly, samples were sealed in Teflon bags and at first treated two times with NaOH 5% solution for 2 hours at 60 °C. Then, samples were treated with NaClO₂ 7% solution, adjusting the pH to 4-5 with acetic acid, for 10 hours at 60 °C. Depending on sample weight and its cellulose content, this phase was repeated for three or four times until samples totally bleached, thus indicating that all fibres but cellulose were completely digested. Dried sub-samples were finally weighted in tin capsules for δ¹³C for analysis using a CHNS elemental analyser (Vario Microcube®, Elementar, Langenselbold, Germany) coupled with a stable isotope ratio mass spectrometer (IRMS Isoprime® 100, Elementar). We also determined δ¹⁸O of a sub-sample of cellulose weighted in silver capsules for each considered time period using a TC/EA elemental analyser (Flash 2000®, ThermoScientific, Waltham, MA, USA) connected to an isotope ratio mass spectrometer (Delta V Advantage®, Thermo Scientific) via a continuous flow interface (ConFlo IV®, Thermo Scientific). Carbon and oxygen isotope ratios were expressed in per mil (‰) relative to the Vienna Pee Dee Belemnite (VPDB) and Vienna Standard Mean Ocean Water (VSMOW) international standards, respectively.

Data analysis

The C isotopic discrimination (Δ¹³C) for each selected period and tree was calculated according to Farquhar et al. ([20] - eqn. 1)

where δ¹³C_a is the isotopic signature of atmosphere and δ¹³C_c is the isotopic signature of the cellulose. As our δ¹³C_c values represent an average value for 10-year intervals, δ¹³C_a were calculated as the mean value for the same temporal interval of the relative δ¹³C_c measurements using values published in McCarroll & Loader ([48]) and those from Mauna Loa Observatory since 2003 (⇒ https:/­/­www.­esrl.­noaa.­gov/­gmd/­).

We used the carbon isotope discrimination model for C3 plants (eqn. 2) proposed by Farquhar et al. ([20]) and Francey & Farquhar ([25]) to compute intercellular CO₂ concentration (C_i, ppm) based on Δ¹³C and atmospheric CO₂ concentration (C_a, ppm; data from Mauna Loa Observatory - ⇒ https:/­/­www.­esrl.­noaa.­gov/­gmd/­):

where a is isotope fractionation during CO₂ diffusion through stomata (a = 4.4‰), b is isotope fractionation during carboxylation processes (b = 27‰). According to the Fick’s law, net photosynthesis (A) can be calculated as follows (eqn. 3):

where g_CO2 is the leaf conductance to CO₂. Knowing that leaf conductance to water vapour (g_H2O) is equal to 1.6 · g_CO2, we calculated iWUE (μmolCO₂ molH₂O^-1), expressed as the ratio between photosynthesis and stomatal conductance to H₂O (A/g_H2O) combining all the previous equations (eqn. 4):

We also calculated the ¹⁸O enrichment of cellulose (Δ¹⁸O_c) according to eqn. 5:

where δ¹⁸O_P is the annual δ¹⁸O of the precipitation estimated following Barbour et al. ([4] - eqn. 6):

where T_a, P_a and E are the mean annual temperature (°C), the annual precipitation (m) and the elevation of the weather station (m), respectively. We assumed that there were no differences among trees in the depth water was taken up and, thus, that the soil water δ¹⁸O reflected precipitation δ¹⁸O. δ¹⁸O_P was computed as the average between the values of the weather station located in Sarajevo (43.8678° N, 18.4228° E - 630 m a.s.l.; data availability: 1901-2019) and the one located Zagreb (45.8167° N, 15.9781° E - 157 m. a.s.l.; data availability: 1881-2019) as these two stations showed similar rain and temperature averages across all the data period. Climate data were downloaded from Global Historical Climatology Network website (GHCND). Long term δ¹⁸O_P trends at the two stations are reported in Fig. S3 (Supplementary material).

We finally estimated the leaf water ¹⁸O evaporative enrichment above the source water (Δ¹⁸O_L) according to Barbour & Farquhar ([5]) and Gessler et al. ([29]):

where p_ex is the proportion of exchangeable oxygen in cellulose formed from sucrose, p_x is the proportion of unenriched (source) water in the developing cambium cell and ε_WC is the isotope fractionation occurring during cellulose synthesis in the stem (27‰). Generally, a fixed value of 0.40 is considered for p_xp_ex.

Yearly basal area increment (BAI_i, cm² yr^-1) was calculated using tree-ring width as follows (eqn. 8):

where R_i+1 is the radius at the end of the year i and R_i is the radius at the beginning of the year i, Δr_i is the tree-ring width at year i.

Statistical analysis

All statistical analyses were performed using R software ver. 3.6.0 (R Foundation for Statistical Computing 2019 - ⇒ https:/­/­www.­r-project.­org/­). We compared iWUE, c_i, Δ¹⁸O_L and BAI in the selected decades at each site by using one-way ANOVA and Tukey’s post-hoc test when significant (p<0.05) differences were detected. All data were eventually log-transformed before analysis to meet the requirements for parametric statistical tests using “powerTransform” and “bcPower” functions of the “car” package ([24]). Linear regressions were used to explain changes in iWUE with BAI; before regression, data normality was checked using the Shapiro-Wilk test. All data throughout the text and in tables and figures are reported as mean ± standard error.

  Results 

The mean stand characteristics at each experimental site are reported in Tab. 1. Total stem density ranged from 412 ± 22 to 489 ± 19 trees ha^-1 at BIO and LOM, respectively, but total basal area followed an opposite trend, being the highest at LOM and the lowest at BIO (60.1 ± 4.4 and 47.1 ± 1.8 m² ha^-1, respectively). European beech was the dominant species in terms of stem density (n ha^-1) at all the three sites, but silver fir was always dominant in terms of basal area (m² ha^-1). The oldest measured tree ring series were 306, 267 and 273-year-old in LOM, PER and BIO, respectively, though it was not always possible to reach the centre of the stem.

C_i significantly increased through time at all sites (LOM: p<0.001; PER: p<0.001; BIO: p=0.01 - Fig. 1; Tab. S1 in Supplementary material) and was significantly related to C_a (LOM: R²=0.57, p<0.001; PER: R²=0.93, p<0.001; BIO: R²=0.64, p<0.001). We did not find any significant difference in the ratio C_i/C_a at BIO (p=0.06), but a significant decrease in C_i/C_a was measured for the periods 1950-1960 and 2000-2010 at LOM (p<0.001) and for the period 2000-2010 at PER (p=0.04 - Fig. S4, Tab. S5 in Supplementary material). Trees also showed a significant increase in their iWUE, as well as in their BAI, at all sites (Fig. 2; Tab. S2, Tab. S3), reaching significant higher values especially at the end of the last century.

Fig. 1 - Intercellular CO₂ concentration (C_i, ppm) by selected decades at the three old-growth forests (PER: Perucica; LOM: Lom; BIO: Biogradska Gora). Means ± standard error are displayed. Different letters indicate significant (p<0.05) differences among periods after post-hoc Tukey’s test.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 2 - Intrinsic water-use efficiency (iWUE, μmolCO₂ molH₂O^-1), basal area increment (BAI, cm² yr^-1) and leaf water evaporative enrichment (Δ¹⁸O_L) by selected decades at the three old-growth forests (PER: Perucica; LOM: Lom; BIO_ Biogradska Gora). Means ± standard error are displayed. Different letters indicate significant (p<0.05) differences among periods after post-hoc Tukey’s test.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

As far as δ¹⁸O_C is concerned, we did not detect any significant difference among the selected periods at all the three sites (data not shown; p>0.05). Moreover, we observed no changes in Δ¹⁸O_L at all the three sites for the periods when δ¹⁸O_P was available (1900-1910; 1950-1960; 2000-2010; p>0.05 - Fig. 2, Tab. S4 in Supplementary material).

Fir trees showed a significant increase in BAI with time (Fig. 2, Fig. 3; Tab. S2 in Supplementary material) at all the three experimental sites, especially after the beginning of the 19^th century. BIO had generally the highest growth rate, though it showed a slight decline in the last two decades, while fir at PER had generally the lowest growth, but did not show any decrease in growth with time.

Fig. 3 - Mean annual basal area increment (BAI, cm² tree^-1 yr^-1 - black lines) of the sampled silver fir trees (PER: Perucica; LOM: Lom; BIO: Biogradska Gora). The grey lines represent the standard error.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Using simple linear regressions, we found that iWUE significantly increase with BAI at all sites (LOM: R²=0.68, p<0.0001; PER: R²=0.42, p=0.0003; BIO: R²=0.45, p=0.0002 - Fig. 4; Tab. S6 in Supplementary material).

Fig. 4 - Basal area increment (BAI, cm² tree^-1 yr^-1) with increasing intrinsic water-use efficiency (iWUE, μmolCO₂ molH₂O^-1). Each point represents a selected decade for each sampled tree. Black line represents the regression line and grey lines represent 95^th confidence intervals.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  Discussion 

In our study, iWUE showed a significant increase in the last two centuries at LOM, PER and BIO, (Fig. 2, Tab. S3). These trends are similar to those measured in a wide range of ecosystems ([55], [67], [26]) either in temperate ([76], [62], [33]) or tropical/sub-tropical forests ([79], [78]), as well as at single experimental sites ([6]) or across the whole Europe ([63]), but none of these studies considered old-growth forests.

Variations in both A and g_s influence the overall change in iWUE ([19], [61]). The increase in atmospheric CO₂ concentration positively affects A by increasing the carboxylation rate of Rubisco ([17], [1], [35]) and positively influences tree growth rate ([40], [71]). In our case, C_i increased proportionally to C_a, C_i/C_a stayed constant over the last two hundred years (BIO) or slightly decreased at the end of last century (PER, BIO). The measured increase in C_i is consistent with what reported by Frank et al. ([26]) for Europe and by Guerrieri et al. ([33]) for North America and might indicate that the increase in CO₂ concentration is sustaining photosynthesis and the growth of the sampled dominant trees (Fig. 3) - the so called “CO₂ fertilization effect on A” ([47], [33]). We were not able to quantify the absolute increase in BAI for the whole stand (m² ha^-1 yr^-1) and to disentangle the role of climate change from small scale stand disturbance dynamics (i.e., release from suppression), because of our sampling strategy. In fact, taller trees may differ physiologically from shorter, younger trees ([58]) and the selection of only dominant trees within the stands might over-estimate absolute growth rates (“slow-grower survivorship bias” and “big tree selection bias” - [15], [11], [52], [18]). On the other, silver fir has the capacity for vigorous growth response after releases related to natural disturbances of the top-canopy layer or to management ([23]) and this makes it difficult to disentangle the CO₂ fertilization effect over time. But, the trend in BAI reported in Fig. 3can be considered a proxy of an overall increase at stand level in the last decades. Looking at Fig. 3, the sampled dominant trees did not show any significant decrease in tree growth after well-known heatwaves, such as that in 2003: these results seem to further confirm the high resistance of fir to drought extremes observed in other regions ([83], [74]) and the possibility of the species to cope with the climate change ([75]). The measured increase in BAI and the resistance to drought could be particularly important for the Dinaric Region, where silver fir has been one of the most ecologically and economically important tree species, mainly managed using the selection system ([7]).

On the other hand, increasing C_a may induce stomata closure and, thus, plants may improve their iWUE by reducing water loss by transpiration ([21], [47]). Fir, in particular, has been reported to be able to efficiently control its transpiration ([53]). In our case, as we did not observe any change in Δ¹⁸O_L, we can hypothesize a constant g_s over the last two centuries similarly to what has been reported by Guerrieri et al. ([33]) for mesic sites in North America in the last thirty years. This and the positive correlation between iWUE and BAI at all our three old-growth forest sites (Fig. 4, Tab. S6 in Supplementary material) further support the hypothesis that the measured increase in iWUE over time is more related to a stimulation of A rather than a reduction in g_s. A widely held assumption is that, after an initial period of increasing growth, the mass growth rate of individual trees declines with increasing tree size ([77]). On the contrary, Stephenson et al. ([70]) have shown that mass growth rate increases continuously with tree size in hundreds of tree species. Our BAI data suggest that large, old trees do not act simply as senescent carbon reservoirs, but actively fix large amounts of C and have been growing even faster in the last decades. These data are in agreement with the substantial increase in silver fir growth rate measured across several European regions probably promoted by warmer climate as well as by the increase in nitrogen depositions and atmospheric CO₂ concentrations ([8], [9], [14], [28]). Thus, the overall C sink potential of old-growth forests might be higher than previously thought as old-living dominant trees are not undergoing any age-related growth decline. The lower juvenile tree growth is also supporting the hypothesis that such a C sink capacity will be maintained for a longer period of time as trees will show a slower ontogeny ([12]).

  Conclusions 

In the last two centuries, we observed a continuous increase in iWUE in silver fir at all the three old-growth forest sites in the Balkans. By estimating leaf water δ¹⁸O from cellulose isotopic analysis in tree rings and using BAI measurements, we were able to separate the relative contribution of assimilation rates (A) and stomatal conductance (g_s). In particular, we observed an increase in iWUE, BAI, but not in δ¹⁸O_L. Thus, our data support some recent findings challenging the hypothesis that iWUE in forests is primarily the result of a CO₂-induced reduction in stomatal conductance. Moreover, our data support the idea that silver fir might cope with the climate change.

  Aknowledgements 

We would like to thank: the team of University of Turin (Department of Agriculture, Forestry and Food Sciences), Italy for helping during sampling in Bosnia-Herzegovina and in Montenegro; the team of WSL (Swiss Federal Institute for Forest, Snow and Landscape Research) in Birmensdorf, Switzerland for helping during cellulose extraction procedures in the lab; Gaia Zanin for sample processing and data analysis of the tree cores taken in Perucica forest in 2014.

  References

  Supplementary Material

Fig. S1 - Location of the three experimental sites.
Fig. S2 - Distribution of the standing trees by diameter and height class within the 1-ha permanent plot in Perucica and relative distribution of the sampled fir trees by diameter and height class in the same plot.
Fig. S3 - δ¹⁸O of the precipitation estimated following Barbour et al. (2001 - eqn. 8) at the weather station in Sarajevo and in Zagreb.
Fig. S4 - Ratio between intercellular CO₂ concentration (C_i) and atmospheric CO₂ concentration (C_a) by selected decades at the three old-growth forests.
Tab. S1 - ANOVA results for intercellular CO₂ concentration C_i by selected decades.
Tab. S2 - ANOVA results for basal area increment (BAI) by selected decades.
Tab. S3 - ANOVA results for intrinsic water use efficiency (iWUE) by selected decades.
Tab. S4 - ANOVA results for leaf water evaporative enrichment (Δ¹⁸O_L) by selected decades.
Tab. S5 - ANOVA results for the ratio between intercellular CO₂ concentration (C_i) and atmospheric CO₂ concentration (C_a) by selected decades.
Tab. S6 - Regression results of Fig. 4.