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.

. Description of the forest stands included in the study. 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, P a , temperature, T a 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 N dep (low vs. high nitrogen deposition) are shown in Fig. 1 and Supplementary Table S1. na data not available.

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 δ 13 C 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).
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 c i /c a ratio and Δ 18 o w .
Similarly to what has been reported already for iWUE, changes in the c i /c a 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 c i /c a 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 Δ 18 O w 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 Δ 18 O w were not correlated with stand-related parameters, except for Rannoch Scots pine, where we observed a significant negative relationship between Δ 18 O w and mean tree height and BA (Supplementary test S1).
Spatial and temporal changes in δ 15 n w at different levels of N dep . Differences in tree ring δ 15 N w values between sites with N dep below (low N dep ) and above (high N dep ) the critical loads (Fig. 1B) did not show a distinctive pattern. In the case of beech, the δ 15 N w values were more positive (on average 1.57‰) at the high N dep than at the low N dep sites (Supplementary Table S3). However, we found the opposite for the oak, Scots pine and Sitka spruce, with trees at higher N dep showing more negative δ 15 N w values (Supplementary Table S3).
No significant trends in δ 15 N w 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 Table 2. Temporal changes in the main environmental parameters at the investigated sites. Values of slope and standard error (in brackets) from linear regression analyses against year of wet NH 4 -N and NO 3 -N, total wet N (NH 4 -N + NO 3 -N) and S deposition, (SO 4 -S), growing season (grs) mean temperature and mean of maximum temperature (T grs and T maxgrs ), vapour pressure deficit (VPD grs ) and precipitation (P grs ), and mean annual temperature (T a ) 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.  (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 c a nor the changes in N dep and S dep 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). 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 N dep and iWUE (Table 3). We found, however, a negative relationship between both sS dep and aS dep and iWUE (Table 3) and the expected significant relationship between c a and iWUE (Table 3). Finally, when performing the analysis by level of N dep (i.e., high and low N dep for sites above and below the critical loads, respectively) and by including or not the Sitka spruce stands, the positive relationship between iWUE and c a turned out to be significant, with the exception of the analysis for the high N dep 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 Δ 18 O w . 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 Δ 18 O w (Table 4). When considering the subset of sites where annual atmospheric deposition data were available, we did find a significant and negative relationship between Δ 18 O w and sS dep and aS dep . This, however, was not the case when Sitka spruce stands were removed from the analyses, to account for possible age effects on changes in Δ 18 O w ( Table 4, Supplementary text S5).

Beech
Oak Scots pine Sitka spruce

