Introduction

In the last six decades, intensification in agricultural practices and fossil fuel combustion have strongly increased both atmospheric nitrogen (N) emission and its deposition to the terrestrial environment on a global scale (Erisman et al. 2011; Wessel et al. 2013). At a European scale, increased N inputs were generally followed by subsequent decreasing trends in N depositions between 1990 and 2010 (Theobald et al. 2019). Increased anthropogenic N inputs can alter the N cycle and affect ecosystem functions. The impact of increased N input to the environment depends on the fate of N in a specific ecosystem (Niu et al. 2016). In N-limited forest ecosystems, increased N deposition might contribute to extra N uptake by trees or accumulation of N in the soil (Nadelhoffer et al. 1999; Kjønaas and Stuanes 2008; Templer et al. 2012). However, long-term chronic N deposition may also exceed the capacity for N uptake by vegetation and microbes and cause N saturation and N loss from the soil (Aber et al. 1989; Lovett and Goodale 2011; Lu et al. 2011). Aber et al. (1989) defined N saturation as the availability of ammonium (NH4+) and nitrate (NO3−) in excess of total combined plant and microbial nitrogen demand. N saturation would be manifested as NH4+ accumulation and/or NO3− leaching (Gundersen and Rasmussen 1995). Nitrogen saturation may lead to soil acidifications and forest decline and increased N leaching losses can cause eutrophication of freshwaters (Aber et al. 1989; Jonard et al. 2015). Furthermore, increased N deposition can contribute to the production of the greenhouse gas N2O.

To investigate the impact of N deposition on forest systems, N deposition manipulation experiments were set up in Europe and North America in the late 1980s and early 1990s (Dise and Wright 1995; Nadelhoffer et al. 1999). A number of these N manipulation experiments included a one-time addition of 15N-enriched NH4+ and/or NO3− on the forest floor. The use of 15N enriched deposition enables monitoring of incoming 15N in various vegetation and soil N pools. Increased N deposition could have an impact on the fate of 15N in these N pools or leaching or gaseous N losses. When increased N deposition would result in increased 15N uptake by vegetation this could imply growth enhancement and increased C storage in tree biomass. If, on the other hand, the majority of the added 15N would be immobilized in the soil, this may generally suggest a lower potential for increased C stocks, considering the lower C/N ratios of soils relative to tree biomass in these forest ecosystems (Nadelhoffer et al. 1999). Furthermore, increased N inputs have been found to cause lower recovery levels of 15N in the soil as a result of 15N leaching (Tietema et al. 1998; Nadelhoffer et al. 1999; Templer et al. 2012). These studies reflect 15N recovery after a relatively short period (3–18 months). Longer time scales are needed to investigate long-term turnover and stabilization processes of 15N in forest systems.

A few studies have investigated 15N recovery in temperate forest ecosystems over longer time scales (Nadelhoffer et al. 2004; Krause et al. 2012; Wessel et al. 2013; Goodale 2017). In most studies, the organic soil layer remained the dominant 15N sink 6–9 years after 15N addition (Nadelhoffer et al. 2004; Krause et al. 2012; Wessel et al. 2013), whereas at one of these sites, the mineral soil became the dominant 15N sink already after 1 year (Goodale 2017). From a resampling of the Ysselsteyn forest in the Netherlands (Wessel et al. 2013) 19 years after 15N addition, it was shown that 15N shifted from the organic to the mineral layer between 8 and 19 years after labelling (Wessel et al. unpublished). The mineral soil might have a larger potential than the organic soil to store organic N on the long-term, due to stabilization of organic N in organo-mineral associations (Moni et al. 2012). Based on the results of Wessel et al. (2013 and unpublished), the 15N recovery in the mineral soil layers increases with time. Yet, there are several other mechanisms contributing to 15N recovery in both the organic and mineral soil layers. N that is stored in the soil is strongly linked to C dynamics (Kopáček et al. 2013; Cheng et al. 2018), which reflects the sensitive interaction between microbes and the internal NO3− cycling (Durka et al. 1994; Tahovská et al. 2020). Most drivers are therefore assumed to be related to C input and C stabilization processes (Wiesmeier et al. 2019). Furthermore, 15N recovery in the soil seems to depend on the N-status of the forest ecosystem and the capability of this system to immobilize extra N inputs (Gundersen et al. 1998).

Due to the limited number of studies on the long-term fate of deposited N in temperate forest systems, we re-investigated several 15N labelling and NITRogen saturation EXperiments (NITREX) experiments approximately 20 years after labelling (Wright and van Breemen 1995). We focused on 15N recovery in the organic and mineral soil layers, as these soil N pools were the dominant 15N sinks in previous studies (Nadelhoffer et al. 1999). Our first hypothesis was that the 15N recovery in the mineral soil layers would become larger than in the organic soil layer ~ 2 decades after labelling, due to a translocation of the 15N tracer in the soil horizons over time as organic matter infiltrated deeper in the mineral soil horizons. Second, based on the short-term results of Templer et al. (2012) we hypothesized that long-term increased N deposition would lead to N-saturation and N-leaching and thus lower 15N recovery levels in the soil. The sites varied in greater or lesser extent from each other related to several soil and site characteristics. We investigated the impact of these characteristics on 15N recovery. Third, we hypothesized that variables such as, C stock, C content and layer thickness, would have a large impact on 15N recovery, due to the strong link with C dynamics. Furthermore, variables such as N-status, NH4+ concentration and C:N ratio would be negatively correlated with 15N recovery, due to the effects of N-saturation and N-leaching on 15 N recovery in the soil. Several other variables were included as described in the methodology below.

