Site and sampling

We selected twelve monoculture tree stands of the most common tree species in Britain, Scots pine (Pinus sylvestris L.), Sitka spruce (Picea sitchensis Bong. Carr.), pedunculate oak (Quercus robur L.) and common beech (Fagus sylvatica L.). The majority of the stands were experimental sites within the Level II- ICP intensive forest monitoring network (, with the exception of Covert Wood, Shobdon and Goyt. The Goyt site was added as a high Ndep site as a contrast to the low Ndep Sitka spruce site in Scotland (Fig. 1, Table 1, Supplementary Table 1). For each species, forests were selected with similar soil type and age, but with contrasting Ndep, Sdep and climate, particularly rainfall and temperature, as described in Fig. 1, Table 1 and Supplementary Table 1. Stand information (mean tree height, mean diameter at the breast height and basal area) as measured for target years and for some of the forest stands are shown in Fig. S4.

At each ICP forest site, a plot of 0.25 ha was established in 1995 to carry out monitoring30 and a similar protocol was followed at the Goyt and Shobdon sites. Within each plot, 10 trees were selected for the collection of 3 wood cores per tree by using a 5 mm diameter increment borer, which were placed in paper straws for transport. Sampling was carried out between November 2010 and March 2011. Once in the laboratory, samples were dried at 70 °C for 48 h. Of the three wood cores sampled, one was kept for dendrochronology, with the other two used for stable isotope analyses.

Climate and atmospheric deposition data

Climate data (Temperature, T, Vapour Pressure Deficit, VPD, Precipitation, P) were obtained from automated weather stations at the sites and/or the nearest local meteorological stations (data were provided by the British Atmospheric Data Centre). Annual mean (Ta) and mean maximum (Tamax) values for temperature were calculated from monthly mean and maximum air temperature, T, respectively, and annual precipitation (Pa) was calculated as the sum of total monthly precipitations. Annual VPD was calculated from averaging monthly values obtained from mean monthly maximum temperature and minimum monthly relative humidity. For all the parameters, mean values were also calculated over the growing season, i.e., from May to September. We also considered the standardised precipitation-evapo-transpiration index, SPEI, relative to August, with 1 (SPEI8_1), 2 (SPEI8_2) and 3 (SPEI8_3) months time-scale and SPEI relative to December, with 1 and 12 months time-scale, the latter providing year-cumulated soil moisture conditions. SPEI values were obtained from the global database with 0.5 degrees spatial resolution available online at:

Yearly wet nitrogen (Ndep) and sulphur deposition (Sdep) were obtained from measured bulk precipitation and throughfall water volumes at the sites and measured elemental concentrations (NO3, NH4+ and SO2–4) as previously described30. For the spatial analyses, we considered mean of annual deposition (sNdep and sSdep), obtained as the sum of total (NH4-N + NO3-N for Ndep) wet and dry deposition. The latter were estimated as difference between throughfall and bulk precipitation N fluxes30. For Rogate only 1 year (2010) of monitoring was available. For Goyt site, atmospheric deposition data collected at Ladybower were considered, as the two sites are 30 km apart. For two sites (i.e., Shobdon and Covert Wood), which were not part of the regular ICP forest sites, the wet deposition obtained from the UK 5 × 5 km grid Ndep and Sdep dataset was used4. The estimate included wet and dry NHx-N (NH4, NH3), NOy-N (NO2, NO3, HNO3) and S (SOx = SO2 and SO4) deposition, modelled using FRAME with 2005 emissions data4. However, only the total wet deposition was included in the analyses here, as we previously reported a good agreement between modelled and measured wet Ndep50.

For the temporal analyses, only wet deposition (as calculated from NO3, NH4+ and SO2–4 concentrations in bulk precipitation) was considered (indicated as aNdep and aSdep), given the uncertainties associated with the quantification of the dry deposition. For instance, when differences between throughfall and bulk precipitation are < 0 it is assumed atmospheric deposition is retained by tree canopies, but this does not necessarily mean that there is no dry deposition. At Rogate only one-year data were available so we considered annual wet deposition data for Alice Holt, which is within 19 km distance. This was also the case for Goyt and Ladybower, which are 30 km apart. Shobdon and Covert Wood were not included in the analyses where annual deposition data were considered (see earlier in the text and Ref. Table S7).

Stable isotope analyses

Wood cores were subjected to removal of mobile N and extractives with a soxhlet apparatus as described in Guerrieri et al.23. After the chemical pre-treatment, wood cores were dated and cross-dated from 2010 back to 1980 and then separated with a scalpel as follow: single annual rings from 2010 to 1995 and then groups of 3 annual rings from 1994 back to 1980. We maintained the annual resolution from 1995 onward because this is the period when the UK-ICP forest network was established and atmospheric deposition was monitored.

To minimise the cost of the stable isotope analyses while including 12 sites and 4 different tree species, the wood materials were pooled from 10 trees (2 wood cores per tree) for each given ring or group of rings. However, for one year (i.e., 2007) and at two sites for each species, carbon and oxygen isotope ratios for each of the 10 trees was measured, so as to assess the variability among trees (Supplementary Table S5). For each core per tree species, each ring was cut with a razor blade under a microscope, milled and homogenized in a centrifugal mill, and then pooled by year. Moreover, for δ15Nw analyses, we only included the sites where long-term atmospheric deposition data were available (i.e., tree species at Shobdon and Covert Wood were excluded).

An amount of 0.4–0.6 mg of extracted wood sample from each given ring (or group of annual rings for the years 1994 back to 1980) was weighed in tin capsules and combusted in the elemental analyzer (NA2500, Carlo Erba) for the determination of δ13Cw by VG Prism III Isotope ratio mass spectrometer at the School of Geosciences (University of Edinburgh, UK). For δ15Nw, 25–30 mg of wood sample was weighed in tin capsules and combusted on a PDZ Europa ANCA-GSL elemental analyzer interfaced to a PDZ Europa 20–20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK). For δ18Ow, 0.8–1 mg of each sample was weighed in silver capsules and analyzed on a PyroCube (Elementar Analysensysteme GmbH, Hanau, Germany) interfaced to an Isoprime VisION (Isoprime Ltd., Stockport, UK, a unit of Elementar Analysensysteme GmbH, Hanau, Germany). Analyses were carried out at the Stable isotopes facilities of the School of GeoSciences (University of Edinburgh, UK) for δ13Cw, and at the Stable Isotope facility of the UC Davis, University of California (USA) for δ18Ow and δ15Nw. Stable isotope abundances are expressed as ratios of 13C/12C, 15 N/14 N and 18O/16O using δ-notation (in per-mil;‰) relative to international standards (VPDB for δ13Cw, atmospheric N2 for δ15Nw and VSMOW for δ18Ow). The standard deviation for internal standards was 0.1‰ for δ13Cw (PACS-2), 0.2, 0.3 and 0.4‰ for δ18Ow (IAEA 600, IAEA 601 and IAEA 602, respectively), and between 0.1 and 0.3‰ for δ15Nw (USGS-41 Glutamic Acid and peach leaves, respectively).

