The long-term fate of deposited nitrogen in temperate forest soils

Increased anthropogenic nitrogen (N) inputs can alter the N cycle and affect forest ecosystem functions. The impact of increased N deposition depends among others on the ultimate fate of N in plant and soil N pools. Short-term studies (3–18 months) have shown that the organic soil layer was the dominant sink for N. However, longer time scales are needed to investigate the long-term fate of N. Therefore, the soils of four experimental forest sites across Europe were re-sampled ~ 2 decades after labelling with 15N. The sites covered a wide range of ambient N deposition varying from 13 to 58 kg N ha−1 year−1. To investigate the effects of different N loads on 15N recovery, ambient N levels were experimentally increased or decreased. We hypothesized that: (1) the mineral soil would become the dominant 15N sink after 2 decades, (2) long-term increased N deposition would lead to lower 15N recovery levels in the soil and (3) variables related to C dynamics would have the largest impact on 15N recovery in the soil. The results show that large amounts of the added 15N remain in the soil after 2 decades and at 2 out of 4 sites the 15N recovery levels are higher in the mineral soil than in the organic soil. The results show no clear responses of the isotopic signature to the changes in N deposition. Several environmental drivers are identified as controlling factors for long-term 15N recovery. Most drivers that significantly contribute to 15N recovery are strongly related to the soil organic matter (SOM) content. These findings are consistent with the idea that much of the added 15N is immobilized in the SOM. In the organic soil layer, we identify C stock, thickness of the organic layer, N-status and mean annual temperature of the forest sites as most important controlling factors. In the mineral soil we identify C stock, C content, pH, moisture content, bulk density, temperature, precipitation and forest stand age as most important controlling factors. Overall, our results show that these temperate forests are capable of retaining long-term increased N inputs preferably when SOM availability is high and SOM turnover and N availability are low.


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 (NH 4 ? ) and nitrate (NO 3 -) in excess of total combined plant and microbial nitrogen demand. N saturation would be manifested as NH 4 ? accumulation and/or NO 3 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 N 2 O.
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 15 Nenriched NH 4 ? and/or NO 3 on the forest floor. The use of 15 N enriched deposition enables monitoring of incoming 15 N in various vegetation and soil N pools. Increased N deposition could have an impact on the fate of 15 N in these N pools or leaching or gaseous N losses. When increased N deposition would result in increased 15 N 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 15 N 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 15 N in the soil as a result of 15 N leaching (Tietema et al. 1998;Nadelhoffer et al. 1999;Templer et al. 2012). These studies reflect 15 N recovery after a relatively short period (3-18 months). Longer time scales are needed to investigate long-term turnover and stabilization processes of 15 N in forest systems.
A few studies have investigated 15 N 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 15 N sink 6-9 years after 15 N 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 15 N sink already after 1 year (Goodale 2017). From a resampling of the Ysselsteyn forest in the Netherlands (Wessel et al. 2013) 19 years after 15 N addition, it was shown that 15 N 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 15 N recovery in the mineral soil layers increases with time. Yet, there are several other mechanisms contributing to 15 N 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 NO 3 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, 15 N 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 .
Due to the limited number of studies on the longterm fate of deposited N in temperate forest systems, we re-investigated several 15 N labelling and NITRogen saturation EXperiments (NITREX) experiments approximately 20 years after labelling (Wright and van Breemen 1995). We focused on 15 N recovery in the organic and mineral soil layers, as these soil N pools were the dominant 15 N sinks in previous studies (Nadelhoffer et al. 1999). Our first hypothesis was that the 15 N 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 15 N 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 15 N 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 15 N recovery. Third, we hypothesized that variables such as, C stock, C content and layer thickness, would have a large impact on 15 N 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 15 N 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.

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 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 ).
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 15 N input under ambient N deposition, except for GD where the control catchment was not labelled with 15 N. 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 (NH 4 NO 3 ) 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 NO 3 and NH 4 ? 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 15 NH 4 and/or 15 NO 3 for a period of 12 months in 1992/1993 at SP, KH, GD and in 1995/1996at AL (high N catchment) and in 2000(K 15 NO 3 ) and 2002/2003 Cl) at the AL control catchment (Table1).