Materials and methods

Site descriptions

The study was carried out at four European coniferous forest sites: Gårdsjön, Sweden (GD), Klosterhede, Denmark (KH), Speuld, The Netherlands (SP) and Alptal, Switzerland (AL). We used one additional site from Ysselsteyn, The Netherlands (YS) for some of our analysis. This site had the same experimental design as our four sites and was re-sampled in 2001 (Wessel et al. 2013) and in 2012 (Wessel et al. unpublished). All of the sites were part of the nitrogen saturation experiment (NITREX) that started in the 1980s (Dise and Wright 1995). The scales of the experiments varied from plot scale (KH, SP and YS) to catchment scale (AL and GD). More details about the scales and replication of the experiments can be found in the site-specific papers (Table 1). The dominant tree species were Norway spruce (Picea abies L.) at GD, KH and AL, Douglas fir (Pseudotsuga menziesii [Mirb.] Franco) at SP and Scots pine (Pinus sylvestris L.) at YS. The soils at GD, KH and SP were classified as Podzols and at AL as Gleysols (Gundersen et al. 1998; Schleppi et al. 1998). The AL catchment is characterized by a topography of depressions and mounds that results in a variation in vegetation and humus types (Providoli et al. 2005). The Umbric Gleysols (IUSS Working Group WRB 2014) at AL has a low permeability with a water table close to the surface throughout the year (Hagedorn et al. 2001). SP and YS are well-drained, flat topography sites with Haplic Podzols (IUSS Working Group WRB 2014) with a mor-moder humus type (Koopmans et al. 1996). KH is a well-drained, flat topography site with a Haplic Podzol (IUSS Working Group WRB 2014) with a mor humus type (Gundersen and Rasmussen 1995). The GD catchment is characterized by a topography of depressions and granitic outcrops, that results in a variation of Haplic Podzols and Gleyic Podzols with typic Leptic Histosols (IUSS Working Group WRB 2014) at the shallow outcrops (Kjønaas et al. 1998).

Table 1 Ambient N deposition, treatment N flux, type of 15N addition, first year of 15N addition, mean annual temperature, mean annual precipitation, forest stand age and scales of the experiment at the start of the experiment (in 1990) for all of the investigated sites (GD, AL, KH, SP) and including Ysselsteyn (YS), which was resampled by Wessel et al. (2013) and Wessel et al. (unpublished) and will be partly included in this paper
Full size table

These five sites covered a wide range (13–58 kg N ha−1 year−1) of ambient N deposition in West-Europe at that time. To investigate effects of N deposition on soil N saturation, N deposition was either decreased or increased, depending on the ambient N deposition (Table 1). These N treatments were compared with a control plot that received 15N input under ambient N deposition, except for GD where the control catchment was not labelled with 15N. At the sites with a relatively low ambient N deposition (< 20 kg N ha−1 year−1; AL, KH and GD), N deposition was increased to a total of approximately 40–50 kg N ha−1 year−1. At SP and YS which already had a high ambient N deposition (50 kg N ha−1 year−1), N input was reduced to ~ 5 kg N ha−1 year−1. At AL KH and GD, the increased levels of N deposition were sprayed with a sprinkler system or backpack sprayer on the forest floor as ammonium nitrate (NH4NO3) (Gundersen et al. 1998; Moldan and Wright 1998; Schleppi et al. 1998). This was done continuously from the start of the experiment until present. At SP and YS, ambient N was intercepted by a roof that was placed beneath the canopy. An artificial throughfall, with the same relative content of NO3− and NH4+ as ambient N, but with a decreased level of N (5 kg N ha−1 year−1) was sprayed beneath the roof with a sprinkler system (Koopmans et al. 1996). The roof was built in 1989 and removed in 1995 (SP), so from that moment onwards the manipulated plots received ambient N deposition. To determine the fate of the deposited N in the forest ecosystems at the various sites, both low and high N depositions were enriched with 15NH4 and/or 15NO3 for a period of 12 months in 1992/1993 at SP, KH, GD and in 1995/1996 at AL (high N catchment) and in 2000/2001 (K15NO3) and 2002/2003 (15NH4Cl) at the AL control catchment (Table1).

Re-sampling and laboratory analysis

