Climate and atmospheric deposition effects on forest water-use efficiency and nitrogen availability across Britain.

Rising atmospheric CO2 (ca) has been shown to increase forest carbon uptake. Yet, whether the ca-fertilization effect on forests is modulated by changes in sulphur (Sdep) and nitrogen (Ndep) deposition and how Ndep affects ecosystem N availability remains unclear. We explored spatial and temporal (over 30-years) changes in tree-ring δ13C-derived intrinsic water-use efficiency (iWUE), δ18O and δ15N for four species in twelve forests across climate and atmospheric deposition gradients in Britain. The increase in iWUE was not uniform across sites and species-specific underlying physiological mechanisms reflected the interactions between climate and atmospheric drivers (oak and Scots pine), but also an age effect (Sitka spruce). Most species showed no significant trends for tree-ring δ15N, suggesting no changes in N availability. Increase in iWUE was mostly associated with increase in temperature and decrease in moisture conditions across the South-North gradient and over 30-years. However, when excluding Sitka spruce (to account for age or stand development effects), variations in iWUE were significantly associated with changes in ca and Sdep. Our data suggest that overall climate had the prevailing effect on changes in iWUE across the investigated sites. Whereas, detection of Ndep, Sdep and ca signals was partially confounded by structural changes during stand development.

Introduction

The ability of forests to mitigate climate change depends on how well they cope and adapt to the rapid forecasted changes in atmospheric conditions, including pollutant emissions from anthropogenic activities , namely CO2, sulphur (S) and reactive nitrogen (N). While sulphur emissions have been successfully regulated in Europe and North America[1], N compounds, particularly ammonia and nitrous oxide, continue at relatively high levels compared to the pre-industrial period[2-4]. These changes in N and S atmospheric concentrations have occurred together with increasing atmospheric CO2 concentration (ca), which recently exceeded 410 ppm[5], with a relative change of 46% compared to pre-industrial levels. Atmospheric N and S compounds are deposited from the atmosphere onto terrestrial ecosystems and together with the increasing ca can have an effect on some of the processes underpinning forest carbon, water and nutrient cycling.

Increases in ca affect leaf gas exchange by increasing photosynthesis (A) and reducing stomatal conductance (g)[6], thus raising intrinsic water use efficiency of trees (iWUE = A/g). Higher iWUE has been commonly reported for conifer compared to broadleaf species[7]. Moreover, the increase in iWUE derived from carbon isotope composition (δ13C) in tree rings[8] has generally been associated with an active response of trees to increasing ca, whereby intercellular CO2 (ci) increased in a proportional way to ca, resulting in a constant ci/ca ratio, due to proportional regulation of A and g[9]. However, changes in N (Ndep) and S (Sdep) deposition can also influence carbon–water relations. High levels of Sdep in the 1980s have been hypothesised to have led to stomata closure[10], thus current Sdep reductions might promote an increase in g[11], counteracting the CO2-induced water-saving effect of stomatal closure, particularly under non-limiting moisture conditions. Conversely, there is considerable evidence that increased Ndep has a fertilization effect on A, thereby contributing to enhancing a ca or climate driven increase in tree iWUE[12,13] and forest C-sinks[14,15] in N-limited forests.

Atmospheric N represents an additional input of N for trees, particularly for N-limited forest ecosystems in the Northern hemisphere. An increase in Ndep, however, has been associated with an acceleration of N cycling, with increase in N losses from ecosystems (through nitrate leaching and denitrification), when N saturation is reached[16]. Stable nitrogen isotope composition in plant materials (δ15N) has been used as a proxy of changes in ecosystem N availability[17]. An increase in tree δ15N has been observed in studies along Ndep gradients[18,19] and soil N manipulation experiments[20,21]. By contrast, a decrease in tree δ15N is expected in the case of a reduction of soil N availability (e.g., under a reduction of Ndep[22] or in case of tree canopy retention of Ndep[23]), or when N supply is insufficient to meet the demand caused by CO2 fertilization effects on A[24].

Combining multiple isotopes in tree rings across climate and Ndep gradients gives a powerful tool to advance our understanding of spatial and temporal patterns of iWUE and its main drivers. In particular, stable oxygen isotope composition (δ18O) provides physiological information in addition to that derived from δ13C, as its variation in plant material depends on the δ18O of the source water and that of the leaf water, the latter affected by transpiration and g[25]. Moreover, including δ15N allows the detection of changes in ecosystem N availability due to changes in the N input from the atmosphere. Studies assessing the effect of Ndep and its interactions with other pollutants and climatic variables by using a multiple isotope approach are predominantly site[26] or species-specific[22]. Analyses at regional scales, and involving different tree species, are paramount to achieving a better understanding of forest functioning in response to global changes.

We measured δ13Cw, δ18Ow and δ15Nw in tree rings from 1980 to 2010 for four of the most common species in Britain, i.e., Scots pine (Pinus sylvestris L.), Sitka spruce (Picea sitchensis Bong. Carr.), pedunculate oak (Quercus robur L.) and European beech (Fagus sylvatica L.) in 12 managed forests. Sites were selected along a gradient of climate and Ndep, while ensuring that within species soil type, forest structure, stand age and management remained similar (Table 1, Fig. 1). Specific goals were to: 1—document the temporal changes and spatial differences in the main climate parameters and Ndep and Sdep along these South–North gradient in Britain; 2—explore the temporal trends in iWUE, ci/ca and oxygen isotope discrimination, Δ18Ow (i.e., difference between δ18Ow and precipitation δ18O, see “Methods”) at the 12 forests and assess the possible species-specific physiological mechanisms (changes in A and/or g) underlying variations in iWUE; 3—evaluate whether sites receiving high Ndep experience increase in ecosystem N availability and N saturation (using δ15Nw as a proxy); 4—elucidate the drivers of spatial and temporal changes in the isotope-derived physiological and ecological processes.

