Edinburgh Research Explorer Legacy effects of nitrogen and phosphorus additions on vegetation and carbon stocks of upland heaths

stocks did not respond to nutrient treatments. We found 40% of the added P had accumulated in the soil. (cid:1) This study showed persistent effects of N and N + P on vegetation composition, whereas effects of P alone were small and showed recovery. We found no indication that P application could mitigate the effects of N on vegetation or increase C sequestration in this system.


Introduction
Terrestrial ecosystems store vast amounts of carbon (C), estimated to range between 1500-2000 Pg C for vegetation and soil combined (Smith et al., 1993). The fate of this terrestrial biosphere C stock will determine the rates of global change via C-cycle climate feedbacks. The magnitude of the terrestrial C stock reflects the balance between C sequestration via net primary production and C losses, which are dominated by decomposition. In systems with low temperatures, wet conditions and litter of poor quality (high C : N), decomposition is slow, resulting in a buildup of C in the soil. Terrestrial nutrient availability, especially nitrogen (N) and phosphorus (P) strongly influences both vegetation productivity and decomposition processes and thus determines how these C stocks, and hence the global scale C-cycle, respond to environmental change (Zaehle, 2013;Wieder et al., 2015). Earth system models now include nutrient C-cycle interactions (Goll et al., 2017) but there are still substantial gaps in our knowledge, particularly around the interactive effect of N and P availability on long-term C-cycle responses.
Plant productivity is dependent on nutrient availability, especially the availability of N and P (Vitousek & Howarth, 1991;Vitousek et al., 2010). Anthropogenic pressures, such as atmospheric N deposition, can change nutrient availability and thus affect C sequestration . Elevated N deposition usually increases plant productivity (Phoenix et al., 2012) and in forests and heathlands, this has been shown to lead to increased C sequestration (de Vries et al., 2009). In nutrient-poor systems such as peatlands, alpine and arctic systems, however, increased C sequestration does not always occur under elevated N (de Vries et al., 2009;Street et al., 2017). The suggested mechanisms for this is that N deposition can cause vegetation to shift from N limitation to limitation by P or NP co-limitation, without increasing productivity, as observed in Calluna vulgaris-dominated upland heaths (Kirkham, 2001;Pilkington et al., 2005a), sand dune grassland (Ford et al., 2016) and mosses within an acid grassland (Arr oniz-Crespo et al., 2008). As a result of this shift in limitation, production rates saturate at lower rates of N deposition.
Carbon loss through decomposition can also be affected by nutrient availability (G€ usewell & Verhoeven, 2006), with impacts differing between systems. The main mechanisms by which nutrient availability could enhance decomposition are through reduction in litter C : N ratio (Bragazza et al., 2012;B€ ahring et al., 2017;Britton et al., 2018), changes in production Ó 2020 The Authors New Phytologist Ó 2020 New Phytologist Trust New Phytologist (2020) 1 www.newphytologist.com This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. of recalcitrant litter and changes in litter quality and quantity related to changes in vegetation composition. Conversely, reduced N limitation of microbial communities could decrease decomposition through the reduced need to mine recalcitrant organic matter for N (Craine et al., 2007;Milcu et al., 2011). In nutrient-poor systems, including the UK uplands, increased N deposition and experimental N addition increased mass loss from litter bags (Britton et al., 2018), increased CO 2 emissions (Bragazza et al., 2006;Kivim€ aki et al., 2013) and decreased the thickness of the litter layer (Armitage et al., 2012a;Street et al., 2017) indicating enhanced rates of decomposition. When an increase in decomposition exceeds the increase in productivity, depletion of ecosystem C stocks will result. By contrast, addition of ≥ 20 kg N ha À1 yr À1 in an upland C. vulgaris heath in Wales increased the litter layer and thus C storage (Field et al., 2017), and in a Scottish low-alpine heathland no effect of N addition on decomposition was found (Papanikolaou et al., 2010).
One potential reason for differing responses of plant productivity and decomposition to high N inputs across ecosystems is variation in P availability. Elevated N availability promotes upregulation of plant P acquisition, demonstrated by increased phosphatase activity in vascular (Johnson et al., 1999;Pilkington et al., 2005b) and nonvascular plants ( Arr oniz-Crespo et al., 2008;Phuyal et al., 2008) and in soil (Johnson et al., 1998). This can delay the shift towards P limitation in plants when N is applied (Marklein & Houlton, 2012), and hence explains why elevated N deposition does not always affect plant N : P ratio (Kirkham, 2001;Britton & Fisher, 2007), sometimes even showing a decrease in N : P ratio with increasing N deposition (Rowe et al., 2008). Increasing P acquisition may allow plants to respond to high N inputs by increasing productivity and thus maintain net C sequestration. Applying this mechanism, it has been suggested that addition of P could be used as a means of promoting C sequestration in ecosystems exposed to high N deposition (Armitage et al., 2012b). The importance of including P in C cycling studies is recognised (Goll et al., 2012), but the role of P is much less studied than the role of N. Furthermore, despite recognition of the importance of long-term experiments for understanding mechanisms of ecosystem change (Silvertown et al., 2010), there is a lack of long-term data on the effects of N and P on C stocks.
The UK uplands, defined as areas lying above the altitudinal limits of enclosed farmland (Ratcliffe & Thompson, 1988) are an example of ecosystems containing large C stocks. These areas experience a cool and wet climate and are dominated by heathlands, bogs and acid grasslands (Averis et al., 2004). Slow rates of decomposition relative to vegetation productivity promote the development and maintenance of large C stocks. The effects of N and P on these upland C stocks have received relatively little attention, however there have been some studies of impacts of combined N and P addition on productivity, both in the UK uplands and in analogous low nutrient ecosystems elsewhere.
Increased moss productivity is one mechanism by which combined N and P addition may increase C stocks in low nutrient ecosystems. Mosses allocate more C to recalcitrant tissue than vascular plants (Woodin et al., 2009), resulting in low decomposition rates (Lang et al., 2009). Increased moss productivity in response to combined N and P application has been demonstrated to increase the depth of organic soil horizons, and their C stocks, 20 yr after nutrient application in the Arctic (Street et al., 2015(Street et al., , 2017. Armitage et al (2012) showed an increase in moss productivity after P addition in a UK alpine heath with high N deposition. Similarly, moss growth was stimulated after combined N and P application in C. vulgaris heath (Pilkington et al., 2007) and in Vaccinium myrtillus heath and Festuca ovina grassland (Stiles et al., 2017). Sphagnum growth increased in response to 2-3 yr N and P treatment in a Scottish raised bog, (Carfrae et al., 2007), Irish raised bogs, Dutch fen bogs (Limpens et al., 2004) and a Patagonian boreal bog (Fritz et al., 2012). If mosses in UK upland plant communities respond similarly to combined N and P application, this mechanism could lead to increased C storage in upland habitats.
Increased vascular plant cover in response to N and P may also increase ecosystem C stocks. In an Arctic heath tundra, combined application of N and P was found to increase vascular plant cover, stimulating photosynthetic activity by 500% and ecosystem respiration by 250% (Arens et al., 2008), suggesting increased C sequestration. However, in UK upland heath, combined N and P addition has been observed to cause changes in vegetation composition, favouring grasses over heather (Hartley & Mitchell, 2005). This could result in reduced C storage, as grass-dominated heath has been found to sequester less CO 2 than C. vulgaris-dominated heath  and to have a lower soil C stock (Quin et al., 2014). Thus, although productivity could increase in response to P where N availability is high, a change in plant species composition may also occur, which could have positive or negative effects on C stocks.
These findings demonstrate the need to better understand the mechanisms by which P interacts with N to affect C sequestration and turnover, particularly over decadal timescales. Such knowledge is critical, not only to understand the basic ecology of upland systems and the response of C stocks to environmental change, but also to provide evidence on the use of P additions to manage C stores. The aim of this study was therefore to assess the long-term interactive effect of P and N on vegetation composition and C sequestration in an upland dry heath ecosystem. We sampled a 23-yr-old experiment on an upland heath where N, P and N + P treatments had been added for 14 yr. We described vegetation composition and estimated soil and vegetation C stocks, 10 yr after the last nutrient application. Given the suggestion that P application might be used to increase C sequestration (Armitage et al., 2012b), and the potential for such a management practice to affect freshwaters if P has any mobility within the system, we quantified the retention of applied P in soil and vegetation within the treatment plots.
We hypothesised that: (1) vegetation composition would have changed following combined N and P addition, favouring fast growing species (e.g. grasses) over the initial dwarf-shrub-dominated vegetation, with effects persisting after treatments ended; (2) changes in vegetation composition would affect vegetation and soil C stocks; (3) plots receiving N would have a smaller total C stock than controls, whereas plots receiving both N and P New Phytologist (2020) Ó 2020 The Authors New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com

Research
New Phytologist would have larger total C stocks as a result of the alleviation of N-induced P limitation of the vegetation; and (4) most of the originally added P would have accumulated in the soil because P is largely immobile.

Experimental setup
In 1993 a fertilisation experiment was set up at four upland sites in north-east Scotland, on dry heath vegetation dominated by C. vulgaris (UK National Vegetation Classification type H12 (Rodwell, 1991), EUNIS F4.2 (Strachan, 2017)). Two of the sites were situated in Glen Clunie and two in Glenshee, the sites within each valley being similar in soil properties and grazing pressure ( Table 1). The sites are at an elevation of c. 500 m above sea level (asl), with an average rainfall of 890 mm and a mean annual temperature of 6°C (Alonso et al., 2001). Total N deposition at these sites is currently c. 12.6 kg N ha À1 yr À1 (APIS, 2017). The vegetation is grazed by sheep and deer; in 1993 grazing levels were classified as high in Glenshee (60% of C. vulgaris shoots grazed) and moderately high in Glen Clunie (40-50% grazed) (Hartley & Mitchell, 2005) but have subsequently declined (R. J. Mitchell, pers. comm.). Soils in Glen Clunie are higher in organic matter and soil moisture, and more acidic than soils in Glenshee, where a more base-rich subsoil is present (Hartley & Mitchell, 2005).
At each of the four sites, four 3 9 5 m plots were set up receiving either no fertiliser (control), N, P or both N and P (N + P). The original experimental design also included grazing and K addition treatments with 16 treatment combinations replicated at each of the four sites giving 64 plots (Hartley & Mitchell, 2005). For our study we resampled only the grazed plots receiving either no fertiliser, N, P, or N + P. Each of these treatments was replicated once at each of the four sites giving 16 plots across all sites. Fertiliser was applied twice per year from April 1993 until 2006 (14 yr). Nitrogen was applied at a rate of 7.5 g N m À2 yr À1 (i.e. 75 kg N ha À1 yr À1 ) in the form of NH 4 NO 3 and P was applied at a rate of 5 g P m À2 yr À1 (i.e. 50 kg P ha À1 yr À1 ) in the form of triple superphosphate (Ca (H 2 PO 4 ) 2 ÁH 2 O), both in solid form. At 10 yr after the last nutrient application, we revisited the plots and sampled vegetation composition and nutrient and carbon stocks.

Vegetation and soil sampling
In 1993 and 1996 cover of C. vulgaris and grasses (all species combined) was estimated in three 1 9 1 m quadrats per plot using point quadrats (Hartley & Mitchell, 2005). Additionally, in 1996 and 1999 full vegetation composition was recorded in these quadrats, and in 2000 one of the three quadrats was harvested for biomass estimation. No data had been collected on full initial vegetation composition or at the time of cessation of the nutrient treatments (2006). In July 2016 we recorded vegetation composition in the two remaining quadrats used in the previous vegetation surveys, and in a third quadrat adjacent to the quadrat harvested in 2000. Vegetation composition was noted by recording species identity and estimating cover to the nearest 5% (or to the nearest 1% if below 10%). We then calculated coverweighted Ellenberg fertility (N) and pH (R) values per quadrat, based on values in PLANTATT and BRYOATT (Hill et al., 2004(Hill et al., , 2007. Shannon species diversity (Shannon H) was also calculated per quadrat.
In August 2016, to quantify aboveground biomass to the soil surface (including moss and litter) and nutrient stocks in the vegetation, we harvested a 50 9 50 cm area in the middle of three different, randomly selected, 1 9 1 m quadrats per plot, excluding quadrats already used in previous analyses. Samples were kept frozen until sorting into five groups: shrubs, bryophytes and lichens, graminoids (including grasses, sedges and rushes), herbs and litter. Oven-dried biomass (48 h at 80°C) was determined per group per quadrat. Material was ground (Retsch cross beater mill, 2 mm sieve size) and subsampled for P analysis. For C and N analysis a subsample was steel ball milled (Retsch MM400; Retsch, Germany).
We sampled soil from the vegetation harvest area by taking three soil cores per quadrat (diameter 5 cm, 15 cm deep). These cores were split into organic and mineral soil horizons, dried (7 d, 30°C), sieved (2 mm sieve) and weighed to determine bulk density. Large roots that were removed by the sieve were weighed, ground (cross beater mill, as for vegetation) and recombined with the soil sample. Soil samples were then bulked per horizon per quadrat and ball milled before nutrient analysis.
Carbon and N contents of soil and vegetation were determined on the milled material using elemental analysis (NA 2500 Series 2; CE Instruments Ltd, Wigan, UK). Total P, including all chemical forms of P, was determined using flow injection analysis (FIAstar 5000; Foss Analystica AB, Hilleroed, Denmark), following sulphuric acid-hydrogen peroxide digestion (Grimshaw, 1987). Soil pH was determined using a 1 : 2.5 ratio of dried soil and deionised water. Phosphorus and C stocks were then calculated per m 2 , including the top 15 cm of the soil, aboveground vegetation and litter.
In September 2017, we collected three additional soil samples per plot for microbial P analysis (diameter 5 cm, 15 cm deep). Fresh soils were sieved on a 4 mm sieve, removing roots and stones. Soils were kept at 4°C until analysis (undertaken within 8 wk) for microbial P following Brookes et al. (1982) using chloroform fumigation, extraction by 0.5 M NaHCO 3 and colorimetric analysis for phosphate (Murphy & Riley, 1962). As some soils were very organic and yielded a coloured extract, 1 g of activated charcoal was added to all extracts to decolour them before phosphate analysis. The estimate for microbial P was then calculated using the difference between the fumigated and nonfumigated sample, correcting for the recovery efficiency in a third, nonfumigated sample spiked with 250 µg of PO 4 -P (Brookes et al., 1982;Voroney et al., 2007). One of the microbial biomass P calculations yielded a negative result and therefore this sample was excluded from the microbial biomass analysis.

Data analysis
Soil and vegetation P and C stocks, P concentrations, microbial biomass P, pH, total vegetation biomass and biomass N : P ratios were analysed using linear mixed effect models in R v.3.4.1, package NLME (Pinheiro et al., 2017;R Core Team, 2017), with treatment N, P and the interaction between N and P (N 9 P) as factors. Change in cover of C. vulgaris and grasses relative to the cover at the start of the experiment was analysed using repeated measures ANOVA over four time points (1993, 1996, 1999 and 2016) with N addition, P addition and year as factors. Similarly, species richness (number of species per quadrat) and diversity were analysed for three time points: 1996, 1999 and 2016 using N addition, P addition and year as factors. Due to the original spatial layout of the plots, within each site, N and P plots, and control and N + P plots were always located together in pairs; this pairing was accounted for in models (defined as 'block'). In all models, a nested error term with the plot, block, site and valley (1|valley/site/block/plot) was included. Additionally, when analysing the changes over time, plot identity was included as a random factor. Normal distribution of the residuals was checked for all models, and response variables were sqrt, inverse sqrt or log e transformed when necessary to meet normality assumptions. For soil C stock and litter C : N ratio, nonparametric Kruskal-Wallis tests were used to test treatment effects, as normality assumptions could not be met.
For analysis of community composition, species per cent cover data were averaged per plot and analysed using multivariate statistical methods in CANOCO 5 (ter Braak & Smilauer, 2012) in order to examine differences in species composition per treatment and the change in species composition over time. Species with a cover < 1% were set to 0.01% cover and species with occurrence in less than three plots were removed from the analysis. Species composition change over time was analysed using principal response curves (PRC). In this analysis, species composition in the control treatment was set as the baseline for all the years, and all fertiliser treatments were compared with this baseline. Significance of the resulting ordination axis was tested using Monte Carlo permutation tests and individual species scores were calculated to show which species changed most in abundance along this ordination axis. To study the resulting difference in species composition between treatments in 2016 we used a partial redundancy analysis (RDA), so that we could constrain the ordination by differences between the valleys.

Vegetation composition
From 1996 to 2016 species richness increased across all treatments, but there were no significant effects of N or P addition on species richness (Fig. 1a; Table 2). In 1999, Shannon H increased in all treatments except for the N + P treatment (Supporting  Information Table S1), but through time the treatment effects converged (significant interaction between N, P and year; Fig. 1b; Table 2). Consequently, species richness and Shannon H in 2016 were not significantly different between treatments (Table S1). In 2016, an average of 18 plant species were found per plot (Fig. 1a), with the most common species being the shrubs C. vulgaris and V. myrtillus, the graminoids Agrostis vinealis and Nardus stricta, the herbs Potentilla erecta and Galium saxatile and the mosses Hylocomium splendens, Pleurozium schreberi and Rhytidiadelphus squarrosus. Cover-weighted Ellenberg N, an indicator of fertility, increased in response to N addition and P addition ( Fig. 1c; Table 2). Ellenberg R, an indicator of soil pH, increased through time in all treatments, indicating a decrease in acid-loving species, which was reinforced by both N addition and P addition ( Fig. 1d; Table 2). The main change in Ellenberg N and R occurred between 1993 and 1999; by 2016 Ellenberg N and R no longer differed significantly between treatments (Tables 2, S1).
Calluna vulgaris cover decreased through time in all treatments, from an average of 66% in 1993 to 28% in 2016. This decrease in C. vulgaris cover was enhanced by addition of N alone and N + P ( Fig. 2a; both N and P were significant factors; Table S2), with the lowest C. vulgaris cover in the N + P treatment. Grass cover was higher with both N and N + P addition; this was consistent in both 1996 and1999 (Fig. 2b; both N and P were significant factors; Tables S1, S2). In 2016, the resulting grass cover was higher in all plots that had received N addition, whereas C. vulgaris cover was lower in these plots ( Fig. 2; Table S1). Phosphorus addition alone had no effect on C. vulgaris or grass cover in 2016, although P was a marginally significant factor in the models, reflecting the large cover changes in N + P treatments ( Fig. 2; Table S1). Moss cover was significantly lower under N addition and increased over the years in all plots (Table 2).
PRC analysis also demonstrated that the effects of the treatments on total vegetation community composition changed over the years; the interaction between year and treatment explained 27% of the total variation in species composition (Table S3). After correcting for year, treatments explained 66% of the variation in species composition on the first ordination axis (PRC.1, Fig. 3; Monte Carlo significance test, pseudo-F = 9.6, P = 0.026; Table S3). Abundance of the shrubs C. vulgaris and V. myrtillus and the mosses H. splendens, P. schreberi and Sphagnum spp. were negatively related to nutrient treatments, whereas abundance of the grasses A. vinealis, F. ovina, N. stricta and the mosses New Phytologist (2020) Ó 2020 The Authors New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com

New Phytologist
Polytrichum commune and R. squarrosus were positively related to nutrient treatments (Fig. 3). Change in vegetation composition compared with the control was greatest in the N + P treatment (higher score on the PRC axis) while the P only plots became more similar to the control treatment over the years.
The strongest influence on species composition in 2016 was the difference between sites (i.e. Glen Clunie and Glenshee). After constraining the ordination for site differences, nutrient treatment explained 36% of the remaining variation in plant species composition (Fig. 4). Ordination of the 2016 data showed a small difference in plant species composition between the control and P only plots; the mosses H. splendens and P. schreberi being more abundant in the P treated plots. Plots receiving N only and N + P were more different from the control; while control plots were strongly associated with C. vulgaris, in N addition plots sedges (Carex sp.) and the herbs Viola riviniana and Polygala serpyllifolia were more common, and in N + P addition plots the shrub Vaccinium vitis-idaea, the mosses Values are means, error bars show AE 1 SE, n = 12, see Table 2 for statistical tests. Numbers in bold indicate significant differences, n = 12. Herb cover and bryophyte cover were sqrt transformed before analysis.  Table S3, for statistical tests.

Carbon stock
Soil carbon stock in the top 15 cm of the soil was on average 9.4 kg C m À2 and was not significantly affected by treatments (Kruskal-Wallis v 2 = 1.91, df = 3, P = 0.589; Fig. 5a). Vegetation C stock was an order of magnitude smaller, on average 0.9 kg C m À2 , and tended to be lower under N addition, although not significantly so ( Fig. 5b; Table S1). This marginal difference could be attributed to lower shrub biomass in the combined N and P addition plots (173 g m À2 compared with 647 g m À2 in the control plots), making up a large part of the vegetation C stocks (Fig. 5b). Litter fractions contributed most to vegetation C stocks, followed by shrubs and bryophytes; graminoid and herb C stock were much smaller. Carbon stocks in each of the different vegetation fractions did not differ significantly between treatments, except for higher herb biomass and associated C stock with N addition (Table 3; Table S2). There was no treatment effect on litter carbon pools, but litter C : N ratio decreased with N addition (Table 3). There was no significant effect of N or P addition on organic layer depth. As soil C stock makes up the majority of the total C stock, there were no effects of N or P addition on total C stock.

Fate of the added phosphorus
In total, 70 g P m À2 was added during the 14-yr experimental period. Average total P stock in the top 15 cm of soil ranged from 58 g P m À2 in the control treatment to 94 g P m À2 in the P only treatment and was significantly greater where P had been added (mean increase 28 g P m À2 compared with plots without P addition, Tables 4, S1). Microbial P stock was on average 14.4 g P m À2 , which is 19% of the total soil P. Microbial P stock showed a similar pattern to total soil P stock, being on average 9.2 g P m À2 higher in P treated plots compared with plots not receiving P (Table 4). Total vegetation P stock averaged 1.9 g P m À2 , only 2.5% the size of the soil P stock, and was higher under P addition (+ 0.6 g P m À2 ) and lower under N addition (À0.4 g P m À2 ), with no interactive effect of N and P (Table 4). The differences Fig. 3 Principal response curve for species composition over time, showing differences in species composition for the first ordination axis (PRC.1, P = 0.026) between nutrient addition treatments. The arrow indicates the last nutrient addition. C, no fertiliser, N, N only addition, P, P only addition and NP, addition of both N and P combined. Vegetation composition in the control plot is set at 0 along the PRC axis. The species score (right panel) depicts the direction and the relative magnitude of species occurrence along the PRC axis compared to the control. Multiplying the exponent of the species score with the value on the PRC axis for a specific treatment gives the approximate difference in abundance for that species compared to the control. A negative species score indicates a reduction in occurrence of that species along the PRC axis, whereas a positive species score means an increase compared to the control.

Research
New Phytologist in vegetation P stocks were not due to differences in biomass but were mainly related to increases in tissue P concentration in response to P and decreases in P concentration in response to N in all functional groups (Tables 3, 4, S4). Tissue N : P ratio in the vegetation was lower under P addition, due to increased P concentration (Table 3).

Discussion
Our findings demonstrate that heathland plant communities have a limited capacity for recovery at decadal timescales following sustained (14 yr) mineral nutrient inputs, and that the community may have surpassed a critical ecological tipping point. Despite these shifts in vegetation composition, we found that heathland ecosystem C stocks were resilient to changes in nutrient inputs, suggesting plant communities and soil processes have limited connectivity, or respond at markedly different timescales, as has been reported in other system using observational analyses (e.g. Stark et al., 2019).
In line with our first hypothesis, vegetation composition initially changed most in response to N addition, and additionally to combined N + P addition, favouring grasses over C. vulgaris (Hartley & Mitchell, 2005). The shrub V. myrtillus showed a similar pattern to C. vulgaris, decreasing with N and P addition. Similar effects have been reported previously (e.g. Heil & Diemont, 1983) where C. vulgaris was replaced by F. ovina after 12 yr of N addition, but did not change under P addition. In pot experiments, N and P additions increased the competitiveness of N. stricta over C. vulgaris (Hartley & Amos, 1999). Other work suggests that the pioneer phase of heathland vegetation may be the sensitive phase when dominance by grasses is able to occur (Friedrich et al., 2011). However, our data show that over 23 yr, both N addition alone and, to an even greater extent, N + P addition increased the cover of grasses and reduced the cover of C. vulgaris, despite the vegetation being in a mature phase of growth. The increase in grass cover in response to nutrient addition was exacerbated by an overall trend for C. vulgaris cover to decrease in all plots, including the controls. This finding may be (a) (b) Fig. 5 Average (a) soil and (b) vegetation C stock per nutrient addition treatment. C, no fertiliser, N, N only addition, P, P only addition and NP, addition of both N and P combined. Soil C stock is split in mineral (grey) and organic (black) layer. Vegetation C stock is split for each fraction. Values are means, error bars show AE 1 SE of total soil C stock and total vegetation C stock, n = 12, ns = not significant, see Supporting Information Table S1 for statistical tests. C, no fertiliser, N, N only addition, P, P only addition and NP, addition of both N and P combined. Significance values: *, P < 0.05; **, P < 0.01; and ***, P < 0.001. Shrub, bryophyte and herb [P], graminoid and herb biomass were square root transformed, bryophyte biomass was loge transformed and litter [P] and C : N were inverse square root transformed before analysis. For the litter biomass a nonparametric Kruskal-Wallis test was used. Values are mean AE 1 SE, n = 12.
Ó 2020 The Authors New Phytologist Ó 2020 New Phytologist Trust New Phytologist (2020) www.newphytologist.com related to the relatively high grazing pressure, as exclosures in the same area showed an increase in C. vulgaris cover from 1993-1999 (Hartley & Mitchell, 2005). Grazing has previously been found to shift vegetation from C. vulgaris to grass domination (Bokdam, 2001) and grazing in combination with N deposition has been shown to decrease soil and vegetation C stock in upland areas across the UK . A change from C. vulgaris-dominated to grass-dominated heath can change soil C stocks, decreasing the amount of recalcitrant C (Quin et al., 2014).
The effects of N and combined N + P addition on vegetation composition were persistent over time. However, vegetation composition recovered in plots that had received P only, becoming more similar to control over time. This could be a legacy of N addition. Nitrogen stock data suggest that N has been largely retained in the soil, but not in the vegetation (Fig. S1). It is likely that the shift in vegetation composition in the N and N + P plots is being maintained by competition between grasses and shrubs, the system having passed a tipping point.
Mosses remain particularly sensitive to N inputs; their reduced cover reflects responses to N inputs seen in other systems such as Arctic and alpine heath (Armitage et al., 2012a;Street et al., 2015), and both H. splendens and P. schreberi have been negatively associated with high N deposition in the UK (Field et al., 2014). Increased moss biomass in response to P addition with high N (Pilkington et al., 2007;Armitage et al., 2012b;Stiles et al., 2017) suggests that P addition may be able to alleviate some of the negative impacts of N deposition on bryophyte cover and biomass and thus on their C storage. Although we observed some increases in nitrophilic mosses with combined N + P addition, we did not see any significant evidence for P alleviating negative effects of N on total moss biomass and C stock.
Shifts in vegetation composition did not translate into changes in vegetation C stocks, contrary to our second hypothesis. The decrease in C. vulgaris cover and biomass in the N + P treatment tended to decrease the vegetation C stock (Fig. 5b), and increases in grass cover only slightly offset this. Herb biomass and C stock was slightly higher with N addition, but this fraction makes up a very small proportion of total vegetation C stock (Fig. 5b). The lack of differences in biomass with changing plant species composition could indicate that functional plant species composition remained diverse enough to maintain productivity (Tilman et al., 1997). We measured standing biomass, rather than productivity, and it is possible that productivity may actually have been increased by the more rapid turnover of graminoid material, without changing standing biomass.
Contrary to our third hypothesis, and despite the changes in vegetation composition, soil C stocks and the depth of the organic layer did not respond to nutrient treatments. The most likely mechanism behind this lack of response is equal change in both of the components determining C stock, that is production and decomposition, or no change in either. Productivity has been found to generally increase with N deposition (Phoenix et al., 2012) and, similarly, C outputs often respond positively to N deposition (Bragazza et al., 2012;Kivim€ aki et al., 2013;Britton et al., 2018), although not in all studies (Papanikolaou et al., 2010;Field et al., 2017). Vegetation in our study shifted towards more productive, fast growing species (e.g. grasses) in response to nutrient addition, especially to N + P treatment. Fast growing species generally decompose faster; we observed a decrease in litter C : N under N addition, indicating more readily decomposable material. Thus, both productivity and decomposition may have increased, resulting in no net changes in C stocks. Our findings contrast with others that have shown both positive and negative responses of C stocks to N and P additions Field et al., 2017;Stiles et al., 2017;Street et al., 2017), suggesting that these responses are context dependent. While some ecosystem models currently assume positive effects of N deposition on C stocks (e.g. Tipping et al., 2017), our findings imply that this may not be ubiquitous.
Our findings highlight the continuing need to better understand the mechanisms regulating P storage and turnover in ecosystems, particularly given the scarcity of rock phosphate (Scholz et al., 2013) and the linkage between global atmospheric N deposition and P demand by plants (Johnson et al., 1999) and soil microorganisms (Johnson et al., 1998). Confirming our fourth hypothesis, the addition of P increased the P stock in both vegetation and soil, reflecting the limited mobility of this element in the soil (e.g. McGill & Cole, 1981). Few studies of N and P addition have reported data on soil P stocks, although, over the short term, added P has been found to be retained in vegetation in alpine Racomitrium heath (Armitage et al., 2012), and in Table 4 Soil, microbial and vegetation total phosphorus, vegetation biomass, depth of the organic (O) layer, pH of the organic and mineral (M) layer by treatment for 2016, including the significance of the ANOVA model factors nitrogen (N), phosphorus (P) and the interaction between N and P (N 9 P).

Treatments
Significance of factors

New Phytologist
Sphagnum bogs (Fritz et al., 2012). In our study, the difference in soil P stocks between plots with and without P addition accounted for 40% of the total added P (28.0/70 g P m À2 added), with a further 1% found in vegetation P stocks. This means that 59% of the added P was unaccounted for.
Despite the small contribution of vegetation to total P stocks, an important finding was that high N inputs reduced vegetation P stocks. This could potentially be explained by reduced availability of P to vegetation as a result of N-induced soil acidification, or by increased microbial demand for, and immobilisation of, P when N is added (Pilkington et al., 2005). However, we did not find any significant change in either soil pH in the organic layer or soil microbial biomass P in N fertilised plots, suggesting that neither of these explanations is correct. An alternative explanation for the reduced vegetation P stock could be that N addition reduces mycorrhizal colonisation (Yesmin et al., 1996), thus reducing the ability of plants to access soil P.
The apparent loss of 59% of added P from the system could be explained by two potential mechanisms, namely biomass removal by grazing and leaching of P beyond the 15 cm soil depth sampled. Herbivore browsing in this system was estimated to affect 40-70% of current years' shoots (Hartley & Mitchell, 2005). If we make a conservative assumption (given that sheep numbers have decreased since the start of the experiment) that this equates to 20% of biomass, it would account for 10.1 g P m À2 loss over 23 yr. Thus, grazing could only account for a small proportion of the unaccounted P, and it is more likely that P was retained at soil depths >15 cm. This would be expected to be due to the immobile nature of P and because average soil depths at both sites are 75-80 cm (Soil Information for Scottish Soils, SIFSS). Tight retention of P is supported by the fact that in a preliminary study (data not shown) we found no evidence of elevated P in the top 15 cm of the soil just 1.5 m beyond the edge of P treatment plots, suggesting limited lateral movement of P.
The amount of P not accounted for in the top 15 cm of soil in our study may reflect treatment application rates. Although high, the rates of N (75 kg ha À1 yr À1 ) and P (50 kg ha À1 yr À1 ) application were comparable with the highest application rates in other studies on bogs (20-64 and 4-64 kg P ha À1 yr À1 ) (Malmer et al., 2003;Carfrae et al., 2007;Phuyal et al., 2008;Currey et al., 2010;Larmola et al., 2013) and on tundra systems (50-100 and 5-50 kg P ha À1 yr À1 ) (Mack et al., 2004;Nowinski et al., 2008;Sundqvist et al., 2014;Street et al., 2015). However, studies in UK ecosystems similar to ours have utilised lower rates of addition, of 10-20 and 5-20 kg P ha À1 yr À1 (Pilkington et al., 2007;Armitage et al., 2012b;Stiles et al., 2017). Despite more than half the P added to our system not being retained in the top 15 cm of soil, 280 kg P ha À1 was retained. This suggests that if P was added at lower rates, even with repeated additions over several years, it would all be retained over decadal time scales, as long as the maximum retention capacity of 280 kg P ha À1 was not exceeded.
Our findings provide valuable insight into the legacy effects of nutrient inputs and limitation in seminatural ecosystems. How ecosystems respond to P limitation is a fundamental question that remains unclear in many ecosystems, and which is only beginning to be considered in land surface models (Reed et al., 2015). Our findings contribute to recent calls for better understanding of how P limits plant productivity and soil C, N and P pools (Jiang et al., 2019), and show that P limitation does not control C accumulation in organic soils. The results also suggest that addition of P to upland vegetation cannot mitigate the effects of N or increase C sequestration. Our study illustrates how shifts in biodiversity (e.g. substantial reduction in cover of the dominant shrub) do not necessarily translate to belowground function even over long timescales, and provide an empirical demonstration of the uncertainty in linking biodiversity change with ecosystem ecology (Adair et al., 2018).

Supporting Information
Additional Supporting Information may be found online in the Supporting Information section at the end of the article. Table S1 ANOVA table for soil and vegetation data in 2016 and for vegetation data of 1996 and 1999 testing the factors nitrogen (N) addition, phosphorus (P) addition and the interaction between N and P addition (N 9 P).
Ó 2020 The Authors New Phytologist Ó 2020 New Phytologist Trust New Phytologist (2020) www.newphytologist.com Table S2 ANOVA table for plant functional group biomass and tissue chemistry testing the factors nitrogen (N) addition, phosphorus (P) addition and the interaction between N and P addition (N 9 P). Table S3 ANOVA table for the changes in C. vulgaris and grass cover over time testing the factors nitrogen (N) addition, phosphorus (P) addition, year and the interaction between these factors.

Table S4
Results of the principal response curve (PRC) analysis with partial redundancy analysis (RDA), showing explained variation by the first four ordination axes.
Please note: Wiley-Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.
New Phytologist is an electronic (online-only) journal owned by the New Phytologist Trust, a not-for-profit organization dedicated to the promotion of plant science, facilitating projects from symposia to free access for our Tansley reviews and Tansley insights.
Regular papers, Letters, Research reviews, Rapid reports and both Modelling/Theory and Methods papers are encouraged. We are committed to rapid processing, from online submission through to publication 'as ready' via Early View -our average time to decision is <26 days. There are no page or colour charges and a PDF version will be provided for each article.
The journal is available online at Wiley Online Library. Visit www.newphytologist.com to search the articles and register for table of contents email alerts.
If you have any questions, do get in touch with Central Office (np-centraloffice@lancaster.ac.uk) or, if it is more convenient, our USA Office (np-usaoffice@lancaster.ac.uk) For submission instructions, subscription and all the latest information visit www.newphytologist.com New Phytologist (2020) Ó 2020 The Authors New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com