More than 20 years after 15N labelling, the sites were re-sampled in March 2015 (SP), in June 2015 (GD and KH) and in May 2016 (AL). The soils were re-sampled based on the original sampling schemes of the individual sites to have the possibility to compare the newly obtained data with the old data. At GD, which had a high spatial variability, the separation of the soil into horizons was done by the same person as for the original sampling. At the two plot-scale sites (SP and KH) 5 soil cores were randomly taken in the low and high N deposition plots. The organic soil layers (LF1 and F2) and the 0–10 cm and 10–25 cm of the mineral soil were sampled with a plastic cylinder. For the organic soil layers a cylinder with a 15.1 cm diameter was used and for the mineral soil a cylinder with 4.3 cm was used. The deeper soil layer (25–50 cm) was sampled with a gauge of 1.2 cm in diameter. At the two catchment-scale sites (AL and GD), soil cores were taken in a grid and separated based on diagnostic horizons differentiation. The catchments at AL were roughly 1500 m2 and at GD the catchment had a size of 5250 m2. At AL 20 cores per catchment were taken with a soil corer with a 5 cm diameter. The cores were 25 cm long and were divided into the organic soil layer, A-horizon and B-horizon (Schleppi et al. 1999). At GD, 51 cores were taken by use of a cylindrical auger with a diameter of 7 cm for the topsoil and an Edelmann auger for the deeper horizons. The soil was sampled down to 1 m soil depth or to bedrock at soil depth < 1 m. The soil samples were divided into diagnostic horizons, following the same sampling scheme of 20 years earlier (Kjønaas et al. 1998). To be able to compare all plots with catchment scale sites, the horizon differentiations at AL and GD were converted to the soil depth differentiations, where the measured soil parameters in the mineral soil horizons of various thicknesses were recalculated based on the following soil depths: 0–10 cm, 10–25 cm and 25–50. In contrast to the depths of the mineral soil layers, the depths of the organic soil layers, were not fixed.

Whereas the soil samples that were sampled in the previous NITREX study was analysed in various labs (Emmett et al. 1998), all soil samples in the current study were analysed at the University of Amsterdam. Small aliquots (5 to 10 g) of the field moist soil samples were dried at 70 °C (organic samples) or 105 °C (mineral samples) to measure moisture content and calculate dry mass and bulk-density (BD). The remaining samples were dried at 40 °C and sieved (2 mm). Small aliquots (5 to 30 g) of sieved soil were milled with a ball mill (mineral samples) or a centrifuge mill (organic samples) to a fine powdery and homogenous material. Small tin capsules were filled with 10 mg (organic) or 50 mg (mineral) soil to measure N and C content with an Elemental Analyser (EA, Vario ISOTOPE cube, Hannau, Germany), which was coupled to a continuous-flow Isotope Ratio Mass Spectrometer (IRMS, Vision Isoprime Manchester, United Kingdom), to measure 15N abundance of the samples. Reference gas (high-purity N2 gas) was calibrated to atmospheric N2 standard (at-air) using certified reference materials (IAEA-N2 and IAEA-NO3 from the Department of Nuclear Applications, International Atomic Energy Agency, Vienna).

Additionally, the sieved soil (< 2 mm) was used for soil chemical analysis to characterize the belowground environment at these sites. Electrical conductivity (EC) and pH were measured with a soil to water ratio that varied between 1:2.5, 1:5, and 1:10. This ratio increased with organic matter content. Cation concentration (K+, Na+, Mg2+, Ca2+, Fe3+, Al3+), effective cation exchange capacity (ECEC) and NH4+ concentrations in the mineral soil samples were measured in barium chloride (BaCl2) extracts of 0.1 M BaCl2 with a soil to BaCl2 ratio of 1:25 (Jaremko and Kalembasa 2014). The cation concentrations were measured on a Perkin Elmer Optima 8000 Optical Emission Spectrometer coupled with a Perkin Elmer S10 auto sampler (ICP-OES, Singapore) whereas NH4+ and NO3− concentrations were measured calorimetrically with a segmented flow SAN++ auto-analyser manufactured by Skalar (Breda, The Netherlands).

Calculations

The N calculated based on total N mass in the dried sample, divided by the area of the cylinder that was used to take the sample. Since we used dried material, and not sieved material for the calculation of the N pool, the coarse fractions were also taken into account. At most sites this cylinder had a fixed size. Therefore, BD was calculated based on the total dry mass divided by the total volume of the cylinder for each soil layer, which may have overestimated the N pools of the soil (Walter et al. 2016). As the lower mineral soil at GD was sampled with an Edelmann auger, Kjønaas et al. (1998) measured the BD in 4 soil profiles situated in representative areas of the dominant vegetation types of the catchment. The diagnostic horizons of those 4 soil profiles were sampled by use of 100 cm3 sample rings. BD was calculated based on sieved soil (< 2 mm) divided by the total volume of the sample ring. The BD dataset of Kjønaas et al. (1998) was used to recalculate the diameters of the cores, belonging to the depths and dry masses that we measured in the field. The N (and C) pools were calculated based on mass per area, which introduced uncertainties related to sampling of deeper mineral soil layers with an Edelmann auger at GD. This, along with uncertainties introduced through recalculating thickness of soil layers to given soil depths, are included as part of the discussion of the results.