Calculation of iWUE and Δ18Ow

The iWUE and the ci/ca ratio were derived from measured δ13Cw values, and based on the well-established theory linking leaf ci/ca with carbon isotope discrimination, Δ13Cw,51 as shown in the equation below:

$$ Delta^{13} C_{w} = a + left( {b – a} right) frac{{c_{i} }}{{c_{a} }} = frac{{left( {delta^{13} C_{a} – delta^{13} C_{w} } right)}}{{left( { 1 + frac{{delta^{13} C_{w} }}{1000}} right)}} $$


where δ13Ca and δ13Cw are the carbon isotope compositions of ca and wood, a is the isotope fractionation during CO2 diffusion through stomata (4.4‰) and b is the isotope fractionation during fixation by Rubisco (27‰). Note that Eq. (1) is the “simplified version” of the Farquhar model describing carbon isotope discrimination in plant material, which does not include effects due to mesophyll conductance and photorespiration. We derived ci from the following equation:

$$ c_{i} = c_{a} frac{{(delta^{13} C_{a} – delta^{13} C_{w} ) – a}}{b – a} $$


ca values were obtained from Mauna Loa records27, and δ13Ca values were obtained from Mauna Loa records5 from 1990 to 2010, while from 1990 back to 1980 we used data published in Ref.52. iWUE (μmol CO2 mol−1 H2O) was then calculated using the following equation:

$$ iWUE = frac{A}{{g_{s} }} = frac{{c_{a} – c_{i} }}{1.6} = frac{{c_{a} }}{1.6} left( {frac{{b – Delta^{13} C_{w} }}{b – a}} right) $$


where 1.6 is the molar diffusivity ratio of CO2 to H2O (i.e., gCO2 = gH2O/1.6). Note that in the Eqs. (2) and (3) we used average of values measured over growing season months (May–September) for both ca and δ13Ca.

Tree-ring oxygen isotope discrimination, Δ18Ow, was calculated according to Eq. (4)53:

$$ Delta^{18 } O_{w} = frac{{delta^{18} O_{w} – delta^{18} O_{s} }}{{1 + left( {frac{{delta^{18} O_{s} }}{1000}} right)}} $$


where δ18Ow is the oxygen isotope composition we measured in each ring, while δ18Os is the oxygen isotope composition of the source water, i.e. the soil water, which we assumed to reflect the δ18O of precipitation (δ18OP). We estimated annual values of δ18OP at each site as described by Barbour et al.54, by considering the following equation:

$$ delta^{18} O_{P} = 0.52T_{a } – 0.006T_{a}^{2} + 2.42P_{a} – 1.43P_{a}^{2} – 0.046sqrt E – 13.0 $$