Re-sampling and laboratory analysis
More than 20 years after 15 N 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 Table 1 Ambient N deposition, treatment N flux, type of 15 N addition, first year of 15 N 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  Koopmans et al. (1996) 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 m 2 and at GD the catchment had a size of 5250 m 2 . 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 . 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 ), 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 15 N abundance of the samples. Reference gas (high-purity N 2 gas) was calibrated to atmospheric N 2 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 ? , Mg 2? , Ca 2? , Fe 3? , Al 3? ), effective cation exchange capacity (ECEC) and NH 4 ? concentrations in the mineral soil samples were measured in barium chloride (BaCl 2 ) extracts of 0.1 M BaCl 2 with a soil to BaCl 2 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 NH 4 ? and NO 3 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 cm 3 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: where d 15 N was expressed as permille delta 15 N (d 15 N %) of total N. R is the molar ratio of 15 N/ 14 N, with atmospheric N 2 used as standard (R standard-= 0.0036764) and R sample was calculated as follows: The mass balance equation of Providoli et al. (2005) where 15 N rec is the recovery (%) of the 15 N tracer in a soil N pool, calculated based on the tracer proportion found in a sample (at. % 15 N sample) related to the tracer that was added to the soil (at. % 15 N tracer). The natural abundance background (at. % 15 N 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 NO 3 leaching measured in the 1990s to indicate that some sites were N-saturated (high NO 3 leaching) and others N-limited (low NO 3 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.

Data analysis
Differences in 15 N 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 15 N 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 ( 15 N 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 15 N 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, NH 4 ? 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 d 15 N, N concentration and mass of the soil material were used to calculate 15 N 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).

Results
15 N recovery changes ''over time'' in the organic and mineral soil Substantial amounts of 15 N were still present in the soil about two decades after 15 N addition at all of the investigated forest sites. Total 15 N 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 The values of the varaibles for SP, KH, YS and GD are from Tietema et al. (1997). The values fof the variables for AL are from Schleppi et al. (1999)  The short-term 15 N recovery levels (1993, 1997 and 2002) were obtained from Tietema et al. 1997 (Table 3). One year after labelling (in 1993at SP, KH and GD or in 1997and 2002 the organic soil layer was the primary 15 N sink across all sites and N deposition levels (Table 3). After approximately two decades (in 2015 or 2016) the 15 N 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 15 N 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 15 N recovery
The long-term N treatment effects on 15 N recovery in various soil layers were investigated at SP, KH and AL based on a comparison of 15 N in the N addition and control plots. Figure 1 compares the 15 N recovery levels between the different N treatments, per soil layer for each site separately. At SP (Fig. 1a) the 15 N 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 15 N additions in the high N treatment catchment.

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 (R 2 = 0.95) between the obtained N-status value and NO 3 leaching (Fig. 2). The NO 3 values for SP, KH, YS, GD are from Tietema et al. (1997) and for AL from Schleppi et al. (1999). The correlation between NO 3 leaching and N-status shows that forest sites with a high N-status value, such as YS and SP, have high NO 3  Table 2. NO 3 values for SP, KH, YS, GD, are from Tietema et al. (1997). NO 3 values for AL are from Schleppi et al. (1999) 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 NO 3 leaching and thus require more N than available (N-limited) or has enough N available but is not saturated (N sufficient).
Drivers of 15 N recovery detected by PLS-R.
The PLS-R coefficients (R 2 ) of the variables that significantly contributed to the prediction of 15 N recovery in the organic and mineral soil layers are given in Table 4. In the organic soil layer, C stock had the highest impact (R 2 = 0.72) on 15 N recovery, followed by layer thickness (R 2 = 0.51), temperature (R 2 = -0.34) and N-status (R 2 = -0.22). The cumulative R 2 value for the predictors (R X 2 ) was 0.89 and the cumulative R 2 value for the response (R Y 2 ) was 0.59 (Table 4). In the mineral soil, C content (R 2 = 0.46) had the highest impact on 15 N recovery, followed by C stock (R 2 = 0.35), moisture content (R 2 = 0.26), mean annual temperature (R 2 = -0.24), bulk density (R 2 = -0.19), forest stand age (R 2 = 0.18) and mean annual precipitation (R 2 = 0.13). The cumulative R X 2 value in the mineral soil layers was 0.82 and the cumulative R Y 2 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 (yloadings) 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 15 N recovery was similar for the organic and mineral soil layers (R 2 = 0.59 and R 2 = 0.52 respectively) (Fig. 4). The corresponding equations based on the z-transformed data in the organic and mineral soil are given by: The coefficients were z-transformed to be able to directly compare the impact of the individual variables on the prediction of the 15 N recovery. The cumulative R 2 is given for the predictors and X and Y

Long-term 15 N recovery in the organic and mineral soil
We hypothesized that the mineral soil would be a larger 15 N 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 15 N recovery levels were higher in the  Fig. 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 15 N 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 15 N 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  (Table 3). This contrasts with the results one year after labelling, where the organic soil layer was the dominant 15 N sink at all of the investigated sites.
The results indicate that 15 N 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 15 N in the organic soil layer as well as the transport processes of 15 N 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 15 N recovery (organic ? mineral) over time at most sites and deposition levels may partly be related to an initial uptake of the added 15 N 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 15 N in the soil is expected to increase due to input of 15 N 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 retranslocation patterns will increase the input of 15 N 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 15 N retention in the soil
The second hypothesis that higher N deposition levels would reduce the total 15 N 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 15 N 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 unpublished). The trees were found to be more competitive for N and had a higher 15 N recovery in the high N plot than in the low N plot on the short term ). This effect was still evident in 2009 (after 17 years) where the 15 N recovery in trees remained higher in the high N plot than in the low N plot (Gundersen unpublished). The long-term 15 N 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 15 N at KH compared to YS. Contrary to our second hypothesis, the results from AL did not show higher 15 N recovery in the low N addition relative to the high N addition, as higher N deposition levels did not change 15 N recovery in any of the soil layers (Fig. 1b). Furthermore, more than 60% of the added 15 N remained in the total soil pool in both catchments. The high soil 15 N 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 15 N losses to occur. Parallel to these high 15 N recovery levels in the organic and mineral soil layers for both N treatments, continuous NO 3 leaching was observed over a period of 20 years by Schleppi et al. (2017). The leaching of 15 NO 3 disappeared, however, a few months after the end of the 15 N addition. This suggested that the leached NO 3 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 15 N addition largely reflected the 15 N level of the incoming N, indicating that the leached NO 3 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 15 N 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 15 N 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 15 N 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% 15 Nref) 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 15 N with high and low N inputs. Additionally, it is likely that the N treatment affects 15 N fractionation during N transformation processes, which in turn affect the 15 N abundance. An error due to fractionation processes would increase with time, leading to an overestimation of 15 N recoveries, because the N lost by leaching or denitrification is lighter 14 N. 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 15 N 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 15 N 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 15 N recovery in the soil. In the organic soil layer (Fig. 3a) C stocks (kg ha -1 ) had the highest impact on 15 N recovery. It is likely that most of the accumulated 15 N 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 15 N 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 15 N recovery levels in the organic soil layer relative to the mineral soil layers at the GD site. The 15 N 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 NO 3 leaching (Kopáček et al. 2013). At AL and GD NO 3 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 Providoli et al. 2005). The negative contribution of temperature to 15 N 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 15 N recovery were directly or indirectly linked to SOM availability. The positive correlations of C concentration and C stock with 15 N 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 15 N recovery in the mineral soil. The positive correlations of 15 N 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 Providoli et al. 2005). Although at GD, the 15 N leached from shallow outcrops as well as a hydrological driven leaching of 15 N was more likely to increase the input of 15 N 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 15 N 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 (R 2 = -0.34) than in the mineral soil (R 2 = -0.19). When decomposition rates would increase, for example due to climate change, the effects on 15 N recovery are expected to be larger in the organic soil than in the mineral soil. Forest stand age had a positive effect on 15 N recovery, which may be due a possibly larger belowground biomass in older forest (Peichl and Arain 2006), potentially contributing to higher 15 N root litter input, as well as higher SOM accumulation and thus higher 15 N 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 15 N recovery in the mineral soil, possibly caused by higher 15 N recovery levels in the upper layers of the mineral soil with a low bulk density and lower 15 N recovery levels in the deeper layers with a high bulk density. However, depth did not significantly correlate with 15 N 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 15 N recovery. In both the organic and mineral soil layers, the drivers that significantly contributed to 15 N 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 longterm fate of 15 N 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 15 N were found to accumulate in the forest soil over * two decades. In half of the investigated sites the mineral soil layer became a larger 15 N sink than the organic soil layer with time. At all of the investigated sites increased N deposition did not contribute to lower 15 N 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 15 N 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 15 N recovery in the organic soil layer. The current study should be extended with measurements of additional environmental drivers to improve the 15 N 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.