The abundance of 15 N in the soil N pools was calculated with the following equation:

$$\delta^{15} {\text{N}} = \left( {\frac{{{\text{ R}}_{{{\text{sample}}}} }}{{{\text{ R}}_{{{\text{standard}}}} }} - 1} \right) \times 100\permille$$
(1)

where δ15N was expressed as permille delta 15N (δ15N ‰) of total N. R is the molar ratio of 15N/14N, with atmospheric N2 used as standard (Rstandard = 0.0036764) and Rsample was calculated as follows:

$$R_{sample} = \, (\delta^{15} N + { 1}) \, \times R_{standard} = {{{}^{15}{\text{N}}} \mathord{\left/ {\vphantom {{{}^{15}{\text{N}}} {{}^{14}{\text{N}}}}} \right. \kern-\nulldelimiterspace} {{}^{14}{\text{N}}}}^{15} \times {1}000\permille$$
(2)

The mass balance equation of Providoli et al. (2005) was used to calculate the 15N recovery (15Nrec) of the applied 15N in each N pool as follows:

$$^{15} {\text{N}}_{rec} = \frac{{at.\%^{15} N sample - at.\%^{15} N ref }}{{at.\%^{15} N tracer - at.\%^{15} N ref}} \times \frac{N pool}{{N tracer }} \times 100\%$$
(3)

where 15Nrec is the recovery (%) of the 15N tracer in a soil N pool, calculated based on the tracer proportion found in a sample (at. %15N sample) related to the tracer that was added to the soil (at. %15N tracer). The natural abundance background (at. %15N ref) was measured in each N pool prior to the labelling and is subtracted from both the sample and the added tracer. The tracer fraction in the sample is multiplied by the N pool of the sample (mmol m−2) related to the amount of N (N tracer) that was added to the soil in mmol m−2. A summary of the variables used in the equation can be found in the online resource 1.

N status

Gundersen et al. (1998) and Tietema et al. (1997) used a principal component analysis (PCA) to summarize several variables into one index value per site, referred to as N-status. The N-status values of the investigated sites were correlated with NO3− leaching measured in the 1990s to indicate that some sites were N-saturated (high NO3− leaching) and others N-limited (low NO3− leaching). The PCA done by Tietema et al. (1997) did not include the site AL, therefore we repeated this PCA for all of our investigated sites (SP, KH, GA and AL) and also included the Ysselsteyn (YS) forest site (Koopmans et al. 1996). The variables used in the PCA were: Needles N content, litter fall N content, forest floor N content and the litter fall N flux (Table 2). These variables were measured in the 1990s, therefore this variable can be considered as the initial N-status.

Table 2 Variables used to calculate the initial N-status of the various sites based on the first principal component (PC1) of a principal component analysis (PCA)
Full size table

Data analysis

Differences in 15N recovery between the soil layers (organic vs mineral) and the N treatments (low vs high) were investigated with a nonparametric two sample Wilcoxon test (significance level set at p < 0.1). Changes over time were not tested for significant differences, because at some of the sites a plot-based sampling was applied, or the samples were bulked based on soil type before analysis at the first year after labelling. Index values for the N status of the various sites were calculated with a principal component analysis (PCA) based on four forest N parameters measured in the 1990s (Table 3) by (Tietema et al. 1997; Gundersen et al. 1998). To identify the most important controlling factors that contributed to 15N retention in the mineral soil, a partial least squares regression (PLS-R) analysis was used. This method combines the principles of a principal component analysis (PCA) and multiple linear regression (MLR) to decompose the matrix of several predictor variables X into components (called latent vectors) to predict the response variable Y (15N recovery). The latent vectors describe the effects that cause the changes in the investigated system, for which a group of variables might be responsible. PLS-R has the ability to analyse data with many, noisy, collinear, and even incomplete variables in both X and Y (Wold et al. 2001). The SIMPLS fit algorithm was used and the model was validated by a tenfold cross-validation procedure. Prior to the PLS-R analysis, the X and Y variables were normalised (z-transformed). To select the predictors of 15N recovery, the variable importance for projection (VIP) statistics were used (Chong and Jun 2005; Wold et al. 2001). The soil parameters for each soil layer that were used in the analysis were: C stock, C concentration, ECEC, moisture content, NH4+ concentration, depth, thickness, pH, BD and C:N ratio. The site parameters that were used in the analysis were: N deposition, N status, mean annual precipitation, mean annual temperature and forest stand age. Since δ15N, N concentration and mass of the soil material were used to calculate 15N recovery, these parameters were excluded from the analysis. Two stages of predictor variable selection were performed based on VIP scores > 0.8 (Online resource 2). The dataset was z-transformed prior to the analysis, which gives coefficients that show the importance of each variable based on latent vector 1 (LV1), latent vector 2 (LV2) as well for the combined (LV1 + LV2). All statistical analyses were performed with R (version 3.4.4).