Sites and species
Alice www.nature.com/scientificreports/ When considering all sites and species together, at high and low N dep sites, a clear pattern for δ 15 N w in tree rings was not seen and no differences were observed when grouping species by plant functional type (Table 4). Changes in δ 15 N w were explained primarily by climate factors, particularly temperature and precipitation, with a marginal (p = 0.1) effect of differences across sites in N dep (Table 4). Soil type and stand age did not appear significant predictors for δ 15 N w 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 Δ 18 O w and positively for δ 15 N w . 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 Δ 18 O w (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. N dep was only marginally (p < 0.1) important for δ 15 N w . A significant positive effect of Δ 18 O w on iWUE was also found, while no relationship was observed between iWUE and δ 15 N w (Fig. 4A). Similar results were found when excluding Sitka spruce stands from the analyses for both Δ 18 O w and δ 15 N w , in addition to the significant positive effect of c a 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 sS dep and aS dep on iWUE (Fig. S5). Temporal PC1 (30 %) Temporal PC3 (19 %) 6 3  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 timescale of one (8_1), two (8_2) and three (8_3) months and for the month of December.

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.  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. Scientific RepoRtS | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w www.nature.com/scientificreports/ While c a 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 S dep and N dep (NH 4 -N + NO 3 -N) was also observed at most sites 1,30,31 , with steeper slopes for S dep than N dep . Decreasing S dep is consistent with national scale reports, whereby emissions and total S dep have greatly declined over the last few decades 4 . However, a significant decline in total N dep (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 N dep and S dep across the environmental gradient may modulate tree species response to c a , as we discuss in the following paragraphs.

Fixed effects Estimate ± SE
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 δ 13 C based analyses 12,13,32 . Across sites, the relationship between iWUE and Δ 18 O w (Fig. 4) suggests that overall there is a tight coupling between A and g s 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.   www.nature.com/scientificreports/ 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 Δ 18 O w ) and a proportional regulation of A and g s, which contributed to a constant c i /c a . Scots pine has a very conservative strategy regarding water use, due to its tight stomatal control under moisture limitation 33 . Whereas, Δ 18 O w for Scots pine at the northernmost sites (Ladybower and Rannoch) significantly changed, though in the opposite direction. The increase in Δ 18 O w for Scots pine at Ladybower would indicate a reduction in g s , which is supported by the positive relationship between iWUE and Δ 18 O w (Table S2). The site showed the highest mean S dep values among all sites. Even though trend in S dep has significantly decreased at the site, we cannot exclude a possible legacy (negative) effect of deposition on g s 23 .
In the case of Rannoch, the physiological signal recorded in the reduction in Δ 18 O w could be partially confounded by tree stand development (as a negative relationship was found between Δ 18 O w and tree height and BA) or the significant reduction in S dep , which would have a positive effect on g s 10 . In contrast, the increase in iWUE and the significant changes in the c i /c a 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 s and lower A, leading to a lower increase in iWUE in the case of Savernake (where both c i /c a and Δ 18 O w increased) compared to Alice Holt (no changes in g s , as suggested by no significant trend in Δ 18 O w , while c i /c a decreased). This result indicated that S dep is still negatively affecting oak trees at Savernake-one of the few sites where S dep has not significantly declined over the last decades. Consistent with our results, a previous study found that exposure to high SO 2 pollution significantly reduced g s 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 δ 13 C 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 CO 2 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  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 R 2 , 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 Δ 18 O w was included as fixed factor in the model for iWUE (Ref. www.nature.com/scientificreports/ 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 g s suggested by the reduction in Δ 18 O w . Moreover, when comparing evenaged 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 Δ 18 O w than oak (Fig. S3). This suggests differences between the two species for g s , which could partially explain the significantly (p < 0.05) higher slope in beech compared to oak for temporal trends in Δ 18 O w .

Changes in tree iWUE and Δ 18 o w 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 g s , likely associated with the increase in VPD, the major component of PC_a3. Indeed, the latter was significantly and positively correlated with Δ 18 O w . 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 N dep or S dep were observed. The detection of atmospheric deposition and c a 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 c a , but also with the reduction in S dep , while there was no relationship with N dep . This is quite important, as S dep and c a may have a similar effect on g s , i.e., an increase in S dep and c a leads to a reduction in g s 10 . However, changes in S dep and c a do not always follow the same direction. Indeed we reported a reduction in S dep at most of the sites, under a general increase in c a . The negative relationship between Δ 18 O w and S dep suggests that the alleviation of the negative effects from S dep on g s 22 may be stronger than the c a effect leading to stomatal closure (which does not find indication in our data). These results on one hand suggest that variations in S dep outweigh changes in N dep and c a . S dep may influence tree water-use directly, by affecting g s 22 or indirectly, by leading to soil acidification, loss of calcium, which play a significant role in controlling g s 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 N dep . Given that some of the investigated forest stands
have reached the critical load in terms of N dep (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 δ 15 N w values did not indicate differences in ecosystem N dynamics between high and low N dep levels. If this was the case, we should have observed an increase in tree-ring δ 15 N w over time and moving from low to high N dep sites, as a consequence of high nitrification and loss (through denitrification and/or leaching) of the 15 N-depleted NO 3 , leaving trees with 15 N-enriched NO 3 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 N dep and NO 3 leaching across some conifer sites in this study, so further N input to conifer forests could cause significantly higher NO 3 leaching 41 .
Our tree-ring δ 15 N w 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 δ 15 N 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 δ 15 N w 44 , particularly in the case where no relationship was found between δ 15 N w and wood %N (Fig. S1).
The strong climate signal explaining spatial variations in the tree-ring δ 15 N values is consistent with a global analysis 45 , which identified temperature and precipitation as the main climatic drivers of changes in the foliar δ 15 N w across different forest ecosystems. An indirect effect of wet N dep -via climate differences-on tree-rings δ 15 N w cannot be fully excluded. Indeed, the significant negative relationship between atmospheric deposition and PC_s1, partially indicating that wet N dep 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 N dep 46 . Our results also indicate that pinpointing variations in δ 15 N w caused by gradual changes in ambient N dep 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 N dep and changes in the N isotopic signatures for NH 4 and NO 3 might increase the uncertainties in the detection of changes in atmospheric N dep by using δ 15 N 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 δ 15 N in tree rings to detect changes in N input from N dep .

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 Scientific RepoRtS | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w www.nature.com/scientificreports/ 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 c a 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 c a . Our results suggest that the effect of increasing c a 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 δ 15 N w analyses did not provide evidence for changes in N availability caused by changes in N dep . In particular, sites receiving high N dep (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 δ 15 N w were mostly explained by differences across sites in temperature and precipitation, rather than changes in N dep .
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-fores ts.net/), with the exception of Covert Wood, Shobdon and Goyt. The Goyt site was added as a high N dep site as a contrast to the low N dep 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 N dep , S dep 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 (T a ) and mean maximum (T amax ) values for temperature were calculated from monthly mean and maximum air temperature, T, respectively, and annual precipitation (P a ) 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 (N dep ) and sulphur deposition (S dep ) were obtained from measured bulk precipitation and throughfall water volumes at the sites and measured elemental concentrations (NO 3 − , NH 4 + and SO 2-4 ) as previously described 30 . For the spatial analyses, we considered mean of annual deposition (sN dep and sS dep ), obtained as the sum of total (NH 4 -N + NO 3 -N for N dep ) 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 N dep and S dep dataset was used 4 . The estimate included wet and dry NH x -N (NH 4 , NH 3 ), NO y -N (NO 2 , NO 3 , HNO 3 ) and S (SO x = SO 2 and SO 4 ) 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 N dep 50 . For the temporal analyses, only wet deposition (as calculated from NO 3 − , NH 4 + and SO 2-4 concentrations in bulk precipitation) was considered (indicated as aN dep and aS dep ), 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 Scientific RepoRtS | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w www.nature.com/scientificreports/ 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 δ 15 N w 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 δ 13 C w by VG Prism III Isotope ratio mass spectrometer at the School of Geosciences (University of Edinburgh, UK). For δ 15 N w , 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 δ 18 O w , 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 δ 13 C w , and at the Stable Isotope facility of the UC Davis, University of California (USA) for δ 18 O w and δ 15 N w . Stable isotope abundances are expressed as ratios of 13 C/ 12 C, 15 N/ 14 N and 18 O/ 16 O using δ-notation (in per-mil;‰) relative to international standards (VPDB for δ 13 C w , atmospheric N 2 for δ 15 N w and VSMOW for δ 18 O w ). The standard deviation for internal standards was 0.1‰ for δ 13 C w (PACS-2), 0.2, 0.3 and 0.4‰ for δ 18 O w (IAEA 600, IAEA 601 and IAEA 602, respectively), and between 0.1 and 0.3‰ for δ 15 N w (USGS-41 Glutamic Acid and peach leaves, respectively).

Calculation of iWUE and Δ 18 o w .
The iWUE and the c i /c a ratio were derived from measured δ 13 C w values, and based on the well-established theory linking leaf c i /c a with carbon isotope discrimination, Δ 13 C w , 51 as shown in the equation below: where δ 13 C a and δ 13 C w are the carbon isotope compositions of c a and wood, a is the isotope fractionation during CO 2 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 c i from the following equation: c a values were obtained from Mauna Loa records 27 , and δ 13 C a 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 CO 2 mol −1 H 2 O) was then calculated using the following equation: where 1.6 is the molar diffusivity ratio of CO 2 to H 2 O (i.e., g CO2 = g H2O /1.6). Note that in the Eqs. (2) and (3) we used average of values measured over growing season months (May-September) for both c a and δ 13 C a .
Tree-ring oxygen isotope discrimination, Δ 18 O w , was calculated according to Eq. (4) 53 : where δ 18 O w is the oxygen isotope composition we measured in each ring, while δ 18 O s is the oxygen isotope composition of the source water, i.e. the soil water, which we assumed to reflect the δ 18 O of precipitation (δ 18 O P ). We estimated annual values of δ 18 O P at each site as described by Barbour et al. 54 , by considering the following equation: (1)  18 O) 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 Δ 18 O results from an increase in transpiration and g s 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 δ 18 O and also δ 13 C 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 Δ 18 O w and parameters describing stand development (diameter at the breast height and height, Supplementary Table S2 text S1); (3) changes in iWUE and δ 15 N w 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, Δ 18 O w , c i /c a, δ 15 N w ) 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., aS dep and aN dep , to assess the effect of temporal changes) and the mean over the monitoring time (i.e., sS dep and sN dep ) to evaluate both the contribution of within-site temporal changes and cross-sites differences on changes in iWUE, c i /c a , Δ 18 O w and δ 15 N w . 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 CO 2 , N dep and S dep ) 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 (R 2 m and R 2 c , respectively) were calculated using package MuMIn 62 . Finally, to examine the joint effects of climatic conditions on Δ 18 O w , δ 15 N w 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 Δ 18 O w 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 (5) δ 18 O P = 0.52T a − 0.006T 2 a + 2.42P a − 1.43P 2 a − 0.046 √ E − 13.0