where Ta, Pa and E are the annual mean air temperature, precipitation (this latter expressed in m) and elevation (m asl), respectively. Mean values of estimated δ18OP from Eq. (5) were in line with estimates from the Online Isotopes in Precipitation Calculator ( and measured δ18OP values at Keyworth (Supplementary Table S6). The modelled δ18OP did not show a significant trend at most of the sites, with the exception of Goyt/Ladybower (slope = − 0.03 ± 0.007‰ per year, p < 0.001) and Covert wood (slope = 0.02 ± 0.009‰ per year, p = 0.05).

We assumed Δ18Ow to reflect the leaf water Δ18O, which is affected by transpiration. Notably, less enriched (in 18O) water from the soil and more enriched (in 18O) water at the leaf evaporative sites continuously mix, as a function of transpiration rates and the pathway of water movement through foliar tissues (i.e., Péclet effect)55 so that lower leaf Δ18O results from an increase in transpiration and gs56.

The physiological signal imprinted in the foliage may be dampened in tree rings, due to post-photosynthetic fractionation during translocation of sucrose and synthesis of cellulose in the tree stem57. This leads to an offset between foliar and tree-ring δ18O and also δ13C values. However, accounting for the offset when interpreting tree-ring isotopes is still challenging, as it is not clear whether the offset is species-specific, if it is maintained over the long-term and what are the mechanisms driving it57,58.

Statistical analyses

Linear regression analyses were initially used to explore whether (1) temporal trends existed between tree ring isotopes, iWUE and environmental data at each site (Fig. 2); (2) there was a relationship between iWUE and Δ18Ow and parameters describing stand development (diameter at the breast height and height, Supplementary Table S2 text S1); (3) changes in iWUE and δ15Nw were correlated with age and soil type (Supplementary text S2 and S3). Subsequently, we considered models jointly allowing explanation of spatial and temporal variation in isotopic data. To explain temporal trends in isotopic data (iWUE, Δ18Ow, ci/ca, δ15Nw) jointly across multiple sites, we considered both explanatory variables that varied yearly for each site and mean climatic data averaged over time for each site. Yearly time series and mean climatic data for each site (mean temperature, maximum temperature, precipitation and vapour pressure deficit VPD) were calculated for each year (from 1980 to 2010) and also separately for each growing season. We considered SPEI for the month of August with a time-scale of one, two and three months and then the SPEI for the month of December with a time-scale of one and twelve months. By definition, being centered around zero, SPEI defines yearly anomalies and cannot be used as a site index. We also included both the annual deposition data (i.e., aSdep and aNdep, to assess the effect of temporal changes) and the mean over the monitoring time (i.e., sSdep and sNdep) to evaluate both the contribution of within-site temporal changes and cross-sites differences on changes in iWUE, ci/ca, Δ18Ow and δ15Nw. Note that we only have data for the aSdep and aNdep at 10 of the 12 sites. To eliminate auto-correlation between mean site variables and yearly variables for each site, site variables were globally centered, whereas yearly data were group-centered59. To eliminate auto-correlation among individual climatic variables of the two groups, we conducted two principal components (PC) analyses, after centering and scaling the variables. The first PC analysis considered across sites long-term averages of environmental variables (PCA_s), the second within sites annual time series (PCA_a). Tables of percentage of variance explained and scree plots were examined to determine how many PC to retain. Linear mixed models (using library nlme in R60) were employed for each of the isotope data series to explain temporal and spatial patterns of variation as a function of climatic and atmospheric deposition conditions. We also ran the linear mixed model without the two Sitka spruce stands, in order to account for other source of variations (e.g., age-effect) for the isotope parameters and iWUE (particularly for the youngest stand at the Goyt site). To explore the significance of systematic differences among the 12 (or 10, when Sitka spruce was excluded) sites occupied by the two evergreen and the two deciduous species, a categorical variable with two levels (combining species into two functional groups) was introduced as a fixed factor. Since multiple species were present at some of the sites, an identification factor for each site × species combination was employed as random factor. The initial models included all PC previously identified as potential explanatory variables for both spatial and temporal variation and also forest stand-related (age and soil type) and anthropogenic (atmospheric CO2, Ndep and Sdep) factors (Ref. Supplementary Table S7 including all the models tested and reference to tables and supplementary text where results are reported). These models were then gradually simplified until the minimal significant model was achieved, i.e., excluding all PC and other factors that were not significant following simultaneous tests for general linear hypothesis testing (package multicomp,61). A correlation structure of order 1 was included in the model for each site × species combination to allow for the temporal dependency of measurements carried out in subsequent years. In the case of nested models, significance was tested using a chi-square test with one degree of freedom. Quality of fit was assessed using residual distribution plots, qqnorm plots of standardised residuals against quantiles of standard normals for both individual points and for the random effects, and auto-correlation function plots of normalized residuals as a function of measurement lags. Marginal (only fixed factors) and conditional (fixed + random factors) percent of explained variance (R2m and R2c, respectively) were calculated using package MuMIn62.

Finally, to examine the joint effects of climatic conditions on Δ18Ow, δ15Nw and iWUE, a mixed-effect confirmatory path model was employed using package piecewiseSEM63. For each isotope and iWUE, the final model from linear mixed effect model analyses (Tables 3, 4, Supplementary Table S7) was considered, with the only modification of excluding PFT as fixed factor and including also Δ18Ow as fixed factor in the model for iWUE and the correlated errors between PC_s1 and sNdep. PiecewiseSEM allows the fitting of hierarchical models with random effects on data with a multivariate structure, allowing for the identification of indirect effects and unresolved covariance among endogenous variables. Model goodness-of-fit was assessed using a chi-square test against Fisher C, based on Shipley’s test of direct separation, which tests for the overall significance of missing relationships among unconnected variables, while the significance of any given missing path was evaluated using individual p-values. A combination of non-significant Fisher-C and individual p-values tests implies that the hypothesised relationships among the variables are consistent with the data without missing any significant relationship63. All statistical analyses were carried out inside RStudio 1.0.143 with R 3.4.064.

Source link

Leave a Reply

Your email address will not be published. Required fields are marked *