Table 3 Recovery of 15N (%) in the organic and mineral soil layers for each site, sampling year and low and high N deposition
Full size table

Results

15N recovery changes “over time” in the organic and mineral soil

Substantial amounts of 15N were still present in the soil about two decades after 15N addition at all of the investigated forest sites. Total 15N recovery in the soil ranged between 19% in the low N treatment plot at SP and 109% in the high N treatment catchment at GD (Table 3). One year after labelling (in 1993 at SP, KH and GD or in 1997 and 2002 et al.) the organic soil layer was the primary 15N sink across all sites and N deposition levels (Table 3). After approximately two decades (in 2015 or 2016) the 15N recovery levels in the organic soil layer were still higher than in the mineral soil at the low (p = 0.2) and high (p = 0.1) N treatments of SP and at GD (p < 0.0001) (Table 3). For the other two sites the 15N recovery had become higher in the mineral soil than in the organic soil, although a significant difference was found only for the high N treatments at AL (p = 0.006) and at KH (p = 0.01).

Long-term N treatment effects on 15N recovery

The long-term N treatment effects on 15N recovery in various soil layers were investigated at SP, KH and AL based on a comparison of 15N in the N addition and control plots. Figure 1 compares the 15N recovery levels between the different N treatments, per soil layer for each site separately. At SP (Fig. 1a) the 15N recovery levels were higher with the high N treatment in comparison with the low N treatment in all soil layers, although only significantly different in the organic layer (p = 0.008) and in the 0–10 cm soil layer (p = 0.06). At AL (Fig. 1b) the N treatment differences were inconsistent and not significant at any soil layer down the profile. At KH (Fig. 1c) the 15 N recovery levels were higher with the low N treatment in comparison with the high N treatment in all soil layers. The differences were significant in the organic layer (p = 0.008) and in the 25–50 cm soil layer (p = 0.06). GD (Fig. 1d) did only have 15N additions in the high N treatment catchment.

Fig. 1
figure 1

Average (± se, n = 5 at SP and KH, n = 21 et al. and n = 52 at GD) recovery of 15N (%) at Speuld (a), Alptal (b), Klosterhede (c) and Gårdsjön (d), ~ 2 decades after 15N tracer addition with low (red) and high (blue) N deposition levels in the organic soil layer and the 0–10 cm, 10–25 cm and 25–50 soil layers of the mineral soil. At Gårdsjön there was only a high N deposition level and the soil depth layers were measured until 100 cm depth. Significant differences between the N treatments are given at levels: ***p < 0.01, **p < 0.05 or *p < 0.1

Full size image

N-status

The variables included in the PCA analysis were chosen because they accounted for 76% of the total standardized variance in the dataset from the 1990s based on PC1, giving a strong correlation (R2 = 0.95) between the obtained N-status value and NO3− leaching (Fig. 2). The NO3− values for SP, KH, YS, GD are from Tietema et al. (1997) and for AL from Schleppi et al. (1999). The correlation between NO3− leaching and N-status shows that forest sites with a high N-status value, such as YS and SP, have high NO3 leaching (> 25 kg N ha−1 year−1) and are considered N saturated as defined earlier by Aber et al. (1989). Sites with a low N-status value, such as GD, AL and KH had low NO3− leaching and thus require more N than available (N-limited) or has enough N available but is not saturated (N sufficient).

Fig. 2
figure 2

Correlation between N-status and NO3− leaching as indicator for N limitation (GD, AL, KH) or N saturation (SP and YS) of a forest ecosystem. N-status values are based on the first principal component (PC1) of a principal component analysis (PCA), calculated with the variables given in Table 2. NO3− values for SP, KH, YS, GD, are from Tietema et al. (1997). NO3− values for AL are from Schleppi et al. (1999)

Full size image

Drivers of 15N recovery detected by PLS-R.