Table 1

Description of the forest stands included in the study.

Site	Species	Lat	Long	Elevation (m asl)	Planting year	BA (m2 ha−1)	LAI (m2 m−2)	Soil type	Organic layer C:N	Top soil pH	Pa (mm)	Ta (°C)
Alice holt	Fagus sylvatica	51.1500	− 0.8600	80	1930	34.0	4.17	Eutric Planosol	21	5.1	662	9.7
Thetford	Fagus sylvatica	52.4100	0.8700	20	1930	31.0	3.78	Ferralic arenosol	19	3.9	528	9.7
Covert Wood	Fagus sylvatica	51.2000	1.1200	20	1950	31.0	5.15	Sceletic Leptosol	20	7.5	527	10.5
Shobdon	Fagus sylvatica	52.2400	− 3.0300	200	1952	31.4	5.49	Cambisol	18	4.1	714	9.6
Tummel	Picea sitchensis	56.7200	− 4.0500	400	1969	87.9	8.5	Ferric Podzol	37	5.1	796	8.1
Goyt	Picea sitchensis	53.2900	− 1.9800	275	1981	59.2	n.a.	Cabric Podzol	22	4.3	697	9.0
Rannoch	Pinus sylvestris	56.6500	− 0.4200	470	1965	46.6	6.13	Gleyic Podzol	38	4.2	1161	8.5
Ladybower	Pinus sylvestris	53.4100	1.7500	275	1952	49.2	2.93	Cabric Podzol	31	4.1	697	9.0
Rogate	Pinus sylvestris	51.0200	− 0.8700	80	1950	37.8 (#)	n.a.	Humic Podzol	25	3.1	662	9.7
Thetford	Pinus sylvestris	52.4100	0.8700	20	1967	36.0	4.08	Ferralic Arenosol	25	5.3	528	9.7
Alice Holt	Quercus robur	51.1500	− 0.8600	80	1935	24.2	6.1	Eutric Vertisol	19	5.5	662	9.7
Savernake	Quercus robur	51.5900	− 1.9200	107	1950	23.1	n.a.	Eutric Vertisol	38	5.2	796	8.1

Main site parameters for the investigated forest stands are given, including Latitude (Lat) and Longitude (Long), stand (planting year, basal area-BA, Leaf Area Index-LAI), soil-related (soil type, C: N ratio and pH) and climate-related (mean of annual precipitation, Pa, temperature, Ta calculated over the investigated years, i.e., 1980–2010) parameters. Basal area was measured in 2010 with the exclusion of Rogate (#), which was measured in 2018. Soil type was classified according to the WRB, 2015[65]. The same climate information was used for Rogate and Alice Holt, as the two sites are very close (within 19 km). This was also the case for Goyt and Ladybower, which are 30 km apart. Note that sites were grouped by tree species to better identify the pairing of sites according to similar age and soil type, but the contrasting levels of Ndep (low vs. high nitrogen deposition) are shown in Fig. 1 and Supplementary Table S1. na data not available.

Figure 1

Sites included in the study. Map showing forested sites along the precipitation gradient in Britain (panel A) and described in the Table 1. Size of the points reflects the levels of Ndep, which are reported in the panel B for each forest stand and species, together with the Sdep (panel C). Each point is the average (± standard deviation) of annual values (kg N ha−1 year−1 and kg S ha−1 year−1) across 1995–2010 for most of the forest stands, with the exception of Rogate (data from only 2010) and Shobdon and Covert Wood (modelled data, Ref. “Methods”). Note that for Goyt, we considered atmospheric deposition data collected at Ladybower, as the two sites are only 30 km apart. Black solid and dashed lines in the panel B indicate the low and high boundaries for the UK critical load for nitrogen deposition.

Results

Spatio-temporal changes in climate and atmospheric deposition across Britain

Across the investigated sites and during the studied years (1980–2010), mean temperature during the growing season significantly increased by approx. 0.02–0.06 °C year−1. Whereas, vapour pressure deficit (VPD) did not show a significant trend, except for an increase at the South-eastern most site of Covert Wood (Table 2). The standardised precipitation-evapotranspiration index (SPEI) relative to August, with three months time-scale (SPEI8_3) increased at seven of the sites, but not at Alice Holt, Rogate, Covert Wood, Goyt and Ladybower (Table 2).

Table 2

Temporal changes in the main environmental parameters at the investigated sites.

Site	Forest stand	Parameter
Name	Species	NH4 –N	NO3-N	Total Ndep	Total Sdep	SPEI	Tgrs	Tmaxgrs	Pgrs	VPDgrs	Ta
		kg (N or S) ha−1 year−2				year−1	°C year−1		cm year−1	kPa year−1	°C year−1
Covert Wood	Beech	n.a	n.a	n.a	n.a	0.03 (0.03)	0.06 (0.01)***	0.08 (0.03)*	0.08 (0.29)	0.008 (0.003)*	0.05 (0.02)*
Alice Holt	Beech	− 0.11 (0.07)	− 0.05 (0.02)	− 0.16 (0.07)*	− 0.28 (0.12)*	0.02 (0.02)	0.02 (0.009)*	0.006 (0.016)	0.008 (0.251)	0.007 (0.003)	0.02 (0.009)*
Alice Holt	Oak
Rogate	Scots pine
Savernake	Oak	− 0.4 (0.12)*	− 0.09 (0.04)*	− 0.49 (0.14)**	− 0.2 (0.2)	0.06 (0.02)*	0.03 (0.016)	0.04 (0.03)	0.29 (0.44)	− 0.002 (0.006)	0.03 (0.01)
Thetford	Beech	− 0.29 (0.10)*	− 0.08 (0.03)*	− 0.37 (0.12)**	− 0.36 (0.08)***	0.05 (0.02)*	0.02 (0.008)*	0.03 (0.01)*	0.64 (0.25)*	0.0008	0.02 (0.009)
Thetford	Scots pine
Shobdon	Beech	n.a	n.a	n.a	n.a	0.06 (0.02)*	0.03 (0.013)*	0.03 (0.02)	0.28 (0.44)	0.005 (0.003)	0.03 (0.01)
Goyt	Sitka spruce	− 0.21 (0.10)	− 0.11 (0.03)**	− 0.31 (0.25)	− 0.52 (0.19)*	0.05 (0.03)	0.03 (0.14)*	0.03 (0.02)	0.68 (0.48)	0.006 (0.004)	− 0.07 (0.01)***
Ladybower	Scots pine
Tummel	Sitka spruce	− 0.15 (0.07)	− 0.06 (0.05)	− 0.21 (0.10)	− 0.09 (0.11)	0.08 (0.03)*	0.02 (0.01)*	0.03 (0.02)	0.49 (0.43)	− 0.0009 (0.003)	0.02 (0.01)
Rannoch	Scots pine	− 0.15 (0.06)*	− 0.09 (0.03)*	− 0.24 (0.08)*	− 0.3 (0.07)***	0.08 (0.03)*	0.03 (0.010)*	0.02 (0.02)	0.52 (0.53)	0.004 (0.003)	0.02 (0.01)

Values of slope and standard error (in brackets) from linear regression analyses against year of wet NH4-N and NO3-N, total wet N (NH4-N + NO3-N) and S deposition, (SO4-S), growing season (grs) mean temperature and mean of maximum temperature (Tgrs and Tmaxgrs), vapour pressure deficit (VPDgrs) and precipitation (Pgrs), and mean annual temperature (Ta) for the 12 forest stands across Britain. The same deposition and climate information were used for Rogate and Alice Holt, as the two sites are very close (within 19 km). This was also the case for Goyt and Ladybower, which are 30 km apart. Slopes significantly different from zero were indicated by stars, according to the p-values: *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001. na long-term data not available.

Principal components analysis (PCA) using across sites long-term averages of environmental variables (PCA_spatial, hereinafter referred to as PCA_s) and within sites annual values (PCA_annual, hereinafter referred to as PCA_a) were conducted to reduce redundancy in the predictors along the climatic gradient. In the case of PCA_s, more than 90% of the variance was explained with two principal components (PC_s). PC_s1 was a general description of site climatic conditions, with highly significant effects from all variables (temperature, VPD and precipitation for the growing season and for the entire year). Strongest (negative) correlations with PC_s1 were with mean and maximum growing season temperature, suggesting that this axis primarily represented a South–North gradient. The PC_s2 correlated most strongly with VPD differences among sites, likely representing the East–West gradient. For the PCA_a, four PCs (PC_a) were required to explain 78% of the total variance. PC_a1 correlated significantly with temporal variability in VPD (negative) and precipitation (positive), PC_a2 with temporal variability in temperature (negative correlations), PC_a3 with growing season (negative) and annual (positive) SPEI indices and VPD (negative), while PC_a4 represented primarily temporal variability in precipitation and VPD (Fig. S1).

Over the 30 year period of the study, the mean annual ca increased by 51 ppm (from 339 to 390 ppm, slope = 1.7 ± 3.9 ppm year−1)[27]. Wet Ndep and wet Sdep decreased between 1995 and 2010 at all sites but Savernake and Tummel for Sdep, and Goyt/Ladybower and Tummel for Ndep (Table 2). Despite the general decreasing trends in wet Sdep and Ndep, half of the investigated sites are above the UK’s critical load for N to woodlands of 10–12 kg N ha−1 year−1;[28] (Fig. 1B).

Species-specific changes in iWUE

Relative changes in iWUE (i.e., value at 2010 minus value at 1980 divided by value at 1980) increased for Scots pine, oak and beech by 17.8 (± 5.6) %, 22.5 (± 11.3) % and 15.3 (± 7.0) %, respectively at all the sites. For Sitka spruce we found instead that iWUE decreased by 22.5 (± 9.9) %, with especially strong reductions in the youngest stand at Goyt that was planted in 1981, when tree-ring δ13C time series began (Table S2).

Trends in iWUE were not consistent across the investigated species. We found that iWUE increased for Scots pine and oak at all the sites (4 and 2 sites, respectively) over the recent 30 years (Fig. 2A, Table S2). Three of the four beech sites (Covert Wood, Shobdon and Thetford) showed no significant changes in iWUE, while a significant increasing trend in iWUE was observed for beech trees at Alice Holt. In contrast to the other species, the two Sitka spruce stands at Goyt and Tummel showed a reduction in iWUE, particularly for the youngest stand at Goyt (Fig. 2, Table S2).

Figure 2

Trends over time in the main physiological and ecological parameters as obtained from the measured tree-ring stable isotopes at the 12 forest stands across Britain. Changes in iWUE (Panel A), ci/ca ratio (Panel B), Δ18Ow (difference between δ18Ow and δ18O of precipitation, panel C) and δ15Nw (Panel D) for four species over the period 1980 to 2010. Each point represents the parameters derived from—or the actual—isotope ratios measured on wood materials pooled from 10 trees for a given year; i.e., δ13Cw—derived iWUE and ci/ca, δ18Ow—derived Δ18Ow and measured δ15Nw. Variability among the 10 trees per species and site was assessed for the year 2007 and data are reported in the Table S5. Regression lines are given for each site grouped by Ndep levels (i.e., low Ndep, light blue; high Ndep, black).

We tested whether the observed trends could be partially explained by changes in stand parameters i.e., diameter at breast height (DBH), mean height, basal area (BA) for six stands for which data was available (no data were available for the youngest stand at Goyt). Overall we did not find significant relationships between iWUE and stand-related parameters, with the exception of the Scots pine at Ladybower, where iWUE was positively correlated to changes in BA (Supplementary text S1).

Species-specific trends in ci/ca ratio and Δ18Ow

Similarly to what has been reported already for iWUE, changes in the ci/ca ratio were not consistent across the species and sites. For both beech and Scots pine at three sites, the parameter did not significantly change over the 30 years (Fig. 2B). The ci/ca ratio increased for the Sitka spruce at Goyt (slope = 0.008, p < 0.001) and Tummel (slope = 0.003, p < 0.001), the oak at Savernake (slope = 0.0007, p < 0.05), while it decreased for the oak at Alice Holt (slope = − 0.002, p < 0.01) and the Scots pine at Rannoch (slope = − 0.0009, p < 0.05).

Trends in Δ18Ow were not significant, except for the decrease for beech in Shobdon (slope = − 0.03 ± 0.013‰ year−1, p < 0.05) and for Scots pine at Rannoch (slope = − 0.021 ± 0.009‰ year−1, p < 0.05) and the increase for Scots pine at Ladybower (slope = 0.026 ± 0.007‰ year−1, p < 0.01) (Fig. 2C). Changes in Δ18Ow were not correlated with stand-related parameters, except for Rannoch Scots pine, where we observed a significant negative relationship between Δ18Ow and mean tree height and BA (Supplementary test S1).

Spatial and temporal changes in δ15Nw at different levels of Ndep

Differences in tree ring δ15Nw values between sites with Ndep below (low Ndep) and above (high Ndep) the critical loads (Fig. 1B) did not show a distinctive pattern. In the case of beech, the δ15Nw values were more positive (on average 1.57‰) at the high Ndep than at the low Ndep sites (Supplementary Table S3). However, we found the opposite for the oak, Scots pine and Sitka spruce, with trees at higher Ndep showing more negative δ15Nw values (Supplementary Table S3).

No significant trends in δ15Nw were observed (Fig. 2D), except for oak stands at Alice Holt (slope = 0.03‰ year−1, p < 0.05), beech stands at Thetford (slope = 0.04‰ year−1, p < 0.01) and Sitka spruce stand at Tummel (slope = 0.08‰ year−1, p < 0.01). A significant positive relationship between δ15Nw and wood %N was observed for 5 sites (Fig. S2), including two beech, one oak and two Scots pine stands, but no relationship was observed in the case of the two Sitka spruce stands.

Contribution of climate and anthropogenic factors on physiological and ecological processes

One of our goals was to elucidate the relative contribution of climate and atmospheric drivers on changes in physiological and ecological parameters obtained from tree-ring stable isotopes. Temperature increase and precipitation decrease (the dominant parameters in PC_s1) enhanced iWUE along the N-S gradient (Fig. 3A–C, Table 3). Whereas temporal increase in iWUE was mostly associated with changes in parameters related to moisture supply and demand conditions (VPD, SPEI and Precipitation, i.e., the dominant variables in PC_a1 and PC_a3) (Fig. 3B–D, Table 3). However, the high correlations among individual variables in many axes both in the PCA_s and the PCA_a precluded the possibility of drawing definitive conclusions about individual climatic drivers. Nonetheless, neither the increase in ca nor the changes in Ndep and Sdep were significant predictors in the model (Table 3). The inclusion of stand-related parameters (age and soil type) did not improve the statistical model fit (Supplementary text S2).

Figure 3

Spatial and temporal changes in iWUE vs. climate. Relationship between iWUE (mean through 1980–2010 for each tree species) and PC_s1 (panels A–C) and between year-by-year iWUE vs. PCA_a1 (panels B–D) across the investigated forest stands as obtained from the linear mixed effect model analyses. For the sake of clarity, we included the results of PCA analyses (panels A and B), to identify the climate factors that mostly affect the changes in iWUE as shown in the panels C and D. Slope and standard error values are reported in the Table 4. We used different symbols and colours for regression lines for each species, but slope presented in Table 4 refer to all observations together. The four investigated species were indicated with different colours (see legend in the panel C for details), while each site was indicated with a different symbol (see main legend for details). In panels A and B, ‘a’ and ‘grs’ indicate mean annual and growing season T, P and VPD, while 'maxgrs' refers to mean of maximum T. As for the SPEI, 8 and 12 indicate values for the month of August with a time-scale of one (8_1), two (8_2) and three (8_3) months and for the month of December.

Table 3

Relationship between physiological and ecological parameters and environmental factors.

Fixed effects	Estimate ± SE	Auto-correlation	Random effect		Rm2	Rc2
		Structure	Standard deviation	Residual variance
(i) All sites—including ca, sNdep and sSdep (n = 251)
Intercept	80.32 ± 2.30***	0.76	3.82	7.15	0.38	0.52
Conifer versus deciduous	8.26 ± 3.40*
PCA_s1	− 2.65 ± 0.67*
PCA_a1	− 0.65 ± 0.12***
PCA_a3	− 0.40 ± 0.19*
(ii) Model (i) without Sitka spruce (n = 210)
Intercept	80.75 ± 2.42***	0.50	5.21	4.93	0.41	0.72
Conifer versus deciduous	7.73 ± 3.85 ¥
PCA_s1	− 2.19 ± 0.79*
PCA_a1	− 0.64 ± 0.14***
ca	0.20 ± 0.04***
(iii) All parameters (n = 140)
Intercept	80.94 ± 2.33***	0.54	3.62	4.78	0.66	0.78
Conifer versus deciduous	9.81 ± 3.20*
PCA_s1	− 4.88 ± 1.12**
PCA_a1	− 0.67 ± 0.17***
PCA_a2	− 0.11 ± 0.22 n.s.
PCA_a3	− 0.46 ± 0.25 ¥
sNdep	− 0.59 ± 0.65 n.s.
sSdep	− 1.15 ± 0.45*
ca	0.04 ± 0.09 n.s.
aSdep	− 0.20 ± 0.19 n.s.
aNdep	0.11 ± 0.18 n.s.
((iv) Model (iii) without the Sitka spruce stands (n = 115)
Intercept	81.94 ± 1.86***	0.47	2.85	4.40	0.71	0.80
Conifer versus deciduous	8.51 ± 2.66*
PCA_s1	− 3.67 ± 0.63**
PCA_a1	− 0.66 ± 0.17***
PCA_a3	− 0.52 ± 0.25*
sSdep	− 0.47 ± 0.18*
aSdep	− 1.77 ± 0.40*

Statistics of the linear mixed effects models for the regression of iWUE as a function of the grouping into conifers and deciduous species (plant functional type, PFT), site (PCA_s) and time (PCA_a) climate variables from PCA analysis, changes in atmospheric CO2 (ca), spatial (sSdep, sNdep) and annual (aSdep, aNdep) changes in sulphur and nitrogen deposition. Output is given separately for the different model tested: (i) all 12 sites including climate, ca and only sSdep and sNdep; (ii) model (i) excluding the Sitka spruce stands; (iii) all parameters (including aSdep and aNdep, which are available only at 10 sites) and (iv) as model in (iii) but excluding Sitka spruce stands. Estimate of slope and intercept and standard error (SE) values are provided for each of the fixed factors. We considered random intercept for the site × species combination. Marginal (only fixed factors) and conditional (fixed + random factors) proportions of the explained variance are indicated as R2m and R2c, respectively. For each model, the number of observations is indicated with ‘n’. Note that we report results from the final model (Ref. “Methods”), but all the different models tested are listed in the Table S7 and results reported in the supplementary text S4. All independent variables were centered prior to analysis. Significance levels for fixed factors referring to simultaneous tests for general linear hypothesis testing are indicated as follows: ¥, p ≤ 0.10; *, p ≤ 0.05; **, p ≤ 0.01; ***, p ≤ 0.001.

Table 4

Relationship between physiological and ecological parameters and environmental factors.

Fixed effects	Estimate ± SE	Auto-correlation	Random effect		Rm2	Rc2
		Structure	Standard deviation	Residual variance
Δ18Ow—(i) All sites including ca, sSdep, sNdep (n = 251)
Intercept	32.88 ± 0.18***	0.25	0.38	0.41	0.31	0.63
Conifer versus deciduous	− 0.06 ± 0.27 n.s.
PCA_s1	− 0.12 ± 0.05*
PCA_a1	− 0.05 ± 0.013***
PCA_a2	0.11 ± 0.016***
PCA_a3	0.05 ± 0.02**
sSdep	− 0.03 ± 0.03 n.s.
Δ18Ow—(ii) All parameters, 10 sites (n = 140)
Intercept	33.07 ± 0.08***	0.12	0.11	0.35	0.72	0.75
Conifer versus deciduous	− 0.05 ± 0.01 n.s.
PCA_s1	− 0.24 ± 0.02***
PCA_a1	− 0.06 ± 0.01***
PCA_a2	0.12 ± 0.02***
PCA_a3	0.05 ± 0.11*
sSdep	− 0.10 ± 0.01***
aSdep	−  0.03 ± 0.01*
Δ18Ow—(iii) All parameters and sites without Sitka spruce (n = 210)
Intercept	33.03 ± 0.19***	0.27	0.45	0.44	0.10	0.56
Conifer versus deciduous	− 0.16 ± 0.30 n.s.
PCA_a1	− 0.05 ± 0.01***
PCA_a2	0.10 ± 0.02***
Δ18Ow—Climate and Ndep 10 sites (n = 204)
Intercept	− 1.25 ± 0.27***	0.14	0.46	0.82	0.43	0.57
Conifer versus deciduous	− 0.04 ± 0.37 n.s.
sNdep	0.14 ± 0.07 ¥
PCA_s1	0.51 ± 0.13**

Statistics of the linear mixed effects models for the regression of for Δ18Ow and δ15Nw as a function of the grouping into conifers and deciduous species (plant functional type, PFT), site (PCA_s) and time (PCA_a) climate variables from PCA analysis, spatial (sSdep, sNdep) and annual (aSdep, aNdep) changes in sulphur and nitrogen deposition. For Δ18Ow output is given separately for the different model tested (Ref. “Methods”): (i) all 12 sites including climate, ca and only sSdep and sNdep; (ii) all parameters (including aSdep and aNdep, which are available only at 10 sites) and (iii) models (i–ii) without Sitka spruce stands (both leading to the iii) as final model (Ref. Table S7 and Supplementary text S5). Estimate of slope and intercept and standard error (SE) values are provided for each of the fixed factors. We considered random intercept for the site × species combination. Marginal (only fixed factors) and conditional (fixed + random factors) proportions of the explained variance are indicated as R2m and R2c, respectively. For each model, the number of observations is indicated with ‘n’. Note that we report results from the final model, but all the different models tested are listed in the Table S7 and results reported in the supplementary text S5 and S6. All independent variables were centered prior to analysis. Significance levels for fixed factors referring to significance in simultaneous tests for general linear hypothesis testing are indicated as follows: ¥, p ≤ 0.10; *, p ≤ 0.05; **, p ≤ 0.01; ***, p ≤ 0.001.

Given that the iWUE trend observed for the Sitka spruce stands could be partially related to the age effect, with particular reference to the youngest one at Goyt, we repeated the analyses excluding Sitka spruce. Results confirmed no relationship between Ndep and iWUE (Table 3). We found, however, a negative relationship between both sSdep and aSdep and iWUE (Table 3) and the expected significant relationship between ca and iWUE (Table 3). Finally, when performing the analysis by level of Ndep (i.e., high and low Ndep for sites above and below the critical loads, respectively) and by including or not the Sitka spruce stands, the positive relationship between iWUE and ca turned out to be significant, with the exception of the analysis for the high Ndep sites and when the Sitka spruce stand at Goyt was included (Supplementary Table S4).

Differences in both spatial and temporal changes in climatic conditions significantly affected Δ18Ow. Three axes from the PCA_a analyses (PC_a1, PC_a2 and PC_a3) were significant, suggesting a combined effect of temperature and moisture variability on changes in Δ18Ow (Table 4). When considering the subset of sites where annual atmospheric deposition data were available, we did find a significant and negative relationship between Δ18Ow and sSdep and aSdep. This, however, was not the case when Sitka spruce stands were removed from the analyses, to account for possible age effects on changes in Δ18Ow (Table 4, Supplementary text S5).

When considering all sites and species together, at high and low Ndep sites, a clear pattern for δ15Nw in tree rings was not seen and no differences were observed when grouping species by plant functional type (Table 4). Changes in δ15Nw were explained primarily by climate factors, particularly temperature and precipitation, with a marginal (p = 0.1) effect of differences across sites in Ndep (Table 4). Soil type and stand age did not appear significant predictors for δ15Nw in the statistical model (Supplementary text S3).

The combined influences of temporal changes and spatial differences in climatic and atmospheric deposition variables across all three isotopes were examined with the help of a mixed-effect path model (Fig. 4A). Consistent with the earlier analysis using linear mixed models, the first axis of PC_s significantly affected all three isotopes, negatively for iWUE and Δ18Ow and positively for δ15Nw. This is consistent with the first axis of PC_s being related to a South–North gradient (i.e., positively with precipitation and negatively with temperature and VPD). Conversely, the first and third axis of PC_a only affected iWUE (both effects negative) and Δ18Ow (negative effect for PC_a1 and positive effect for PC_a3). This is also consistent with PC_a1 being positively related to precipitation and SPEI and negatively related to VPD and temperature. Ndep was only marginally (p < 0.1) important for δ15Nw. A significant positive effect of Δ18Ow on iWUE was also found, while no relationship was observed between iWUE and δ15Nw (Fig. 4A). Similar results were found when excluding Sitka spruce stands from the analyses for both Δ18Ow and δ15Nw, in addition to the significant positive effect of ca on iWUE and (Fig. 4B). When including in the analyses only sites where annual atmospheric deposition data were available, we also found a negative effect of sSdep and aSdep on iWUE (Fig. S5).

Figure 4

Direct and indirect effects on tree-ring stable isotopes and derived physiological parameters. Result from the structural equation modelling analysis, which included all sites (panel A) and all sites excluding the two Sitka spruce stands (panel B). Details of equations are reported in the Table S7. Continuous arrows indicate relationships significant at the p ≤ 0.05 or greater (depending on number of stars, i.e., *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001). Dashed arrows (and the ‘¥’ symbol) indicate relationships significant at p ≤ 0.10. Thickness of lines reflects level of significance. Notice that degrees of freedom vary between site-based and time-based analyses. Black and blue arrows indicate negative and positive relationships, respectively. Double-headed arrow indicates correlated errors between variables. PC_s1 indicates the first component from the PCA_s (using long-term average climate variables across sites) and PC_a1, PC_2 and PC_a3 the first three principal components from the PCA_a (using annual climate parameters). Number next to the paths indicate standardized path coefficients. Numbers below each isotope-related parameters indicate marginal and conditional R2, respectively. They are not always the same as those reported in the Tables 3 and 4 due to slightly different equations considered in the SEM analyses. In particular, PFT was not included as fixed factor in any model, while Δ18Ow was included as fixed factor in the model for iWUE (Ref. Table S7 for all details).

Discussion

Climate and atmospheric deposition changes across the investigated sites

The first goal of this study was to document the changes in climate and atmospheric deposition across the investigated sites over the period 1980–2010, so as to better assess their effects on iWUE and other isotope-related parameters. Trees at most of the investigated sites experienced increasingly wetter and warmer growing conditions, particularly in Scotland, corroborating general trends for Britain[29]. PC analyses showed that temperature and precipitation are the main climate factors describing the South–North gradient, while parameters related moisture conditions (SPEI, VPD and precipitation) were more relevant in long-term climate variations.

While ca has continued to increase by 14% from 2010 to 1980, differences were observed in magnitude and directionality of changes in atmospheric deposition. A significant reduction of total wet Sdep and Ndep (NH4-N + NO3-N) was also observed at most sites[1,30,31], with steeper slopes for Sdep than Ndep. Decreasing Sdep is consistent with national scale reports, whereby emissions and total Sdep have greatly declined over the last few decades[4]. However, a significant decline in total Ndep (including oxidised and reduced form) has not yet been reported nationally despite reduction of N emissions of both nitrogen oxides and ammonia[4]. Taken all together, these results suggest that directionality and magnitude of changes in Ndep and Sdep across the environmental gradient may modulate tree species response to ca, as we discuss in the following paragraphs.

Tree species differed in iWUE trends and underlying mechanisms

Our second goal was to assess trends in iWUE and elucidate potential physiological mechanisms underlying changes in iWUE. Three of the four investigated tree species showed that iWUE increased over 30-years, corroborating global tree-ring δ13C based analyses[12,13,32]. Across sites, the relationship between iWUE and Δ18Ow (Fig. 4) suggests that overall there is a tight coupling between A and g across the investigated species. However, the within-sites relationships were significant only at 5 of the 12 forest stands (Table S2), suggesting physiological strategies differ amongst—and likely also within—tree species.

Scots pine and oak showed an increase in iWUE, but the underlying mechanisms were likely different. Scots pine at southernmost sites (Rogate and Thetford) showed no changes in transpiration losses (no significant trend for Δ18Ow) and a proportional regulation of A and g which contributed to a constant ci/ca. Scots pine has a very conservative strategy regarding water use, due to its tight stomatal control under moisture limitation[33]. Whereas, Δ18Ow for Scots pine at the northernmost sites (Ladybower and Rannoch) significantly changed, though in the opposite direction. The increase in Δ18Ow for Scots pine at Ladybower would indicate a reduction in gs, which is supported by the positive relationship between iWUE and Δ18Ow (Table S2). The site showed the highest mean Sdep values among all sites. Even though trend in Sdep has significantly decreased at the site, we cannot exclude a possible legacy (negative) effect of deposition on gs[23]. In the case of Rannoch, the physiological signal recorded in the reduction in Δ18Ow could be partially confounded by tree stand development (as a negative relationship was found between Δ18Ow and tree height and BA) or the significant reduction in Sdep, which would have a positive effect on gs[10].

In contrast, the increase in iWUE and the significant changes in the ci/ca ratio (though in opposite direction) for the two oak stands was mostly related to a more dynamic adjustment of leaf gas exchanges to environmental conditions, with likely reductions in g and lower A, leading to a lower increase in iWUE in the case of Savernake (where both ci/ca and Δ18Ow increased) compared to Alice Holt (no changes in g, as suggested by no significant trend in Δ18Ow, while ci/ca decreased). This result indicated that Sdep is still negatively affecting oak trees at Savernake—one of the few sites where Sdep has not significantly declined over the last decades. Consistent with our results, a previous study found that exposure to high SO2 pollution significantly reduced gs leading to an increase in iWUE for oak trees in southern England[10].

The reduction in iWUE in the fast-growing[34] Sitka spruce trees, particularly the youngest stand at Goyt), most likely reflect age or management history effects. Height growth during tree development may affect δ13C and hence iWUE estimates via changes in light availability[35] and hydraulic conductance[36], but also through increasing LAI as canopy develops and contributes to respired CO2[37]. Note, however, that the trend in iWUE at Goyt cannot be explained by increasing hydraulic constraints with tree height[36], which would be in the opposite direction (i.e., an increase in iWUE in taller trees). It is most likely that the observed trends may reflect common historical upland conifer establishment practices, e.g., fertilisation on tree establishment until approximately 10 years of age and thinning between 15 and 20 years. This may have contributed to reduce tree competitions for light and soil water, thus accelerating growth but reducing iWUE in the early stage of development until canopy closure[38].

Beech was the least responsive species, as no changes in iWUE were observed in general, even though the four selected sites were along a clear precipitation gradient (Fig. 1). However, the wettest site (Shobdon) showed the lowest % increase in iWUE (only 5% compared to 16–19% changes at the remaining three sites), which could partially be explained by the increase in gs suggested by the reduction in Δ18Ow. Moreover, when comparing even-aged beech and oak stands at the same site (i.e., Alice Holt), with similar soil and climate conditions, we found that despite the two species showing similar iWUE, beech had a higher Δ18Ow than oak (Fig. S3). This suggests differences between the two species for gs, which could partially explain the significantly (p < 0.05) higher slope in beech compared to oak for temporal trends in Δ18Ow.

Changes in tree iWUE and Δ18Ow across Britain were associated with climate and atmospheric S deposition

Changes (over time) and differences (across sites) in moisture conditions and temperature along the North–South gradient in Britain were overall more important than atmospheric deposition in explaining variations in iWUE. Increase in temperature, both moving North–South (PC_s1) across Britain and over time (PC_a2), significantly increased iWUE. This could be explained by a positive effect of temperature on A, but also a reduction in gs, likely associated with the increase in VPD, the major component of PC_a3. Indeed, the latter was significantly and positively correlated with Δ18Ow. This finding is in agreement with results for two broadleaf species in the Northeastern US[26], whereby changes in iWUE were predominantly associated with soil moisture conditions and no effects of either Ndep or Sdep were observed.

The detection of atmospheric deposition and ca signals in our study, however, was possibly confounded by likely age effects and/or associated tree structural and functional changes (as we discussed above). Interestingly, when Sitka spruce stands were excluded from the analyses, the increase in iWUE was associated with rising ca, but also with the reduction in Sdep, while there was no relationship with Ndep. This is quite important, as Sdep and ca may have a similar effect on gs, i.e., an increase in Sdep and ca leads to a reduction in gs[10]. However, changes in Sdep and ca do not always follow the same direction. Indeed we reported a reduction in Sdep at most of the sites, under a general increase in ca. The negative relationship between Δ18Ow and Sdep suggests that the alleviation of the negative effects from Sdep on gs[22] may be stronger than the ca effect leading to stomatal closure (which does not find indication in our data). These results on one hand suggest that variations in Sdep outweigh changes in Ndep and ca. Sdep may influence tree water-use directly, by affecting gs[22] or indirectly, by leading to soil acidification, loss of calcium, which play a significant role in controlling gs[39]. On the other hand, they indicate that it is essential to account for changes caused by stand development and management history, if the goal is to disentangle climatic and anthropogenic drivers of change in iWUE.

No changes in ecosystem N availability due to Ndep

Given that some of the investigated forest stands have reached the critical load in terms of Ndep (Ref. Fig. 1B), we expected to find significant changes in ecosystem N due to possible N losses from the forest ecosystems. Yet, spatial and temporal variations in δ15Nw values did not indicate differences in ecosystem N dynamics between high and low Ndep levels. If this was the case, we should have observed an increase in tree-ring δ15Nw over time and moving from low to high Ndep sites, as a consequence of high nitrification and loss (through denitrification and/or leaching) of the 15 N-depleted NO3, leaving trees with 15 N-enriched NO3[40]. This was only observed for the beech at Thetford, which is already confirmed as an N saturated site by detailed gradient studies[41]. Repeated soil analysis suggests N accumulation in soil organic layers under the conifer species, including the Tummel Sitka spruce site over the last 15 years, but no signs of N saturation[41]. However, a significant positive relationship has been observed between Ndep and NO3 leaching across some conifer sites in this study, so further N input to conifer forests could cause significantly higher NO3 leaching[41].

Our tree-ring δ15Nw data suggest that N availability has not changed for most of the forest stands over the investigated 30-years period, which is in contrast to declining δ15N trends reported both in the USA[22,42] and globally[43], whereby studies suggest on going oligotrophication. However, we cannot exclude that N recycling processes within trees could also have contributed to the observed lack of a trend in δ15Nw[44], particularly in the case where no relationship was found between δ15Nw and wood %N (Fig. S1).

The strong climate signal explaining spatial variations in the tree-ring δ15N values is consistent with a global analysis[45], which identified temperature and precipitation as the main climatic drivers of changes in the foliar δ15Nw across different forest ecosystems. An indirect effect of wet Ndep-via climate differences- on tree-rings δ15Nw cannot be fully excluded. Indeed, the significant negative relationship between atmospheric deposition and PC_s1, partially indicating that wet Ndep is strongly related to the amount of precipitation. The predominantly westerly airflow across the UK brings less polluted air from the Atlantic. Orographic cloud formation in the more mountainous regions of the NW however leads to a substantially higher rainfall and consequently higher Ndep[46].

Our results also indicate that pinpointing variations in δ15Nw caused by gradual changes in ambient Ndep is more challenging than in N manipulation experiments, where abrupt and high N doses over relatively short periods of time predominate. Moreover, the lack of information regarding both temporal variability in dry Ndep and changes in the N isotopic signatures for NH4 and NO3 might increase the uncertainties in the detection of changes in atmospheric Ndep by using δ15N in tree rings[47]. Monitoring changes in the isotopic composition of N-specific compounds in rainfall over time can greatly improve our ability to use δ15N in tree rings to detect changes in N input from Ndep.

Conclusions

Long-term forest monitoring systems, such as the Level II—ICP forests programme, provide unique near-natural systems for assessing the effect of climate change on ecophysiological responses of different tree-species at a regional scale and elucidating interactions among environmental forcing factors and forest ecosystem response. Our results showed that the increase in iWUE was not uniform across sites and that species-specific underlying physiological mechanisms were likely affected by the interactions between climate and atmospheric drivers (for oak and Scots pine), but also an tree structural changes during stand development (for Sitka spruce).

Spatial and temporal changes in temperature and moisture conditions overrode the effect of atmospheric deposition and ca on changes in iWUE for the investigated forests in Britain. This is remarkable since such an increase is widely predicted to occur in response to increasing ca. Our results suggest that the effect of increasing ca on temporal changes in iWUE could be over estimated, if concomitant changes in atmospheric deposition or ontogenic effects (such as structural changes during tree development) are not accounted for. The tree-ring δ15Nw analyses did not provide evidence for changes in N availability caused by changes in Ndep. In particular, sites receiving high Ndep (now considered above the critical loads) did not show evidence of N saturation, with the exception of the beech site at Thetford. Spatial differences in tree ring δ15Nw were mostly explained by differences across sites in temperature and precipitation, rather than changes in Ndep.

The multiple-species and regional analysis indicate that climate change may affect the most common native and introduced species in British woodlands. Lower summer rainfall and high temperature and VPD are likely to become more frequent in South-eastern Britain, thus affecting the future site suitability of beech woodland, as the species is more susceptible to drought[48], while oak and Scots pine could cope better (see e.g., results at Alice Holt). However, even though oak may be physiologically plastic in response to future climate change, widespread oak decline across Britain has been observed related to a number of biotic and abiotic factors, including climate and pollution[49]. Sitka spruce, the major upland planted timber species in the UK, native of coastal Northwestern America, is likely to maintain continued good growth in British northern uplands where water stress is less pronounced. However, the young ages of these stands and intensity and frequency of management interventions make it more difficult to disentangle development from environmental effects on iWUE trends.

Methods

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 (http://icp-forests.net/), 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 monitoring[30] 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: https://sac.csic.es/spei/.

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 described[30]. 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 fluxes[30]. 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 used[4]. 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 data[4]. However, only the total wet deposition was included in the analyses here, as we previously reported a good agreement between modelled and measured wet Ndep[50].

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: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:ca values were obtained from Mauna Loa records[27], and δ13Ca values were obtained from Mauna Loa records[5] 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: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]: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: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 (https://wateriso.utah.edu/waterisotopes/pages/data_access/oipc.html) 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 gs[56].

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 stem[57]. 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 it[57,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-centered[59]. 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 R[60]) 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 MuMIn[62].

Finally, to examine the joint effects of climatic conditions on Δ18Ow, δ15Nw and iWUE, a mixed-effect confirmatory path model was employed using package piecewiseSEM[63]. 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 relationship[63]. All statistical analyses were carried out inside RStudio 1.0.143 with R 3.4.0[64].

Supplementary file1 (PDF 1557 kb)