www.nature.com/scientificreports OPEN Climate and atmospheric deposition effects on forest water‑use efficiency and nitrogen availability across Britain Rossella Guerrieri 1,5*, Elena Vanguelova2, Rona Pitman2, Sue Benham 2 , Michael Perks3, James I. L. Morison 2 & Maurizio Mencuccini1,4 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. 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 C ­ O 2, sulphur (S) and reactive nitrogen (N). While sulphur emissions have been successfully regulated in Europe and North ­America1, N compounds, particularly ammonia and nitrous oxide, continue at relatively high levels compared to the pre-industrial p ­ eriod2–4. These changes in N and S atmospheric concentrations have occurred together with increasing atmospheric C ­ O2 concentration (­ ca), which recently exceeded 410 p ­ pm5, with a rela- tive 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 (gs)6, thus raising intrinsic water use efficiency of trees (iWUE = A/gs). Higher iWUE has been commonly reported for conifer compared to broadleaf ­species7. Moreover, the increase in iWUE derived from carbon isotope compo- sition (δ13C) in tree r­ ings8 has generally been associated with an active response of trees to increasing c­ a, 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 gs9. 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 ­closure10, thus current ­Sdep reductions might promote an increase in gs11, counteracting the C ­ O2-induced water-saving effect of stomatal closure, particularly under non-limiting moisture conditions. Conversely, there is considerable evidence that 1 Centre for Ecological Research and Forestry Applications, CREAF, c/o Universidad Autonoma de Barcelona, Edificio C, 08290 Cerdanyola, Barcelona, Spain. 2Forest Research, Alice Holt Lodge, Farnham, Surrey GU10 4LH, UK. 3Forest Research, Northern Research Station, Roslin EH25 9SY, Midlothian, Scotland, UK. 4ICREA, Barcelona, Spain. 5Present address: Department of Agricultural and Food Sciences, University of Bologna, 40127 Bologna, Italy. *email: rossellaguerrieri@gmail.com Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 1 Vol.:(0123456789) www.nature.com/scientificreports/ Elevation Planting LAI Organic Site Species Lat Long (m asl) year BA ­(m2 ­ha−1) ­(m2 m−2) Soil type layer C:N Top soil pH Pa (mm) Ta (°C) Fagus Eutric Alice holt 51.1500 − 0.8600 80 1930 34.0 4.17 21 5.1 662 9.7 sylvatica Planosol Fagus Ferralic Thetford 52.4100 0.8700 20 1930 31.0 3.78 19 3.9 528 9.7 sylvatica arenosol Fagus Sceletic Covert Wood 51.2000 1.1200 20 1950 31.0 5.15 20 7.5 527 10.5 sylvatica Leptosol Fagus Shobdon 52.2400 − 3.0300 200 1952 31.4 5.49 Cambisol 18 4.1 714 9.6 sylvatica Picea Tummel 56.7200 − 4.0500 400 1969 87.9 8.5 Ferric Podzol 37 5.1 796 8.1 sitchensis Picea Cabric Goyt 53.2900 − 1.9800 275 1981 59.2 n.a. 22 4.3 697 9.0 sitchensis Podzol Pinus Gleyic Rannoch 56.6500 − 0.4200 470 1965 46.6 6.13 38 4.2 1161 8.5 sylvestris Podzol Pinus Cabric Ladybower 53.4100 1.7500 275 1952 49.2 2.93 31 4.1 697 9.0 sylvestris Podzol Pinus Humic Rogate 51.0200 − 0.8700 80 1950 37.8 (#) n.a. 25 3.1 662 9.7 sylvestris Podzol Pinus Ferralic Thetford 52.4100 0.8700 20 1967 36.0 4.08 25 5.3 528 9.7 sylvestris Arenosol Quercus Eutric Alice Holt 51.1500 − 0.8600 80 1935 24.2 6.1 19 5.5 662 9.7 robur Vertisol Quercus Eutric Savernake 51.5900 − 1.9200 107 1950 23.1 n.a. 38 5.2 796 8.1 robur Vertisol Table 1.  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, ­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, ­201565. 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. increased ­Ndep has a fertilization effect on A, thereby contributing to enhancing a c­ a or climate driven increase in tree ­iWUE12,13 and forest C-sinks14,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 N ­ dep, 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 ­reached16. Stable nitrogen isotope composition in plant materials (δ15N) has been used as a proxy of changes in ecosystem N a­ vailability17. An increase in tree δ15N has been observed in studies along N ­ dep ­gradients18,19 and 20,21 15 soil N manipulation ­experiments . By contrast, a decrease in tree δ N is expected in the case of a reduction of soil N availability (e.g., under a reduction of N ­ dep22 or in case of tree canopy retention of N ­ dep23), 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 gs25. 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 ­site26 or species- specific22. 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 differ- ences in the main climate parameters and N ­ dep 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 gs) 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. Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 2 Vol:.(1234567890) www.nature.com/scientificreports/ 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 N ­ dep, 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 sig- nificantly increased by approx. 0.02–0.06 °C year−1. Whereas, vapour pressure deficit (VPD) did not show a sig- nificant 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). 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 gen- eral 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 represent- ing 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 c­ a increased by 51 ppm (from 339 to 390 ppm, slop e = 1.7 ± 3.9 ppm year−1)27. Wet ­Ndep and wet S­ dep decreased between 1995 and 2010 at all sites but Savernake and Tummel for S­ dep, and Goyt/Ladybower and Tummel for N ­ dep (Table 2). Despite the general decreasing Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 3 Vol.:(0123456789) www.nature.com/scientificreports/ Site Forest stand Parameter NH4 –N NO3-N Total ­Ndep Total ­Sdep SPEI Tgrs Tmaxgrs Pgrs VPDgrs Ta Name Species 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 Alice Holt Oak − 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)* 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) 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 ­NH4-N and N ­ O3-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. trends in wet S­ dep and N ­ dep, 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). 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‰ y­ ear−1, p < 0.01) (Fig. 2C). Changes in Δ18Ow were not cor- related 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 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 δ15Nw 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 ­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 Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 4 Vol:.(1234567890) www.nature.com/scientificreports/ A) Beech Oak Scots pine Sitka spruce B) Beech Oak Scots pine Sitka spruce 0.8 115 100 0.7 ) 1 iWUE(µmol mol 85 ci ca 0.6 70 Ndep level 0.5 55 High Low 40 0.4 Sites and species Alice Holt-beech 80 90 00 19 0 80 90 00 19 0 80 90 00 19 0 80 90 00 10 80 90 00 19 0 80 90 00 19 0 80 90 00 19 0 80 90 00 10 1 1 1 1 1 1 Alice Holt-oak 19 19 20 20 19 20 20 19 20 20 19 20 20 19 19 20 20 19 20 20 19 20 20 19 20 20 Year Year Covert Wood-beech Goyt -Sitka spruce C) Beech Oak Scots pine Sitka spruce D) Beech Oak Scots pine Sitka spruce Ladybower-Scots pine 36 4.0 Rannoch-Scots pine Rogate-Scots pine Savernake-oak 2.5 Shobdon-beech Thetford-beech 34 1.0 Thetford-Scots pine ) ) Tummel-Sitka spruce Ow ( Nw ( 0.5 18 15 32 2.0 3.5 30 5.0 80 90 00 19 0 80 90 00 19 0 80 90 00 19 0 80 90 00 10 80 90 00 19 0 80 90 00 19 0 80 90 00 19 0 80 90 00 10 1 1 1 1 1 1 19 19 20 20 19 20 20 19 20 20 19 20 20 19 19 20 20 19 20 20 19 20 20 19 20 20 Year Year 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 c­ i/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). (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 pro- cesses. 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 gradi- ent (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 predic- tors 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 ­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 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 criti- cal 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 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 Δ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 s­ Sdep and a­ Sdep. 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). Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 5 Vol.:(0123456789) www.nature.com/scientificreports/ A) B) 0.6 0.6 a T 0.4 T SPEI PC3 (19 %) a 0.4 Spatial PC2 (24 %) grs 12 a_max 0.2 grs_max 0.2 0.0 PP T SPEI, P a a 0.0 SPEI 8_3 grs a_max 0.2 Temporal Pa 0.2 Tgrs Pgrs p 0.4 VPDgrs Tgrs_max VPDa 0.6 VPD VPD a 0.4 VPD grs T, VPD 8_2 0.8 0.6 8_1 SPEI 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.6 0.4 0.2 0.0 0.2 0.4 0.6 Spatial PC1 (70 %) Temporal PC1 (30 %) C) D) 120 beech 120 oak Scots pine 110 Sitka spruce 110 )1 iWUE (µmol mol 100 100 90 90 80 80 70 70 60 60 6 3 0 3 6 6 3 0 3 6 Spatial PC1 Temporal PC1 Alice Holt beech Shobdon beech Ladybower Scots pine Rogate Scots pine Alice Holt oak Goyt Sitka spruce Sites and species Covert Wood beech Thetford beech Rannoch Scots pine Thetford Scots pine Savernake oak Tummel Sitka spruce 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. 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). Con- sistent 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 (nega- tive 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 δ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 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 s­ Sdep and ­aSdep on iWUE (Fig. S5). Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 6 Vol:.(1234567890) www.nature.com/scientificreports/ Auto-correlation Random effect Fixed effects Estimate ± SE Structure Standard deviation Residual variance R2m R2c (i) All sites—including ca, sNdep and sSdep (n = 251) Intercept 80.32 ± 2.30*** Conifer versus deciduous 8.26 ± 3.40* PCA_s1 − 2.65 ± 0.67* 0.76 3.82 7.15 0.38 0.52 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*** Conifer versus deciduous 7.73 ± 3.85 ¥ PCA_s1 − 2.19 ± 0.79* 0.50 5.21 4.93 0.41 0.72 PCA_a1 − 0.64 ± 0.14*** ca 0.20 ± 0.04*** (iii) All parameters (n = 140) Intercept 80.94 ± 2.33*** 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 ¥ 0.54 3.62 4.78 0.66 0.78 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*** Conifer versus deciduous 8.51 ± 2.66* PCA_s1 − 3.67 ± 0.63** PCA_a1 − 0.66 ± 0.17*** 0.47 2.85 4.40 0.71 0.80 PCA_a3 − 0.52 ± 0.25* sSdep − 0.47 ± 0.18* aSdep − 1.77 ± 0.40* Table 3.  Relationship between physiological and ecological parameters and environmental factors. 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. 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 ­Britain29. 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. Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 7 Vol.:(0123456789) www.nature.com/scientificreports/ Auto-correlation Random effect Fixed effects Estimate ± SE Structure Standard deviation Residual variance R2m R2c Δ Ow—(i) All sites including ca, sSdep, sNdep (n = 251) 18 Intercept 32.88 ± 0.18*** Conifer versus deciduous − 0.06 ± 0.27 n.s. PCA_s1 − 0.12 ± 0.05* PCA_a1 − 0.05 ± 0.013*** 0.25 0.38 0.41 0.31 0.63 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*** Conifer versus deciduous − 0.05 ± 0.01 n.s. PCA_s1 − 0.24 ± 0.02*** PCA_a1 − 0.06 ± 0.01*** 0.12 0.11 0.35 0.72 0.75 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*** Conifer versus deciduous − 0.16 ± 0.30 n.s. 0.27 0.45 0.44 0.10 0.56 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*** Conifer versus deciduous − 0.04 ± 0.37 n.s. 0.14 0.46 0.82 0.43 0.57 sNdep 0.14 ± 0.07 ¥ PCA_s1 0.51 ± 0.13** Table 4.  Relationship between physiological and ecological parameters and environmental factors. 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 R ­ 2c, 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. 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 S­ dep and N ­ dep ­(NH4-N + NO3-N) was also observed at most s­ ites1,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 ­decades4. 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 ­ammonia4. Taken all together, these results suggest that directionality and magnitude of changes in ­Ndep and S­ dep across the environmental gradient may modulate tree species response to c­ a, 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 ­analyses12,13,32. Across sites, the relationship between iWUE and Δ18Ow (Fig. 4) suggests that overall there is a tight coupling between A and gs 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. Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 8 Vol:.(1234567890) www.nature.com/scientificreports/ 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 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 Δ18Ow was included as fixed factor in the model for iWUE (Ref. Table S7 for all details). 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 gs, 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 l­ imitation33. 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 S­ dep has significantly decreased at the site, we cannot exclude a possible legacy (negative) effect of deposition on g­ s23. 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 S­ dep, which would have a posi- tive effect on ­gs10. 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 gs 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 gs, as suggested by no significant trend in Δ18Ow, while ­ci/ca decreased). This result indicated that ­Sdep is still negatively affecting oak trees at Sav- ernake—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 S­ O2 pollution significantly reduced ­gs leading to an increase in iWUE for oak trees in southern ­England10. The reduction in iWUE in the fast-growing34 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 ­availability35 and hydraulic ­conductance36, but also through increasing LAI as canopy develops and contributes to respired ­CO237. Note, however, that the trend in iWUE at Goyt cannot be explained by increasing hydraulic constraints with tree ­height36, 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 ­closure38. 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 Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 9 Vol.:(0123456789) 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 Δ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 g­ s, 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 explain- ing 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 ­US26, whereby changes in iWUE were predominantly associated with soil moisture conditions and no effects of either ­Ndep 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 ­ca, but also with the reduction in S­ dep, while there was no relationship with N ­ dep. This is quite important, as S­ dep and ­ca may have a similar effect on g­ s, i.e., an increase in S­ dep and c­ a leads to a reduction in g­ s10. 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 c­ a. The negative relationship between Δ18Ow and S­ dep suggests that the alleviation of the negative effects from S­ dep on g­ s22 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 ­Sdep outweigh changes in ­Ndep and ­ca. ­Sdep may influence tree water-use directly, by affecting ­gs22 or indirectly, by leading to soil acidification, loss of calcium, which play a significant role in controlling g­ s39. 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 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 δ15Nw 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 δ15Nw 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 ­NO3, leaving trees with 15 N-enriched ­NO340. This was only observed for the beech at Thetford, which is already confirmed as an N saturated site by detailed gradient ­studies41. 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 ­saturation41. 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 ­leaching41. 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 ­USA22,42 and ­globally43, 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 δ15Nw44, 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 ­analysis45, 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 N ­ dep-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 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 ­Ndep46. Our results also indicate that pinpointing variations in δ15Nw 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 ­Ndep and changes in the N isotopic signatures for ­NH4 and ­NO3 might increase the uncertainties in the detection of changes in atmospheric N ­ dep by using δ15N in tree r­ ings47. 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 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 10 Vol:.(1234567890) 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 ­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 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 δ15Nw were mostly explained by dif- ferences 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 ­drought48, 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 ­pollution49. Sitka spruce, the major upland planted timber species in the UK, native of coastal Northwest- ern 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 Brit- ain, 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, Shob- don and Goyt. The Goyt site was added as a high ­Ndep site as a contrast to the low ­Ndep Sitka spruce site in Scotland (Fig. 1, Table 1, Supplementary Table 1). For each species, forests were selected with similar soil type and age, but with contrasting ­Ndep, ­Sdep and climate, particularly rainfall and temperature, as described in Fig. 1, Table 1 and Supplementary Table 1. Stand information (mean tree height, mean diameter at the breast height and basal area) as measured for target years and for some of the forest stands are shown in Fig. S4. At each ICP forest site, a plot of 0.25 ha was established in 1995 to carry out ­monitoring30 and a similar pro- tocol 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 mete- orological 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 grow- ing 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 condi- tions. 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 previ- ously ­described30. For the spatial analyses, we considered mean of annual deposition ­(sNdep and ­sSdep), obtained as the sum of total ­(NH4-N + NO3-N for ­Ndep) wet and dry deposition. The latter were estimated as difference between throughfall and bulk precipitation N ­fluxes30. For Rogate only 1 year (2010) of monitoring was available. For Goyt site, atmospheric deposition data collected at Ladybower were considered, as the two sites are 30 km apart. For two sites (i.e., Shobdon and Covert Wood), which were not part of the regular ICP forest sites, the wet deposition obtained from the UK 5 × 5 km grid N ­ sed4. The estimate included wet ­ dep and S­ dep dataset was u and dry ­NHx-N ­(NH4, ­NH3), ­NOy-N ­(NO2, ­NO3, ­HNO3) and S (­ SOx = SO2 and S­ O4) deposition, modelled using FRAME with 2005 emissions d ­ ata4. 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 ­ dep50. For the temporal analyses, only wet deposition (as calculated from N ­ O3 , ­NH4+ and S­ O2–4 concentrations − in bulk precipitation) was considered (indicated as a­ Ndep 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 Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 11 Vol.:(0123456789) 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 δ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 N ­ 2 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 c­ i/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:  13 δ Ca − δ 13 Cw  13 ci � Cw = a + (b − a) =  13 C  (1) ca 1 + δ1000 w 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: (δ 13 Ca − δ 13 Cw ) − a ci = ca (2) b−a ca values were obtained from Mauna Loa r­ ecords27, and δ13Ca values were obtained from Mauna Loa r­ ecords5 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: ca b − 13 Cw   A ca − ci iWUE = = = (3) gs 1.6 1.6 b−a where 1.6 is the molar diffusivity ratio of C ­ O2 to H ­ 2O (i.e., g­ CO2 = 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: δ 18 Ow − δ 18 Os �18 Ow =  18  (4) Os 1 + δ1000 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: Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 12 Vol:.(1234567890) www.nature.com/scientificreports/ √ δ 18 OP = 0.52Ta − 0.006Ta2 + 2.42Pa − 1.43Pa2 − 0.046 E − 13.0 (5) 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​://water​iso.utah.edu/water​isoto​pes/pages​/data_acces​s/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 ­gs56. The physiological signal imprinted in the foliage may be dampened in tree rings, due to post-photosynthetic fractionation during translocation of sucrose and synthesis of cellulose in the tree ­stem57. This leads to an offset between foliar and tree-ring δ18O and also δ13C values. However, accounting for the offset when interpreting tree-ring isotopes is still challenging, as it is not clear whether the offset is species-specific, if it is maintained over the long-term and what are the mechanisms driving i­ t57,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 relation- ship 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 tem- perature, 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., a­ Sdep and a­ Ndep, to assess the effect of temporal changes) and the mean over the monitoring time (i.e., ­sSdep and s­ Ndep) to evaluate both the contri- bution of within-site temporal changes and cross-sites differences on changes in iWUE, ­ci/ca, Δ18Ow and δ15Nw. Note that we only have data for the aSdep and aNdep at 10 of the 12 sites. To eliminate auto-correlation between mean site variables and yearly variables for each site, site variables were globally centered, whereas yearly data were group-centered59. To eliminate auto-correlation among individual climatic variables of the two groups, we conducted two principal components (PC) analyses, after centering and scaling the variables. The first PC analysis considered across sites long-term averages of environmental variables (PCA_s), the second within sites annual time series (PCA_a). Tables of percentage of variance explained and scree plots were examined to deter- mine 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 anthro- pogenic (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 correla- tion 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 measure- ment lags. Marginal (only fixed factors) and conditional (fixed + random factors) percent of explained variance ­(R2m and ­R2c, respectively) were calculated using package MuMIn62. Finally, to examine the joint effects of climatic conditions on Δ18Ow, δ15Nw and iWUE, a mixed-effect con- firmatory path model was employed using package piecewiseSEM63. For each isotope and iWUE, the final model from linear mixed effect model analyses (Tables 3, 4, Supplementary Table S7) was considered, with the only modification of excluding PFT as fixed factor and including also Δ18Ow as fixed factor in the model for iWUE and the correlated errors between PC_s1 and sNdep. PiecewiseSEM allows the fitting of hierarchical models with random effects on data with a multivariate structure, allowing for the identification of indirect effects and unresolved covariance among endogenous variables. Model goodness-of-fit was assessed using a chi-square test against Fisher C, based on Shipley’s test of direct separation, which tests for the overall significance of missing Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 13 Vol.:(0123456789) www.nature.com/scientificreports/ relationships among unconnected variables, while the significance of any given missing path was evaluated using individual p-values. A combination of non-significant Fisher-C and individual p-values tests implies that the hypothesised relationships among the variables are consistent with the data without missing any significant ­relationship63. All statistical analyses were carried out inside RStudio 1.0.143 with R 3.4.064. Data availability All isotope data used in this paper are available on Zenodo: https​://doi.org/10.5281/zenod​o.39078​49. Received: 25 July 2018; Accepted: 8 June 2020 References 1. Lajtha, K. & Jones, J. Trends in cation, nitrogen, sulfate and hydrogen ion concentrations in precipitation in the United States and Europe from 1978 to 2010: a new look at an old problem. Biogeochemistry 116, 303–334 (2013). 2. Galloway, J. N. et al. Nitrogen cycles: past, present, and future. Biogeochemistry 70, 153–226 (2004). 3. Jia, Y. et al. Global inorganic nitrogen dry deposition inferred from ground- and space-based measurements. Sci. Rep. 6, 19810 (2016). 4. RoTAP. Review of Transboundary Air Pollution: Acidification, Eutrophication, Ground Level Ozone and Heavy Metals in the UK. Contract Report to the Department for Environment, Food and Rural Affairs. Centre for Ecology & Hydrology (2012). 5. Scripps 2018. Scripps ­CO2 programme—Keeling, C. D., et al.. Atmospheric ­CO2 and 13CO2 exchange with the terrestrial biosphere and oceans from 1978 to 2000: observations and carbon cycle implications, pages 83–113, in A History of Atmospheric CO2 and its effects on Plants, Animals, and Ecosystems (editors, Ehleringer, J.R., T. E. Cerling, M. D. Dearing). Springer, New York, 2005. https​ ://scrip​psco2​.ucsd.edu/data/atmos​pheri​c_co2/mlo (2018). 6. Ainsworth, E. A. & Rogers, A. The response of photosynthesis and stomatal conductance to rising ­[CO2]: mechanisms and envi- ronmental interactions. Plant Cell Environ. 30, 258–270 (2007). 7. Franks, D. C. et al. Water-use efficiency and transpiration across European forests during the anthropocene. Nat. Clim. Chang. 5, 579–584 (2015). 8. Ehleringer, J. R., Hall, A. E. & Farquhar, G. D. Water use in relation to productivity. In Stable Isotopes and Plant Carbon-Water Relations (eds Ehleringer, J. R. et al.) 3–8 (Academic Press, New York, 1993). 9. Saurer, M., Siegwolf, R. T. W. & Schweingruber, F. H. Carbon isotope discrimination indicates improving water-use efficiency of trees in northern Eurasia over the last 100 years. Glob. Change Biol. 10, 2109–2120 (2004). 10. Rinne, K., Loader, N., Switsur, V., Treydte, K. & Waterhouse, J. Investigating the influence of sulphur dioxide (­ SO2) on the stable isotope ratios (δ13C and δ18O) of tree rings. Geochim. Cosmochim. Acta 74, 2327–2339 (2010). 11. Thomas, R. Q., Canham, C. D., Weathers, K. C. & Goodale, C. L. Increased tree carbon storage in response to nitrogen deposition in the US. Nat. Geosci. 3, 13–17 (2009). 12. Peñuelas, J., Canadell, J. G. & Ogaya, R. Increased water-use efficiency during the 20th century did not translate into enhanced tree growth. Glob. Ecol. Biogeogr. 20, 597–608 (2010). 13. Leonardi, S. et al. Assessing the effects of nitrogen deposition and climate on carbon isotope discrimination and intrinsic water-use efficiency of angiosperm and conifer trees under rising ­CO2 conditions. Glob. Change Biol. 18, 2925–2944 (2012). 14. Magnani, F. et al. The human footprint in the carbon cycle of temperate and boreal forests. Nature 447, 849–851 (2007). 15. Fernández-Martínez, M. et al. Atmospheric deposition, C ­ O2, and change in the land carbon sink. Sci. Rep. 7, 1–13 (2017). 16. Aber, J. D. et al. Is nitrogen deposition altering the nitrogen status of northeastern forests?. Bioscience 53, 375 (2003). 17. Craine, J. M. et al. Ecological interpretations of nitrogen isotope ratios of terrestrial plants and soils. Plant Soil 396, 1–26 (2015). 18. Pardo, L. H. et al. Regional assessment of N saturation using foliar and root δ15N. Biogeochemistry 80, 143–171 (2006). 19. Emmett, B. et al. Natural abundance of δ15N in forests across a nitrogen deposition gradient. For. Ecol. Manag. 101, 9–18 (1998). 20. Elhani, S., Guehl, J.-M., Nys, C., Picard, J.-F. & Dupouey, J.-L. Impact of fertilization on tree-ring δ15N and δ13C in beech stands: a retrospective analysis. Tree Physiol. 25, 1437–1446 (2005). 21. Balster, N. J., Marshall, J. D. & Clayton, M. Coupling tree-ring δ13C and δ15N to test the effect of fertilization on mature Douglas-fir (Pseudotsuga menziesii var. glauca) stands across the Interior northwest, USA. Tree Physiol. 29, 1491–1501 (2009). 22. Mathias, J. M. & Thomas, R. B. Disentangling the effects of acidic air pollution, atmospheric C ­ O2, and climate change on recent growth of red spruce trees in the Central Appalachian Mountains. Glob. Change Biol. https​://doi.org/10.1111/gcb.14273​ (2018). 23. Guerrieri, R. et al. The legacy of enhanced N and S deposition as revealed by the combined analysis of δ13C, δ18O and δ15N in tree rings. Glob. Change Biol. 17, 1946–1962 (2010). 24. Craine, J. M. et al. Isotopic evidence for oligotrophication of terrestrial ecosystems. Nat. Ecol. Evol. 2, 1735–1744 (2018). 25. Farquhar, G. D., Cernusak, L. A. & Barnes, B. Heavy water fractionation during transpiration. Plant Physiol. 143, 11–18 (2007). 26. Levesque, M., Andreu-Hayles, L. & Pederson, N. Water availability drives gas exchange and growth of trees in northeastern US, not elevated ­CO2 and reduced acid deposition. Sci. Rep. 7, 46158 (2017). 27. NOAA ESRL Global Monitoring Division. 2016, updated annually. Atmospheric Carbon Dioxide Dry Air Mole Fractions from quasi-continuous measurements at Mauna Loa, Hawaii. Compiled by K.W. Thoning, D.R. Kitzis, and A. Crotwell. National Oce- anic and Atmospheric Administration (NOAA), Earth System Research Laboratory (ESRL), Global Monitoring Division (GMD): Boulder, Colorado, USA. Version 2017–2018 at https​://doi.org/10.7289/V54X5​5RG (2016) 28. Hall, J., Emmett, B., Garbutt, A., Jones, L., Rowe, E., Sheppard, L., Vanguelova, E., Pitman, R., Britton, A., Hester, A., Ashmore, M., Power, S. and Caporn, S. UK Status Report July 2011: Update to empirical critical loads of nitrogen. Report to Defra under contract AQ801 Critical Loads and Dynamic Modelling. NERC/Centre for Ecology & Hydrology, 57 p. (2011). 29. Murphy, J. M. et al. UK Climate Projections Science Report: Climate change projections (Met Office Hadley Centre, Exeter, 2009). 30. Vanguelova, E. et al. Chemical fluxes in time through forest ecosystems in the UK—soil response to pollution recovery. Environ. Pollut. 158, 1857–1869 (2010). 31. Waldner, P. et al. Detection of temporal trends in atmospheric deposition of inorganic nitrogen and sulphate to forests in Europe. Atmos. Environ. 95, 363–374 (2014). 32. Saurer, M. et al. Spatial variability and temporal trends in water-use efficiency of European forests. Glob. Change Biol. 20, 3700–3712 (2014). 33. Irvine, J., Perks, M. P., Magnani, F. & Grace, J. The response of Pinus sylvestris to drought: stomatal control of transpiration and hydraulic conductance. Tree Physiol. 18, 393–402 (1998). 34. Cameron, A. D. Building Resilience into Sitka Spruce (Picea sitchensis (Bong.) Carr.) Forests in Scotland in Response to the Threat of Climate Change. Forests 6, 398–415 (2015). 35. Brienen, R. J. W. et al. Tree height strongly affects estimates of water-use efficiency responses to climate and ­CO2 using isotopes. Nat. Commun. https​://doi.org/10.1038/s4146​7-017-00225​-z (2017). Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 14 Vol:.(1234567890) www.nature.com/scientificreports/ 36. McDowell, N. G., Bond, B. J., Dickman, L. T., Ryan, M. G. & Whitehead, D. Relationships between tree height and carbon isotope discrimination. In Size- and Age-Related Changes in Tree Structure and Function (eds Meinzer, F. C. et al.) 255–285 (Springer, New York, NY, 2011). 37. Buchmann, N., Kao, W.-Y., & Ehleringer, J. Influence of stand structure on carbon-13 of vegetation, soils, and canopy air within deciduous and evergreen forests in Utah, United States. Oecologia 110, 109–11 (1997). 38. Manrique-Alba, A. et al. Long-term thinning effects on tree growth, drought response and water use efficiency at two Aleppo pine plantations in Spain. Sci. Total Environ. 728, 138536 (2020). 39. Lanning, M. et al. Intensified vegetation water use under acid deposition. Sci. Adv. 5, 7. https​://doi.org/10.1126/sciad​v.aav51​6 (2019). 40. Högberg, P. 15N natural abundance in soil–plant systems. New Phytol. 137, 179–203 (1997). 41. Vanguelova, E. I. & Pitman, R. Nutrient and carbon cycling along nitrogen deposition gradients in broadleaf and conifer forest stands in the east of England. For. Ecol. Manag. 447, 180–194 (2019). 42. Mclauchlan, K. K. et al. Centennial-scale reductions in nitrogen availability in temperate forests of the United States. Sci. Rep. 7, 1–7 (2017). 43. Craine, J. et al. Isotopic evidence for oligotrophication of terrestrial ecosystems. Nat. Ecol. Evol. 2, 1735–1744 (2018). 44. Tomlinson, G. et al. Can tree-ring δ15N be used as a proxy for foliar δ15N in European beech and Norway spruce?. Trees 30, 627–638 (2016). 45. Craine, J. M. et al. Global patterns of foliar nitrogen isotopes and their relationships with climate, mycorrhizal fungi, foliar nutrient concentrations, and nitrogen availability. New Phytol. 183, 980–992 (2009). 46. Dore, A. J., Choularton, T. W. & Fowler, D. An improved wet deposition map of the United Kingdom incorporating the seeder- feeder effect over mountainous terrain. Atmos. Environ. Part A Gener. Top. 26(8), 1375–1381 (1992). 47. Gerhart, L. M. & Mclauchlan, K. K. Reconstructing terrestrial nutrient cycling using stable nitrogen isotopes in wood. Biogeo- chemistry 120, 1–21 (2014). 48. Cavin, L. & Jump, A. S. Highest drought sensitivity and lowest resistance to growth suppression are found in the range core of the tree Fagus sylvatica L. not the equatorial range edge. Glob. Change Biol. 23, 362–379 (2017). 49. Brown, N., Vanguelova, E., Parnell, S., Broadmeadow, S. & Denman, S. Predisposition of forests to biotic disturbance: predicting the distribution of acute oak decline using environmental factors. For. Ecol. Manag. 407, 145–154 (2018). 50. Guerrieri, R., Vanguelova, E. I., Michalski, G., Heaton, T. H. E. & Mencuccini, M. Isotopic evidence for the occurrence of biological nitrification and nitrogen deposition processing in forest canopies. Glob. Change Biol. 21, 4613–4626 (2015). 51. Farquhar, G. D., Ehleringer, J. R. & Hubick, K. T. Carbon isotope discrimination and photosynthesis. Annu. Rev. Plant Physiol. Plant Mol. Biol. 40, 503–537 (1989). 52. McCarroll, D. & Loader, N. J. Stable isotopes in tree rings. Quat. Sci. Rev. 23, 771–801 (2004). 53. Barbour, M. M. Stable oxygen isotope composition of plant tissue: a review. Funct. Plant Biol. 34, 83–94 (2007). 54. Barbour, M. M., Andrews, T. J. & Farquhar, G. D. Correlations between oxygen isotope ratios of wood constituents of Quercus and Pinus samples from around the world. Funct. Plant Biol. 28(5), 335–348. https​://doi.org/10.1071/PP000​83 (2001). 55. Barbour, M. M., Roden, J. S., Farquhar, G. D. & Ehleringer, J. R. Expressing leaf water and cellulose oxygen isotope ratios as enrich- ment above source water reveals evidence of a Péclet effect. Oecologia 138(3), 426–435 (2004). 56. Barbour, M. M., Fischer, R. A., Sayre, K. D. & Farquhar, G. D. Oxygen isotope ratio of leaf and grain material correlates with stomatal conductance and yield in irrigated, field-grown wheat. Aust. J. Plant Physiol. 27, 625–637 (2000). 57. Gessler, A. et al. Stable isotopes in tree rings: towards a mechanistic understanding of isotope fractionation and mixing processes from the leaves to the wood. Tree Physiol. 34, 796–818 (2014). 58. Cernusak, L. A. et al. Why are non-photosynthetic tissues generally 13C enriched compared with leaves in C3 plants? Review and synthesis of current hypotheses. Funct. Plant Biol. 36, 199–213 (2009). 59. Gelman, A. & Hill, J. Data Analysis Using Regression and Multilevel/Hierarchical Models (Cambridge University Press, Cambridge, 2007). 60. Pinheiro, J. C. & Bates, D. M. Mixed-Effects Models in S and S-PLUS (Springer, New York, 2000). 61. Hothorn, T., Bretz, F. & Westfall, P. Simultaneous inference in general parametric models. Biomet. J. 50(3), 346–363 (2008). 62. Bartoń, K. MuMIn: multi-model inference. Model selection and model averaging based on information criteria (AICc and alike) (2016). https​://CRAN.R-proje​ct.org/packa​ge=MuMIn​. 63. Lefcheck, J. S. PiecewiseSEM: piecewise structural equation modeling in R for ecology, evolution, and systematics. Methods Ecol. Evol. 7(5), 573–579. https​://doi.org/10.1111/2041-210X.12512​ (2015). 64. R Studio Team.RStudio: Integrated Development for R. RStudio, Inc., Boston, MA (2016). https​://www.rstud​io.com/. 65. IUSS Working Group WRB. World Reference Base for Soil Resources 2014, Update 2015 International Soil Classification System for Naming Soils and Creating Legends for Soil Maps. World Soil Resources Reports No. 106. FAO, Rome (2015). Acknowledgements R.G. thanks the follow-up support (2015–2017) from the Royal Society within the Newton International Fellow- ship and the MSCA individual fellowship (NITRIPHYLL, Project No. 705432), funded within the EU-HORIZON 2020. The sites are part of the European long term Intensive Forest Monitoring Network (ICP Level II) and the long-term research of the Soil Sustainability programmes at Forest Research. We would like to thank the Forest Research and the Forestry Commission for funding this research. The staff from the Technical Support Units and the Chemical Laboratory at Forest Research are also acknowledged for collecting and analysing the soil, plant and water samples for over 15 years at the Level II network in Britain. We thank anonymous reviewers for constructive criticisms and suggestions on earlier versions of the manuscript. Author contributions R.G., E.V., M.M., R.P. planned the study, with input from M.P., J.M. and S.B.. R.G. carried out the sampling of wood cores across the investigated sites with the support of E.V. and R.P. and with technical help from Forest Research. E.V., R.P. and S.B. were responsible for collecting soil site information and atmospheric deposition at the investigated sites. R.G. carried out sample preparation for isotope analyses and was responsible for data analyses with the support of M.M.. R.G. led paper writing, with contribution of all the co-authors. Competing interests The authors declare no competing interests. Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 15 Vol.:(0123456789) www.nature.com/scientificreports/ Additional information Supplementary information is available for this paper at https​://doi.org/10.1038/s4159​8-020-67562​-w. Correspondence and requests for materials should be addressed to R.G. Reprints and permissions information is available at www.nature.com/reprints. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creat​iveco​mmons​.org/licen​ses/by/4.0/. © The Author(s) 2020 Scientific Reports | (2020) 10:12418 | https://doi.org/10.1038/s41598-020-67562-w 16 Vol:.(1234567890)