The PLS-R coefficients (R2) of the variables that significantly contributed to the prediction of 15N recovery in the organic and mineral soil layers are given in Table 4. In the organic soil layer, C stock had the highest impact (R2 = 0.72) on 15N recovery, followed by layer thickness (R2 = 0.51), temperature (R2 = − 0.34) and N-status (R2 = − 0.22). The cumulative R2 value for the predictors (RX2) was 0.89 and the cumulative R2 value for the response (RY2) was 0.59 (Table 4). In the mineral soil, C content (R2 = 0.46) had the highest impact on 15N recovery, followed by C stock (R2 = 0.35), moisture content (R2 = 0.26), mean annual temperature (R2 = -0.24), bulk density (R2 = -0.19), forest stand age (R2 = 0.18) and mean annual precipitation (R2 = 0.13). The cumulative RX2 value in the mineral soil layers was 0.82 and the cumulative RY2 value was 0.52 (Table 4). The bi-plots in Fig. 3 visualize the formation of clusters of the individual sample points (x scores) in combination with the direction of the variables (y-loadings) based on the PLS-R results. In the organic layer, LV1 separated the sites with negative values (YS, SP, KH and AL) from GD (and a part of AL) with predominantly positive values. The addition of LV2 contributed to some extend to a separation of all individual sites, although the variation of the GD samples along LV2 was large, encompassing all the samples of the other sites. Altogether, the separation was based on a combined (LV1 + LV2) positive impact of C stock and layer thickness and a negative impact of temperature and N-status. In the mineral soil layer, there was a larger number of variables that contributed to the LV’s than in the organic soil layer (Fig. 3 and Table 4). In the mineral soil layers there were three main clusters formed based on LV1 that separated AL, GD and a cluster of YS, SP and KH. AL had the highest positive values based on LV1. LV2 did not contribute to a further separation of the sites, it showed a large sample variation within each site. The fit of the prediction of the 15N recovery was similar for the organic and mineral soil layers (R2 = 0.59 and R2 = 0.52 respectively) (Fig. 4). The corresponding equations based on the z-transformed data in the organic and mineral soil are given by:

  • 15N recovery organic (%) = 0.72 °C stock (kg ha−1) + 0.51 Organic layer thickness (cm) − 0.34 Temperature (°C) − 0.22 N-status.

  • 15N recovery mineral (%) = 0.46 °C concentration (%) + 0.35 °C stock (kg ha−1) + 0.26 Moisture content (g/g soil) − 0.24 BD (g cm−3) − 0.19 Temperature (°C) + 0.18 Stand age (years) + 0.13 Precipitation (mmyr−1).

Table 4 PLS-R coefficients (R2) that significantly contributed to the prediction of 15N recovery, based on 2 latent vectors and the combined effects of these two latent vectors in the organic and mineral soil
Full size table
Fig. 3
figure 3

Bi-plots with x scores and y loadings. The scores are colour labelled by site and shown with different symbols for low N deposition (o) and high N deposition (+). The selected loadings are the predictors that significantly contributed to 15N recovery. The results are given for LV 1 and LV 2 for the organic soil (a) and the mineral soil (b), which includes 2 (AL), 3 (SP, KH, YS) or up to 4 (GD) soil depth layers per core. The closer the x-scores are in the direction of the y-loadings, the more important the variables are for the prediction of 15N recovery at these sample points. The angle between the vectors of the y-loadings show how the variables are correlated with each other. If the angle between the vectors is close to 90°, they are not likely to be correlated. At an angle of 180°, the variables are negatively correlated

Full size image
Fig. 4
figure 4

Partial least squared regression (PLS-R) model output and 1:1 line. The predicted values of 15N recovery are given on the y-axis and the measured 15N recovery values on the x-axis for the organic (a) and mineral (b) soil. The 15N recovery (%) data is z transformed prior to the analysis

Full size image

Discussion

Long-term 15N recovery in the organic and mineral soil

We hypothesized that the mineral soil would be a larger 15N sink than the organic soil on the long-term. Our results were only partly consistent with this hypothesis because approximately 20 years after labelling, the 15N recovery levels were higher in the mineral soil layer than in the organic soil layer in 2 out of 4 sites (Table 3). This contrasts with the results one year after labelling, where the organic soil layer was the dominant 15N sink at all of the investigated sites. The results indicate that 15N may move from the organic soil layer towards the mineral soil layer over time via downward transport processes (Cheng et al. 2018). The immobilization of 15N in the organic soil layer as well as the transport processes of 15N from the organic soil layer to the mineral soil layer will be influenced by specific soil and site characteristics. The increase in the total soil 15N recovery (organic + mineral) over time at most sites and deposition levels may partly be related to an initial uptake of the added 15N in the vegetation pools. This uptake depended on the competition between the soil microbes and the trees for N, and the availability of incoming N relative to mineralized organic N. As time passes, the amount 15N in the soil is expected to increase due to input of 15N from tree litter. Additionally, with an increasing N availability in the soil, the re-translocation of N from litter to needles during senescence tends to decline (Kjønaas and Stuanes 2008). This change in re-translocation patterns will increase the input of 15N to the soil pool both from above-and below ground litter. At GD, the increases over time in the total soil pool exceeded a recovery of 100%, suggesting a substantial overestimation of the recovery rates. This may be due to the re-calculation of diagnostic horizons to the given soil depths and/or the re-calculations of BD. Additionally, these overestimations are most likely related to a varying topography combined with a large catchment size, which could have reduced the exact pool size estimations. The measurement uncertainties were expected to be smaller at the other sites, due to smaller catchment/plot sizes with less variability.

Long-term N treatment effects on 15N retention in the soil

