Our study site is located in the Lüneburg Heath nature reserve in NW Germany (Lower Saxony, Schneverdingen, 53°15′N, 9°58É, 105 m a.s.l.), the site of the largest complex of heathlands (about 5,000 ha) in NW Germany. The study sites is characterised by Pleistocene sandy deposits with a weak relief and nutrient-poor podzol soils. Topsoil-pH_H2O values range between 3.3 and 3.5 (Härdtle et al., 2006). The climate is of a humid sub-oceanic type (mean precipitation: 811 mm year⁻¹, annual mean temperature: 8.4°C; Niemeyer et al., 2005)., Heathlands of the nature reserve are dominated by the dwarf shrub Calluna (cover >80%). The soil surface is covered by mosses (mean cover 93%), with Hypnum judlandicum and Pleurozium schreberi being the prevailing species. The current background atmospheric N deposition in the study area amounts to 23 kg ha⁻¹ year⁻¹ (von Oheimb et al., 2010 for more detailed information on the quantity of different N-forms deposited in NW Germany and in the study area see Gauger et al. (2000) and Niemeyer et al. (2005), respectively).

The present study made use of the experimental setup of Friedrich et al. (2011b) established in 2007. The experiment comprises a total of seven randomly selected plots, each of which was divided into two subplots (2 × 4 m in size). At the beginning of the experiment (June 2007), in one subplot (per plot) a ¹⁵N tracer was added in liquid form (pulse labeling with 100 mg ¹⁵NH415NO₃ (98 at.%) using a special spray bottle that allowed for a fine nebulisation of the liquid applied), and the corresponding subplot was used to determine the natural ¹⁵N abundance (these subplots are henceforth referred to as labelled and non-labelled subplots, respectively). Labelled and non-labelled subplots were separated by a buffer zone 1 m in width to avoid cross contamination after ¹⁵N tracer addition (for a detailed description of the ¹⁵N application procedure see Friedrich et al., 2011b).

All analyses performed during the present study refer to a time span of 7 years, with all samples collected between 2009 and 2015 (with the exception of one plot, which was sampled from 2008 to 2015, as well as litter and leachate samples; see description below). To present a complete overview of N flows and partitioning since the setup of the experiment in 2007, we included analyses of Friedrich et al. (2011b) in the present study (i.e., analyses from 2007 to 2008). In this way, we were able to document compartment-related shifts in ¹⁵N abundance, enrichment and recovery for a 9-year timespan (2007–2015).

To quantify N allocation patterns in the biomass and soil we sampled the following ecosystem compartments annually: above-ground biomass of Calluna (current year and 1- to 2-year-old shoots), moss layer, and the soil O-, A- and B-horizons (with the latter two representing the albic and spodic horizons, respectively, both of which are typical of podzol soils; nomenclature according to the IUSS Working Group WRB, 2015) in labelled and non-labelled subplots in each August from 2008 to 2015, corresponding to the flowering period of Calluna. Twenty randomly chosen current year shoots as well as 1- to 2-year-old shoots of Calluna from the top of randomly chosen Calluna plants were cut with scissors and were pooled to one sample per subplot. Six randomly chosen squares of the moss layer (2 × 2 cm in size) were sampled and bulked to one sample per subplot. Organic soil samples (O-horizon; mean thickness: 2.7 cm) were cut with knives (six squares of 2 × 2 cm of the entire horizon depth per subplot) and separated from fresh fallen litter. Mineral soil samples were taken from the A- and B-horizon with a soil auger (Pürckhauer, EcoTech, Bonn, Germany) from six randomly chosen spots per subplot and bulked to one sample each mean thickness of A- and B-horizons: 7.1 cm and 7.5 cm, respectively). The annual litter production of Calluna was quantified in labelled and non-labelled subplots by means of fine-mesh nets stretched under one randomly selected plant per subplot and sampled during winter between the end of February and the beginning of April (depending on the duration of snow cover) from 2008 to 2016.

To obtain an overview of the quantities of ¹⁵N accumulated in older Calluna twigs and stem wood (age >2 years) and Calluna below-ground biomass (i.e., roots), but to avoid an interference of plant growth in plots, these compartments were sampled only twice, in August 2014 and 2015. To this end, six randomly chosen stems of Calluna per subplot were cut with clippers at ground level (excluding current year and 1- to 2-year-old shoots). For root analysis, we sampled the O-horizon at six randomly selected sites per subplot (3 × 3 cm samples from across the entire organic layer).

Soil and root samples were stored in a freezer (−18°C) until analysis. Prior to analysis, above-ground biomass and soil samples were air-dried, and the latter one sorted and sieved (2 mm). Root samples were washed using a 0.2 mm sieve to remove soil particles and then dried (48 h at 80°C). All samples were ground with a mixer mill (MM 400, Retsch, Haan, Germany) and re-dried at 105°C for weighing. Total N contents and δ¹⁵N were determined with a continuous flow elemental analyser-isotopic ratio mass spectrometer (vario EL cube, Elementar, Hanau, Germany, coupled to an Isoprime IRMS, Isoprime Ltd., Cheadle Hulme, UK; note that total N analyses of soils hence included microbial biomass N).

N losses via leaching were quantified by means of four pairs of lysimeters, one of which was labelled with 100 mg ¹⁵NH415NO₃ m⁻² in 2007 (analogous to the labelling of subplots), and one which received a corresponding amount of deionised water (non-labelled lysimeter; for a detailed description of the setup of lysimeters see Friedrich et al., 2011b). Since the collection and analysis of leachate was very laborious and ¹⁵N leaching losses (in terms of inorganic N and DON) proved to be negligible at the beginning of the experiment, leaching losses of ¹⁵N were only quantified on 20 occasions in the period from June 2014 to May 2015 (weekly in June 2014, once every 2 weeks from July to August 2014 and monthly from September 2014 to June 2015). The leachate which accumulated between the sampling occasions was extracted with a temporarily connected vacuum pump. The quantity of leachate per lysimeter was recorded and a subsample of 500 ml was filtered and stored in a freezer (−18°C) for analysis. Samples were carefully defrosted and analysed with a UV/Vis scanning spectrophotometer (UV-3100PC, VWR, Radnor, PA, USA) using the method of Hagedorn and Schleppi (2000). In order to determine δ¹⁵N values, samples were prepared following the diffusion method of Sebilo et al. (2004) modified according to Friedrich et al. (2011b).

To calculate ¹⁵N tracer recovery, the mass of N pools in each ecosystem compartment must be known. Therefore, the total N content of a compartment in labelled subplots was multiplied by its mean N pool mass. Pool masses for the moss and soil compartments were taken from Friedrich et al. (2011b). Calluna plants were 10–12 years old at the beginning of the experiment (i.e., early mature phase of life-history; Gimingham, 1972). The continuous increment of Calluna above-ground biomass (i.e., total above-ground biomass and total biomass of current year and 1- to 2-year-old shoots) was calculated on the basis of age-dependent growth rates of Calluna in the nature reserve, derived from our own measurements (based on biomass harvests at six randomly chosen sites adjacent to the experimental plots and comparisons with data from Friedrich et al., 2011b) in the Lüneburg Heath and complemented by literature data (Gimingham, 1972; Milne et al., 2002). Calluna below-ground biomass was also quantified in mature stands adjacent to our plots (at six randomly selected sites, 15 soil samples per site taken with a soil auger 8 cm in diameter and to a soil depth of 15 cm; following the procedure for root analyses described above).

N leaching losses were calculated from the total amount of leachate multiplied by the total N contents (N_inorg and DON) of the leachate from labelled lysimeters (Friedrich et al., 2011b).

¹⁵N abundance, ¹⁵N enrichment and ¹⁵N tracer recovery were calculated according to Fry (2006). ¹⁵N abundance and the natural ¹⁵N abundance in labelled and non-labelled subplots, respectively, are represented in the δ notation:

where R_sample and R_standard are the ¹⁵N/¹⁴N ratio of the sample and the N₂ of the atmosphere, which is used as standard, respectively.

¹⁵N enrichment [‰] represents the isotopic enrichment of a sample from a labelled subplot (δ¹⁵N_sample) compared to a sample in a non-labelled subplot (δ¹⁵N_ref):

¹⁵N recovery [g N m⁻²] is the mass of ¹⁵N tracer recovered in the different ecosystem compartments (biomass, soil horizons, and leachate) in labelled subplots and was calculated as follows:

where m_pool is the amount of N in a respective compartment [g N m⁻²], at.%¹⁵N_pool and at.%¹⁵N_ref are the at.%¹⁵N in the N pool of the samples (or in the leachate) of labelled and non-labelled subplots, respectively, and at.%¹⁵N_tracer is the at.%¹⁵N of the added ¹⁵N tracer (Nadelhoffer et al., 2004).

Differences between ¹⁵N abundances of labelled and ¹⁵N natural abundances from non-labelled subplots were tested by means of pairwise comparisons of δ¹⁵N values of each compartment and for each sampling occasion using a Student's t-test (α = 0.05). We applied generalised additive mixed models (GAMM) with a log link function and gamma distribution to assess the effects of ecosystem compartments and time (2007–2015) on ¹⁵N enrichment. Ecosystem compartments were stratified into five main compartments: above-ground biomass of Calluna (including current year shoots, 1- to 2-year-old shoots, and current year litter), moss layer, and O-, A- and B-horizons of the soil. To address the skewed response of the ¹⁵N enrichment data, a gamma probability distribution was preferred, because it retains the structure of the data while accounting for a heteroscedastic error structure (Zuur et al., 2014). We added a value of 0.001 to the ¹⁵N enrichment data to ensure convergence of the fitting algorithm. The non-linear time effect was modelled using a cubic regression spline with four degrees of freedom to allow some complexity in the ¹⁵N enrichment function, while avoiding over-fitting the data (Wood, 2006). We used plot as random factor to account for between-plot variability, and a first-order autocorrelation structure (AR-1) nested within plot to account for temporal correlation (Zuur et al., 2009). The resulting GAMM is:

where EN_ij is the ¹⁵N enrichment of sample j in plot i and EN_ij is assumed to follow a Gamma distribution with a log link. α is the intercept, f is a smoothing function (cubic regression spline) of the time effect and β is a parametric coefficient of the compartment effect. b_i is a normally distributed random intercept for the plot effect accounting for any dependency between samples from the same plot. The residuals at time s are modelled as a function of the residuals of time s−1 along with noise (as ε_s = ρ ε_s−1 + η_s, in which ρ is the temporal correlation parameter to be estimated; Zuur et al., 2009). The importance of variables for ¹⁵N enrichments was assessed by fitting various candidate models accounting for the main effects and their interaction (as f(Time_ij) x Comaprtment_ij) using maximum likelihood estimation (ML). Model selection was based on the Akaike Information Criterion (AIC) and the model with the smallest AIC and highest Akaike weights (w_i) was identified as the minimum-adequate model (Burnham and Anderson, 2002). Parameter estimates of the minimum-adequate model were based on the restricted maximum likelihood (REML) method (Zuur et al., 2009). Model assumptions were visually assessed following Zuur et al. (2010). We found no trends in the residuals. To avoid biased parameter estimates, we omitted extreme values (outliers, <2% of the data) for the analyses. Outliers were visually identified by Cleveland dotplots (Zuur et al., 2010). All analyses were conducted in R 3.3.1 (http://www.R-project.org) with the packages nlme (Pinheiro et al., 2016), mgcv (Wood, 2006), and MuMIn (Bartón, 2016).

δ¹⁵N values (natural abundance) varied between −8.0 and +8.7‰ across compartments and were lowest in Calluna biomass (range between −8.0 and −5.3‰). The moss layer showed δ¹⁵N values between −6.8 and −5.9‰. In the soils, δ¹⁵N values increased continuously with soil depth (organic layer: between −4.9 and −3.5‰, A-horizon: between 3.5 and 5.3‰ and B-horizon: between 6.1 and 8.7‰; Table 1).

As was to be expected, most compartments showed a significant (P < 0.05) increase in δ¹⁵N values (enrichment) after tracer addition (exceptions were the B-horizon, some measurements of the A-horizon and the litter samples in the first year of the experiment). The highest enrichment values were found for the moss layer and Calluna shoots (Table 1).

The direction and magnitude of temporal changes in ¹⁵N enrichment distinctly varied among compartments (Figure 1), as indicated by the strong support for interacting effects of time and compartment (Akaike model weight of 1.0; Table 2). 86 % of the variation in ¹⁵N enrichment was explained by the interaction between time and compartment, whereas compartment or time effects explained 45 and 7% of the variance, respectively (Table 2). The partial effects of the above-ground biomass of Calluna and the moss layer on ¹⁵N enrichment showed a significant decrease over time (P < 0.001 for both compartments; Table 3; Figure 1), while a significant increase of ¹⁵N enrichment was observed in the soil (O- and B-horizon: P < 0.001; A-horizon: P < 0.01; Table 3; Figure 1). ¹⁵N enrichment in the A-horizon also increased within the course of the experiment, but decreased in 2015 (Table 1).

Total ¹⁵N recovery was highest in the first year (about 93%); it then (almost linearly) decreased until 2010 to a value of about 54% (Figure 2; Table 4). Interestingly, recovery rates remained largely constant in subsequent years and achieved values between 49 and 54% until 2015. The recovery rates calculated increase by about 7% if the compartments which were not sampled annually are taken into account: “Calluna biomass older than 2 years” (about 2.5%), and “Calluna roots” (about 4.2%; means across sampling years; Table 4).

The moss layer proved to be the most important short-term sink and longer-term source for applied ¹⁵N, and showed a distinct decrease in recovery rates within the first 5 years (from 66% in 2007 to 14% in 2011; Figures 2, 3). This loss of ¹⁵N from the moss layer was not compensated for by a corresponding increase in recovery rates in any other compartment (Figure 3). We did find increasing recovery rates in the O- and A-horizons in the course of the experiment (e.g., from 18 to 40% for the O-horizon). This increase corresponded to about 50% of the total ¹⁵N losses from the moss layer between 2007 and 2015 (Figure 2). For the B-horizon, we found a slight but continuous increase in ¹⁵N recovery, with a maximum of 2.1% in the last year of the experiment.

Recovery rates in both the current year and 1- to 2-year-old shoots of Calluna were negligibly low in all study years, with values below 2%, and in most years even below 1% (Figures 2, 3, Table 4). Accordingly, we also found low recovery rates for the annual litter of Calluna, which amounted to ≤1%. In addition, leaching losses of ¹⁵N proved to be very low, as indicated by an annual recovery of 0.03% of the total ¹⁵N applied. Thus, leaching losses were even lower than at the beginning of the experiment (0.05% in 2007).

Our findings support our first hypothesis that the moss layer acts as both an important short-term sink and a longer-term source of captured ¹⁵N tracer. The moss layer thus represents a crucial link between airborne N inputs and the allocation of N to the heathland's N cycle (in term of its capture-release function). The moss layer's ability to store incoming N in the short term may be related to the fact that mosses represent ectohydric plants that absorb dissolved nutrients at the cell surface (Proctor, 2008). The decadal-scale availability and release of N could be attributable to a longer-term sequestration of N due to N assimilation by mosses and the microbial community associated with the moss phytomass (Gordon et al., 2001; Tye et al., 2005).

An important finding in terms of the N balance of the ecosystem was that, during the first 4 years of the experiment, ¹⁵N losses from the moss layer were not compensated for by a corresponding increase in recovery rates in any other compartment. These losses accounted for more than 30% of the total ¹⁵N applied and suggest that N losses from the moss layer in gaseous forms are very likely (Startsev and Lieffers, 2007; for a further discussion of this aspect see paragraph ‘Leaching losses and putative flows’ below). In contrast, the largely constant recovery rates found from 2010 onwards point to a conservative N cycling with marginal ¹⁵N losses from the N pool of the ecosystem in the long term (Tye et al., 2005).

Our second hypothesis was not confirmed, since the O-horizon (and not the B-horizon; see discussion below) proved to be an important long-term sink for sequestered ¹⁵N. In this horizon, recovery rates more than doubled in the course of the experiment. Several mechanisms may explain the N-flow to this compartment: Firstly, a passive downward transport of N from the moss to the organic layer via leaching might have taken place (Startsev et al., 2008), particularly at the beginning of the experiment. Secondly, decomposition processes of the lower and dead segments of the bryophytes, which sequestered ¹⁵N in biomass after tracer addition in 2007, are likely to have contributed to an increase in ¹⁵N recovery in the O-horizon over time (Aldous, 2002). Thirdly, processes such as mineralisation of Calluna litter and leaking bryophyte cells (as a result of desiccation-rehydration processes) could have contributed to increasing ¹⁵N recovery rates in the O-horizon (Bates, 2008; Friedrich et al., 2011b). Downward allocation processes of N (both as particular and dissolved organic matter) in the soil profile also explain an increasing recovery of ¹⁵N in the A-horizon. These findings coincide with observations in grasslands and forests, which showed the upper soil horizons (and the O-horizons in particular) to be important long-term sinks for airborne N loads (Gerzabek et al., 2004; Nadelhoffer et al., 2004; Templer et al., 2012; Gurmesa et al., 2016). High N retention in the O-horizon was also reported from heathland studies (although not based on tracer analyses; Kristensen and McCarty, 1999; Kristensen, 2001; Schmidt et al., 2004; Pilkington et al., 2005), and mainly attributed to an immobilisation of N by the microbial biomass present in these horizons (Jonasson et al., 1996; Nordin et al., 2004). A reverse trend in N recovery in the O- and A-horizons observed in the last 2 years of the experiment suggests that upward transports of N in the soil profile are also possible, probably due to bioturbation via root growth or soil invertebrates (Gabet et al., 2003).

Despite a continuous increase in ¹⁵N enrichment and recovery in the podzol's B-horizon, we found no clear support for our assumption that considerable quantities of N were sequestered in this horizon. We consider a decadal-scale recovery of about 2% too low to cause an ecologically relevant N withdrawal from the N cycle (which would be related to an N storage in insoluble and not readily bio-available forms; Hagedorn et al., 2005). Since the formation of a podzol spodic horizon is assumed to take centuries (Brady and Weil, 2008), an accumulation of N up to one ton per ha in this horizon (as found for our study area; Friedrich et al., 2011b) might require timespans far beyond the observation time of the present study. This assumption would be backed by the observed ¹⁵N accumulation rates, given that these could be extrapolated to a centennial scale.

Recovery of ¹⁵N in above- and below-ground biomass of Calluna was very low compared to recovery rates in the moss and organic layer. We hypothesise that soil microbes in the organic layer were the superior competitors for N (Andresen et al., 2008), and only small amounts of ¹⁵N remained available for plant uptake (Schimel and Bennett, 2004). This was demonstrated by Larsen et al. (2012) in a subarctic heathland, where 65% of the applied ¹⁵N was sequestered by soil microbes. In this context, it is important to note that high immobilisation rates by soil microbes are also indicative of a low N saturation level (Curtis et al., 2005).

¹⁵N losses via leaching proved to be negligibly low. This is an important indication of conservative N cycles (Tye et al., 2005) and contradicts the notion that the heathlands of the study area have achieved an advanced stage of N saturation (Aber et al., 1998), despite decades of N inputs beyond critical load thresholds. Obviously, the sites analysed still have a high capacity to sequester pollutant N (Choudhary et al., 2016), suggesting that a constant proportion of incoming N enters the heathland's N cycle and contributes to a slow but steady increase in the N level (De Schrijver et al., 2011). This confirms our third hypothesis.

In contrast to the observed accumulation of N, the first 4 years of the experiment (2007–2010) also exhibited a considerable loss of ¹⁵N. Tracer mass balances suggest that this loss was attributable to a leak of N from the moss layer in gaseous forms (Startsev and Lieffers, 2007). High N losses from bryophyte layers due to denitrification have been found for different species and ecosystems, and increase with increasing humidity (e.g., with water saturation of the moss layer after precipitation events), temperature, and N deposition (Laverman et al., 2000; Opelt and Berg, 2004; Lenhart et al., 2015). As a consequence, gaseous N losses from moss carpets are considered to have significant impacts on ecosystem N budgets (Fang et al., 2015), and also contribute to global nitrous oxide emissions from terrestrial ecosystems (Lenhart et al., 2015). For example, Calvo-Fernández et al. (2015) found high N losses as a result of denitrification in Cantabrian Mountain heathlands, when soils were water-saturated during winter. Fang et al. (2015) concluded that the effects of denitrification losses on ecosystem N balances are often underestimated. NH₃ volatilisation might also have contributed to ¹⁵N losses from the moss layer (Mahendrappa and Ogden, 1973; Startsev et al., 2008). Further mechanisms of N losses could be lateral N transport by mycorrhiza or soil fauna. However, since the natural ¹⁵N abundances in non-labelled subplots remained unchanged in the course of the experiment, we consider these losses (including conceivable lateral transports due to surface run-off) to be of minor importance. Low lateral transports also might be related to high absorption rates of ammonium, which is the prevailing N_inorg form at podzol sites with strongly acidic soil condition (cf. study site description). Due to the avoidance of destructive sampling, N contents of roots below a soil depth of 15 cm were not quantified. Since Calluna is considered a shallow-rooted plant (with more than 90% of root biomass being located in the humus horizons; Gimingham, 1972), and CN ratios of deep rooting woody roots are exceptionally high (>120; own unpublished data), we consider the N pool related to this compartment to be very low. We also rule out a significant effect of N deposition on ¹⁵N recovery patterns, since N deposition is evenly distributed across the study area, and thus between-plot differences in N deposition are low (including differences in N forms and isotope signatures; cf. deposition measurements of Niemeyer et al., 2005).

In conclusion, we found evidence that about 60% of the added ¹⁵N remained in the N cycle of the ecosystem at least for about one decade. In this context moss layers may play a crucial functional role for the N cycle of ecosystems: On the one hand, moss layers may contribute to high N losses by fostering denitrification processes (particularly under oceanic to sub-oceanic climates), and thus may decelerate N eutrophication processes in ecosystems. On the other hand, mosses represent a considerable short-term sink for incoming N which is released to the N cycle of an ecosystem in the long term. As a result, a considerable proportion of airborne N inputs can accumulate in different ecosystem compartments. Since the O-horizon proved to be an important long-term sink, it is likely that N captured by this horizon remains bio-available in the long term.

This has important implications for both critical load estimates and the development of management strategies to counteract airborne N inputs. A continuous accumulation of bio-available N (even under low N inputs) would implicate that critical load revisions need to account for cumulative effects of (low) N additions into ecosystems (Payne et al., 2011; Phoenix et al., 2012; Bähring et al., 2017). Mitigation strategies (in terms of management) should comprise high-intensity measures (e.g., choppering and sod-cutting; Härdtle et al., 2006; Frouz et al., 2009) as an important means to achieve balanced N budgets in the long term, because only these measures affect upper soil horizons with high contents of accumulated N.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.