The second hypothesis that higher N deposition levels would reduce the total 15N recovery in the soil was not supported by the current study. The mechanism behind the hypothesis was that especially over longer time scales, higher N deposition would promote N saturation and N loss from the soil via leaching or gaseous emission (Aber et al. 1989; Gundersen et al. 1998; Templer et al. 2012). At KH, the 15N recovery levels in the organic layer and in the deepest (25–50 cm) mineral soil layer were significantly lower in the high N plot (55 kg N ha−1 year−1) than in the low N plot (20 kg N ha−1 year−1). At first sight these results seemed to support the second hypothesis. However, this hypothesis was primarily based on the results from YS (Wessel et al. 2013 and Wessel et. Unpublished) that showed higher N leaching losses with higher N deposition levels. For KH the mechanism behind these results seemed to differ from our second hypothesis. There were almost no N losses measured at KH (Gundersen et al. 1998; Gundersen unpublished). The trees were found to be more competitive for N and had a higher 15N recovery in the high N plot than in the low N plot on the short term (Gundersen et al. 1998). This effect was still evident in 2009 (after 17 years) where the 15N recovery in trees remained higher in the high N plot than in the low N plot (Gundersen unpublished). The long-term 15N recovery levels in the soil at KH (L:52% H:33%) were higher than the recovery levels in the soil than at YS forest (L:27% H:18%). This indicated higher leaching losses at YS than at KH (Wessel et al. 2013; Wessel et al. unpublished). The lower N-status at KH (with minor N leaching) as compared to YS (major N leaching loss) at the start of the experiment (Fig. 2), showed that these forests were at different stages along an N saturation range at the onset of the experiments, which is supported by the higher retention of incoming 15N at KH compared to YS.

Contrary to our second hypothesis, the results from AL did not show higher 15N recovery in the low N addition relative to the high N addition, as higher N deposition levels did not change 15N recovery in any of the soil layers (Fig. 1b). Furthermore, more than 60% of the added 15N remained in the total soil pool in both catchments. The high soil 15N recovery combined with the significant lower soil C/N ratio in the high N catchment indicates that the N status at AL is increasing, but the capacity of the soil sink is not yet saturated for increased 15N losses to occur. Parallel to these high 15N recovery levels in the organic and mineral soil layers for both N treatments, continuous NO3− leaching was observed over a period of 20 years by Schleppi et al. (2017). The leaching of 15NO3− disappeared, however, a few months after the end of the 15N addition. This suggested that the leached NO3− was mainly derived from newly deposited N and that this leaching was hydrologically driven (Schleppi et al. 2004). The same situation was observed in the runoff from the GD catchment, where the dissolved inorganic N in runoff during the year of the 15N addition largely reflected the 15N level of the incoming N, indicating that the leached NO3− came predominantly from sprinkled N input and precipitation. At GD, only 1.1% of the incoming N was lost during the year of the tracer addition, whereas the cumulative loss of tracer N over a 10-year period was 3.9% as DIN and 1.1% as DON (Kjønaas and Wright 2007). These examples of kinetic N saturation (Lovett and Goodale 2011) indicate that the “old” labelled 15N remained in the forest system while newly deposited unlabelled N was lost from the system.

The results of SP were not supporting our second hypothesis, because the 15N recovery levels in the organic soil layer (p = 0.008) and the 0–10 cm soil layer (p = 0.06) significantly decreased after decreasing N deposition levels. SP differed from the other sites regarding the duration of the low N treatment, which was much shorter. As the roof was removed in 1995, the N deposition was actually high in the last 20 years prior to sampling. Furthermore, the N and C concentrations in the mineral soil at SP were lower than at the other sites (Online resource 1). This reflects a generally low potential for 15N incorporation in the mineral soil, which was also found to be lower at SP relative to AL, KH and GD.

There are several reasons why support for our second hypothesis may not be evident, potentially due to differences between the N status of these sites. Additionally, uncertainties in the results are related to the use of a reference (atom%15Nref) natural abundance background level which was kept constant with time as well as with N treatment. This may introduce a bias in the estimates of the long- term recovery of 15N with high and low N inputs. Additionally, it is likely that the N treatment affects 15N fractionation during N transformation processes, which in turn affect the 15N abundance. An error due to fractionation processes would increase with time, leading to an overestimation of 15N recoveries, because the N lost by leaching or denitrification is lighter 14N. The exact effects are unknown but are assumed to have had a smaller impact on the sites with low N losses (AL, GD and KH) than at sites with higher N losses (SP and YS). A better reference would potentially be the use of an adjacent long-term N treated but unlabelled plot or catchment, which was not available. Further, to be able to understand the treatment effects at the sites more fully, whole ecosystem studies including all N pools, N transformation processes and leaching losses seem necessary.

Drivers of 15N recovery detected by PLS-R

The PLS-regression analysis was used to identify the soil and site characteristics that had the highest impact on the prediction of the 15N recovery in the organic and the mineral soil layers. These results were generally in line with our third hypothesis, which explained that variables related to C dynamics would have the largest impact on 15N recovery in the soil. In the organic soil layer (Fig. 3a) C stocks (kg ha−1) had the highest impact on 15N recovery. It is likely that most of the accumulated 15N in the organic soil layer was immobilized by microbial and plant biomass (Kopáček et al. 2013). The thickness of the organic layer was clearly correlated with C stock and the immobilization processes. Furthermore, a thick organic layer enables re-immobilization of N in deeper parts of the organic soil layer instead of reaching the mineral soil. The highest 15N recovery levels in the organic soil layers were found at GD, which on average had the thickest organic layer amongst the sites (Online resource 1). This, along with the low N status, which is characterized by a tight N cycling (Fig. 2), most likely contributed to the high 15N recovery levels in the organic soil layer relative to the mineral soil layers at the GD site. The 15N recovery levels in the organic soil layers of sites with a relatively low (initial) N-status (GD, AL and KH) were higher compared to the sites (SP and YS) that were (initially) N-saturated (Aber et al. 1989; Tietema et al. 1997; Gundersen et al. 1998). These results point towards biologically driven processes of NO3− leaching (Kopáček et al. 2013). At AL and GD NO3− leaching was also expected to be hydrologically driven. The low permeability of the Gleysols at AL or the shallow soils and bedrock outcrops at GD, combined with a high rainfall intensity and a varying topography most likely contributed to that (Kjønaas et al. 1998; Providoli et al. 2005). The negative contribution of temperature to 15N recovery in the organic soil layer was in line with the temperature sensitivity of SOM decomposition and respiration (Davidson and Janssens 2006).

In the mineral soil layers (Fig. 3b) all of the variables that significantly contributed to 15N recovery were directly or indirectly linked to SOM availability. The positive correlations of C concentration and C stock with15N recovery were obviously a direct effect of SOM accumulation in the mineral soil. Higher input of organic material, lower decomposition and/or higher retention of incoming SOM in the mineral soil increased its C and N concentration. The fixed soil layer depths across the sites excluded thickness as a predictor for 15N recovery in the mineral soil. The positive correlations of 15N recovery with moisture content and mean annual precipitation might be related to a decrease in soil respiration in some of the frequently waterlogged depressions at the AL and GD catchments (Kjønaas et al. 1998; Providoli et al. 2005). Although at GD, the 15N leached from shallow outcrops as well as a hydrological driven leaching of 15N was more likely to increase the input of 15N to these depressions. The moisture content data are representing a high level of uncertainty as the data comes from single point measurements and thus will strongly depend on the season and weather conditions prior to and during the sampling period. Mean annual temperature had a negative impact on 15N recovery potentially related to respiration rates, as previously mentioned. This negative impact is expected to be higher for systems with higher decomposition rates. The negative impact of temperature was larger in the organic soil (R2 = − 0.34) than in the mineral soil (R2 = − 0.19). When decomposition rates would increase, for example due to climate change, the effects on 15N recovery are expected to be larger in the organic soil than in the mineral soil. Forest stand age had a positive effect on 15N recovery, which may be due a possibly larger belowground biomass in older forest (Peichl and Arain 2006), potentially contributing to higher 15N root litter input, as well as higher SOM accumulation and thus higher 15N recovery levels. Additionally, older forests have had more time to recover from harvesting disturbances which may influence the SOM accumulation positively (Nave et al. 2010). Bulk density was negatively correlated with 15N recovery in the mineral soil, possibly caused by higher 15N recovery levels in the upper layers of the mineral soil with a low bulk density and lower 15N recovery levels in the deeper layers with a high bulk density. However, depth did not significantly correlate with 15N recovery due to the concurrent increase in thickness (0–10 cm, 10–25 cm, 25–50 cm).

The PLS-R model gave useful insights into the relative impact of the available soil and site variables on the 15N recovery. In both the organic and mineral soil layers, the drivers that significantly contributed to 15N recovery were strongly related to C availability and seemed mostly biologically driven. However, the absence of data on expected key driving factors, as well as limited numbers of field sites limit the model results. Soil parameters such as aggregate stability, soil texture, soil redox potential, N mineralization, SOM decomposition and soil microbial community structure (Lal et al. 2015; Wiesmeier et al. 2019) should be included to improve the explanatory value of the PLS-R analysis. Furthermore, we think that the larger sample sizes of the catchments (AL = 21, GD n = 52) in comparison with the plots (n = 5) might have had an impact on the pls-r analysis by driving the results in the direction of AL and GD.

Conclusion

The main goal of this study was to investigate the long-term fate of 15N in coniferous forest soils. The results of this research show that forest soils are important pools for longer-term N retention. Between 19 and 100% of the once added 15N were found to accumulate in the forest soil over ~ two decades. In half of the investigated sites the mineral soil layer became a larger 15N sink than the organic soil layer with time. At all of the investigated sites increased N deposition did not contribute to lower 15N recovery levels in the soil, possibly due to limited N leaching losses that was predominantly hydrologically driven. The drivers that contributed to the prediction of 15N recovery were mostly related to SOM accumulation and turnover processes (indicated with C stock, organic layer thickness, C concentration, forest stand age, temperature, precipitation and moisture content). The initial N-status of the sites seemed to have an impact on the longer-term 15N recovery in the organic soil layer. The current study should be extended with measurements of additional environmental drivers to improve the 15N recovery prediction model. However, it is clear that long-term increased N deposition had a minor impact on N leaching, due to the high buffer capacity of these forests to store N in both the organic and mineral soil.