Introduction

Nitrogen enrichment is considered one of the most prominent threats to biodiversity worldwide (Sala et al. 2000). In the terrestrial environment nitrogen is an important limiting nutrient for primary producers over a wide variety of ecosystems (LeBauer and Treseder 2008). Worldwide, but most prominently in Western Europe, atmospheric deposition of reactive nitrogen originating from agriculture and combustion processes has led to a reduction of botanical diversity (Bobbink and Hettelingh 2011). The most popular explanation is that addition of nitrogen results in an increase of biomass production and a shift from nutrient limitation to light limitation and thus to the dominance of a few species that most efficiently convert light into biomass (Stevens et al. 2004; De Schrijver et al. 2011). Considerable effort is put into a reduction of nitrogen emissions (both as ammonia and as nitrogen oxides) (Erisman et al. 2003), and in local mitigation of the effects of excess nitrogen (Stevens et al. 2011). The ‘critical load’ approach (Bobbink and Hettelingh 2011) provides a widely accepted, quantitative basis to the nitrogen emission reduction policy.

The production increase and loss of diversity will only occur if nitrogen is indeed limiting. Locally, or even on a larger scale nutrients other than nitrogen may be limiting (Wassen et al. 2005; Elser et al. 2007). Of these, phosphorus is the most probable candidate although potassium (Van Duren and Pegtel 2000) and iron (Jickells et al. 2005) have also been proposed. In contrast to nitrogen or potassium, phosphorus is highly immobile in soils and plants have different strategies to actively acquire phosphorus (Vance et al. 2003; Olde Venterink 2011), for example through exudation of phosphatase (Kroehler and Linkins 1988) or mycorrhizal infection (Furlan and Bernier-Cardou 1989). This may lead to higher species richness under P-limited than under otherwise comparable N-limited conditions. Fujita et al. (2014) found that P-limited conditions favour species that invest little in sexual reproduction and are thus vulnerable to local extinction. Olde Venterink et al. (2003) found indications that under P-limited conditions endangered plant species persist at higher production levels than under N-limited conditions. P enrichment may therefore be even more detrimental to plant species diversity than N enrichment (Wassen et al. 2005). Ceulemans et al. (2014), in a large-scale gradient study, showed a clear negative relationship between species number in grassland plots and P availability, independent of the level of atmospheric N deposition. Like the global carbon (Falkowski et al. 2000) and nitrogen (Galloway and Cowling 2002) cycles, the phosphorus cycle has been accelerated by human activities during previous centuries, and probably even more strongly so (Wassen et al. 2005), which would be a serious threat to biodiversity if P-limitation is a widespread phenomenon.

In a large-scale meta-analysis of N and P fertilisation experiments, Elser et al. (2007) showed that over a wide range of marine, freshwater and terrestrial ecosystems, production is most often co-limited by nitrogen and phosphorus (this does not exclude co-limitation by other nutrients but these were not included in their study). These authors show that the occurrence of N + P co-limitation is not only common in the three biomes, but within the terrestrial environment it is also common over a wide range of ecosystem types (forest, grassland, wetland, tundra). Harpole et al. (2011) analysed a subset of Elser et al.'s (2007) database including only factorial addition experiments. When the effect threshold is set to a production increase of 35 % compared to blank treatments they report a form of co-limitation in c. 50 %, single-nutrient (either N or P) limitation in c. 15 %, and no or a negative effect in the remainder of their 641 investigated studies. Although co-limitation is apparently the most common situation, its form may vary widely, from purely additive to multiplicative or ‘serial’ (i.e. one nutrient has an effect and the second only if the first is also added). These results again seem to be surprisingly similar for the three biome types (terrestrial, freshwater, marine).

In a meta-analysis of nitrogen addition experiments LeBauer and Treseder (2008) arrived at similar conclusions for nitrogen limitation in terrestrial ecosystems. They also found that the ‘geophysical hypothesis’ (Walker and Syers 1976), which predicts a poleward increase of nitrogen limitation because glaciation increases phosphorus availability, is not generally true. Although co-limitation is apparently the common situation, single-nutrient limitation may locally occur for nitrogen and phosphorus but possibly also for potassium (Van der Woude et al. 1994). If nitrogen is not a limiting nutrient in an area with a high nitrogen deposition reduction of deposition may be ineffective and efforts should be directed towards reduction of inputs of phosphorus or potassium if the aim is to decrease biomass production in order to protect threatened species.

There is a question related to the effect on species composition of adding limiting nutrients. By definition the addition of a limiting nutrient will result in an increase in production. Usually, this increase is accompanied by a shift in species composition, which most often entails a reduction of the overall species number, but most notably of the number of threatened species (Wamelink et al. 2003; Olde Venterink et al. 2003; Hecjman et al. 2010). Studies on the effect of nutrient additions usually concentrate on species numbers rather than on species composition (cf. Grace 1999; Mittelbach et al. 2001), and often have the aim of testing Grime’s (1979) ‘hump-back hypothesis’ (cf. Grace et al. 2016 for a critical re-assessment of this hypothesis). In considerations of species composition two hypotheses may be put forward: (a) addition of a limiting nutrient will lead to a production increase, and consequently to a shift towards more productive species (of which there are less than of low-productive species); or (b) addition of a limiting nutrient will lead to a shift in the relative availabilities of nutrients, and thus to a shift in species composition, as species differ in their abilities to exploit various resources (the ‘resource-ratio hypothesis’, Tilman 1982). In the former the species composition would only depend on the production, while in the latter a different limiting nutrient would lead to a shift in species composition, even if production remains unchanged.

Here we present the results of a long-term factorial fertilisation experiment in grassland on drained peat, where the following treatments were applied: blank, addition of potassium, addition of potassium and phosphorus, and of potassium, phosphorus and nitrogen. This experiment was originally set up in order to find optimal fertiliser applications for the agricultural use of drained peat. It attracted our attention because: (a) a part of the plots is rich in species compared to normal agricultural grassland; and (b) the experiment exhibits a strong co-limitation of P and K, but not of P and N (Kamínsky and Szymanowski 2007). Our aims are: (1) to investigate the effect of the fertilisation treatments on above-ground biomass production as a test of Harpole et al. (2011) hypothesis of co-limitation as the most common situation (as opposed to Von Liebig’s (1840) Law of the Minimum that predicts single-nutrient limitation); (2) to investigate the effect of the treatments on the overall number of species and presence of endangered species as a test of Olde Venterink et al.'s (2003) hypothesis of endangered species persistence under P-limited but N-rich conditions; and (3) to investigate whether, in a production gradient, species composition only depends on production or also on the nature of the limiting nutrient(s) as a test of Tilman’s (1982) resource-ratio hypothesis. An important limitation of our experiment is that it is not full-factorial, i.e. it does not include application of phosphorus or nitrogen alone. This means that, except for potassium, we cannot determine the effect of adding a single nutrient, and, depending on the outcome of the experiment, it will not be possible to discriminate between co-limitation and single-nutrient limitation in all cases.

Materials and methods

The long-term agro-technical experiment at the Research Station for Land Reclamation and Grassland Farming (IMUZ) Biebrza (Fig. 1) was installed in the period 1957–1959 (Gotkiewicz et al. 1980; Gotkiewicz 1987; Okruszko 1989). Its original aim was to determine optimal fertilisation treatments for increasing agricultural productivity on drained peat in Eastern Poland. The experiment is located at Biebrza, near Grajewo, Poland (53°38′27″ N, 22°35′17″ E), and consists of 216 adjacent plots of 14 × 6 m that form a factorial management x fertilisation experiment. The experiment is fenced as a protection against large herbivores.

Fig. 1
figure 1

Overview of the experiment during harvest 1st of July 1996

The management treatments are: (1) permanent grassland; (2) arable land; and (3) ‘ley’ (i.e. alteration of grassland and arable land); however, the present analysis confines itself to the permanent grassland plots. The fertilisation treatments are: (1) blank; (2) addition of potassium; (3) addition of potassium and phosphorus; and (4) addition of potassium, phosphorus and nitrogen. These treatments are denoted as O, K, PK and NPK, respectively. These nutrients were added yearly, their forms and quantities are given in Table 1. In the permanent grassland plots some species were sown at the start of the experiment (Appendix I).

Table 1 Treatments

There are three blocks, but in 1992 scrub encroachment began to take place when treatment of one of them was discontinued. The present analysis therefore confines itself to the two remaining blocks. Each block contains 6 replicates per treatment, placed in a randomised design; the total number of plots analysed is therefore 48 (4x6x2 for the fertilisation treatments, replicates and blocks, respectively). The layout of the experiment is shown in Appendix II. Each plot was harvested twice a year, in the last weeks in June and August. In 1996 the harvested biomass was oven-dried and weighed; in the analysis of the data the sum of both harvests was used. In June 1996 vegetation relevés were made before the harvest in all plots. A net surface area of 12 × 3 m was used to avoid edge effects. The relevés consisted of estimates of ground cover percentages per species (phanerogams, mosses and lichens). In our analyses we only used data collected in 1996 as the analysis of time series was not our aim in this study.

For each plot an index for nature conservation value was computed according to Van Dobben et al. (2015). This index is computed on the basis of the actual rarity and decline of each species over the past 20 years, and can be seen as an estimate of the probability of finding species that are listed in the International Union for the Conservation of Nature (IUCN) Red List (IUCN 2012) in this or a comparable vegetation type. We used data derived from the Dutch Red List (Van der Meijden et al. 2000) because it gives detailed quantitative information on rarity and decline (also compare Olde Venterink et al. 2003).

We used multivariate statistics to relate species composition per plot to the treatments. Total number of species, the nature conservation index and production were related to the treatments using ANOVA and multiple regression. We used the program package CANOCO (Ter Braak and Šmilauer 2012) for the multivariate analyses, and the program package GENSTAT V.17.1 (VSN International 2014) for all other statistical computations. Nomenclature of vascular plants follows Flora Europaea (Tutin et al. 1964–1980), in some cases supplemented by the Polish checklist of Mirek et al. (2002). For lichens we followed Smith et al. (2009), and for bryophytes Siebel and During (2006). For identification in the field we used the flora by Szafer et al. (1986).

Results

Production

Fertilization has a highly significant effect on production (ANOVA, P < 0.001) (Fig. 2). Production increases in the order O < K < PK ≤ NPK; the production levels in the PK and NPK treatments do not differ significantly (P > 0.05). The response ratios (production of fertilised plot / production of unfertilised plot; Harpole et al. 2011) are 3.35 for the K vs. O contrast and 2.78 for the PK vs. K contrast and thus fall well within Harpole et al.'s criterion for ‘biological significance’. The value of 1.08 for NPK vs. PK falls outside this criterion. The production in the O treatment is low even compared to unfertilised natural grasslands, and the production in the PK and NPK treatments is low compared to agricultural standards (compare e.g. Tallowin and Jefferson 1998).

Fig. 2
figure 2

Effect of the treatments on mean production (upper), number of species per plot (middle), and nature conservation value (lower). Error bars are 95 % confidence intervals; different letters indicate significant (P < 0.05) differences based on a t-test. Nature conservation value is on an arbitrary scale with interpretation: >15: high probability to find red-list species; 12–15: red-list species may be present; <12: very low probability to find red-list species

Species number and nature conservation value

Both the species number and the nature conservation value (determined by the presence of rare or declining species) are significantly influenced by the treatments. For both these variables the difference between the O and K treatments is not significant, and neither is the contrast between PK and NPK (Fig. 2). The value of both variables increases in the order NPK ≤ PK < K ≤ O. Both species number and conservation value decrease highly significantly (P < 0.001) with increasing production (Fig. 3, Table 2). In a regression model, species number and conservation value are almost equally well explained by either production, the fertilisation treatments, or their combination; the only exception being that the conservation value is slightly (but significantly, P ≈ 0.01) better explained by the treatments (Table 2).

Fig. 3
figure 3

Relation between production and number of species (upper) or nature conservation value (lower) per plot. Symbols indicate treatments. Regression equations are: N species =34.6*** + 1.05ns X (Block dummy) - 1.32*** X Production (62 % explained variance) Cons value =13.5*** + 0.43ns X (Block dummy) - 0.55*** X Production (71 % explained variance). Block effects are included as a dummy variable with value 1 in Block 1 and 0 in block 2. Superscripts indicate significance of regression coefficients determined by a t-test, ns: P > 0.05, ***: P < 0.001. See Fig. 2 for the interpretation of conservation value

Table 2 Number of species and nature conservation value explained by production, fertilisation and their combination. The figures in each row give the effect of the change indicated in the first column (with +, addition of a term to the model in the row above, for the first row this is a null model of only block effects; and -, deletion of a term from the model in the row above). F = (difference in regression mean square due to given change) / (error mean square), P = probability of a higher F-value under the null hypothesis determined by ANOVA

Multivariate approach

A preliminary analysis showed that the heterogeneity in the vegetation data was rather limited (‘gradient length’ sensu Ter Braak and Šmilauer (2012) is 2.5 for the first axis in DCA) and we therefore used linear techniques (PCA, RDA) for our analyses. The general pattern that emerges from a PCA ordination of all plots is that the first axis (the most important one in terms of explained variance) represents the contrast O and K vs. PK and NPK, the second axis represents the O vs. K contrast, and the third axis the PK vs. NPK contrast (Figs. 4 and 5). No ecological interpretation could be found for the fourth axis.

Fig. 4
figure 4

Ordination results: sample plot of first and second PCA axis. The blocks were used as a covariable. Data points correspond to individual plots, colours to treatments, and symbol size to production (upper), number of species (middle) and nature conservation value (lower; see Fig. 2 for its interpretation)

Fig. 5
figure 5

Ordination results: species plot. The upper plot gives the first and second PCA axis, the lower the first and third axis. Eigenvalues are 0.42, 0.16, and 0.04 for the first, second and third axis, respectively, i.e. the two plots together explain 63 % of the variance in the species data (after accounting for block effects). The treatments are indicated as their centroids (triangles); a notional arrow can be drawn from the origin to each species name, the projection of a centroid on an arrow is the fitted abundance of that species in that treatment, with scaling: origin = mean, head of arrow = mean + standard deviation (and mirror-image of arrow head with respect to origin = mean - standard deviation). Note that the centroids for the treatments are those of the respective sample scores in PCA and not in RDA (which was only used to determine the strengths of the treatment effects and not to make a picture of the effects themselves). Only species with at least three occurrences are displayed, some species names in the left part of the lower plot have been removed to increase readability. Explanation of abbreviated species names (phanerogams first letter uppercase, cryptogams first letter lowercase, * indicates that a species has been sown at the start of the experiment) Achilmil, Achillea millefolium; Agrossto, Agrostis stolonifera*; Anthrsyl, Anthriscus sylvestris; Arrheela, Arrhenatherum elatius; Betulpub, Betula pubescens; bractalb, Brachythecium albicans; Bromuine, Bromus inermis*; Campapat, Campanula patula; Capsebur, Capsella bursa-pastoris; Cardnare, Cardaminopsis arenosa; Carducri, Carduus crispus; Carexnig, Carex nigra; Carexpan, Carex panicea; Carexser, Carex serotina; ceradpur, Ceratodon purpuraeus; Cerasfon, Cerastium fontanum; Cerassem, Cerastium semidecandrum; Cirsiarv, Cirsium arvense; cladochl, Cladonia chlorophaea; cladocoi, Cladonia ochrochlora; cladofim, Cladonia fimbriata; cladopyx, Cladonia pyxidata; cladorei, Cladonia rei; climaden, Climacium dendroides; Convoarv, Convolvulus arvensis; Dactyglo, Dactylis glomerata*; Deschces, Deschampsia cespitosa; Elymurep, Elymus repens; Epiloros, Epilobium roseum; Erigecan, Conyza canadensis; Festuaru, Festuca arundinacea*; Festuo-O, Festuca ovina; Festupra, Festuca pratensis*; Festur-C, Festuca rubra*; Fragaves, Fragaria vesca; Galeotet, Galeopsis tetrahit; Galiumol, Galium mollugo; Heracsos, Heracleum sosnowskyi; Hierapil, Hieracium pilosella; Inulabri, Inula britanica; Lamiup;P, Lamium purpureum; Leontaut, Leontodon autumnalis; Leucavul, Leucanthemum vulgare; Linarvul, Linaria vulgaris; Lotuscor, Lotus corniculatus*; Luzulcam, Luzula campestris; Luzulmul, Luzula multiflora; Lysimvul, Lysimachia vulgaris; Lythrsal, Lythrum salicaria; Mediclup, Medicago lupulina; Phalaaru, Phalaris arundinacea; Phleupra, Phleum pratense*; plagmcus, Plagiomnium cuspidatum; Plantlan, Plantago lanceolata; Poa pra, Poa pratensis*; Poa pal, Poa palustris*; polymjun, Polytrichum juniperinum; Polynavi, Polygonum aviculare s.l.; Polynper, Polygonum persicaria; Potenans, Potentilla anserina; Potennor, Potentilla norvegica; Prunevul, Prunella vulgaris; Ranunacr, Ranunculus acris; Ranunaur, Ranunculus auricomus; Ranunfla, Ranunculus flammula; Ranunrep, Ranunculus repens; Rhamncat, Rhamnus catharticus; Rumexace, Rumex acetosa; Rumexact, Rumex acetosella; Saginnod, Sagina nodosa; Silenl-A, Melandrium album; Soncharv, Sonchus arvensis; Stellaqu, Myosoton aquaticum; Stellgra, Stellaria graminea; Stellmed, Stellaria media; tarax-sp., taraxacum officinale; torturur, Syntrichia ruralis; Tragopra, Tragopogon pratensis; Triforep, Trifolium repens*; Urticdio, Urtica dioica; Veronarv, Veronica arvensis; Veroncha, Veronica chamaedrys; Viciacra, Vicia cracca; Violaarv, Viola arvensis; Violacan, Viola canina.

The strength of the effect of the treatments on the species composition of the plots was determined by forward selection in RDA (the ‘canonical’ form of PCA, cf. Ter Braak and Šmilauer 2012). Again the effect of all contrasts is highly significant except PK vs. NPK which is nearly significant (P ≈ 0.06) (Table 3).

Table 3 Significance of the effect of treatment on species composition: result of forward selection of treatment variables in RDA. F = (difference in regression mean square due to given contrast) / (error mean square); P = probability of a higher F-value under the null hypothesis as determined on the basis of 9999 bootstrap samples; % explained variance = increase in explained variance on inclusion of the given contrast in the model

Figure 4 is the sample plot for the first two axes resulting from PCA, and shows the production, number of species, and conservation value of each plot. The data points belonging to the O, K, and [N]PK treatments form three very distinct clusters, thus illustrating the strong effect of the treatments on the species composition. These plots also clearly show the increase of production in the direction K < PK < [N]PK, and the concomitant decrease in both species number and nature conservation value. For conservation value the decrease is more outspoken than for species number, but for both indices the difference between the O and K treatments is only small.

Response of individual species to the treatments

In the PCA biplots there appear to be three rather distinct groups of species that are associated with the O, K and [N]PK treatments (Fig. 5). The O treatment is mainly characterised by the absence of grasses (except Festuca spp.) and the presence of Carex nigra, Rumex acetosa and a high number of cryptogamic species (e.g. Cladonia spp., Polytrichum juniperinum or Ceratodon purpureus). The plots of the other treatments mainly contain meadow species and weedy species. Those of the K treatment are often species of oligotrophic or mesotrophic grassland (e.g. Phleum pratense, Luzula campestris, L. multiflora, Plantago lanceolata, Hieracium pilosella, Carex panicea or Stellaria graminea); while in the PK and NPK treatments species of more eutrophic grassland occur (e.g. Arrhenatherum elatius, Poa pratensis, Bromus inermis or Urtica dioica). Grasses like Lolium perenne, which is characteristic of highly productive eutrophic grassland, are however absent from our plots. The most obvious floristic difference between the PK and NPK treated plots is that the NPK plots contain more weedy species (e.g. Stellaria media, Conyza canadensis, Viola arvensis, Lamium purpureum) and the PK plots contain more grasses (e.g. Poa pratensis, Arrhenatherum elatius). However, the floristic difference between the latter two treatments is small (Table 3).

Discussion

Production

The most conspicuous features of our experiment are the virtual absence of any effects of nitrogen addition, and the very high nitrogen availability due to mineralisation of the drained peat (Okruszko et al. 1988). If nitrogen and phosphorus are considered, our experiment seems to belong to the less common case of ‘single-nutrient limitation’, although sub-additive independent co-limitation in the sense of Harpole et al. (2011) cannot be ruled out (in that case the production of hypothetical K + N plots would be equal to the production of K + P and K + P + N plots). However, if potassium and phosphorus are considered there is co-limitation or at least serial limitation, as potassium alone significantly increases production; the combination of phosphorus and potassium still further (and again significantly) increases production.

Species number and nature conservation value

Both species number and nature conservation value appear to decrease almost linearly with increasing production, but conservation value is slightly better explained by the treatments than by production (Table 2). This is an indication that the presence of endangered species is at least partly determined by the individual nutrients and not exclusively by production. The negative relationship between species number and production has been found in many other studies (cf. Huston 2014), however the ‘hump-back’ relationship between species number and production postulated by Grime (1973, 1979) is only weakly present in our data, even though our relationship can be seen to begin at extremely low production values. In our plots the highest species number (37 species in a plot of 36 m2) was found at a low (1.74 t.ha−1), but not the lowest (0.31 t.ha−1) production. Our results strongly disagree with Laliberté et al. (2014) who, in a study on a soil chronosequence, found soil factors (and most notably pH) to be the key factor that determines species richness, independent of production.

Species composition

The effect of the treatments on the species composition of the plots appears to be very large. The eigenvalues (which in our scaling can be seen as fractions explained variance) are 0.42, 0.16 and 0.04 for the first, second and third axis respectively, which is high for ecological data that usually have percentages explained variance in the order of 10 % for the first axis. However, such values are not uncommon for designed experiments (Van Dobben et al. 1999). The total percentage explained variance (59 %) is also high for ecological data but not uncommon in this type of experiment.

The slightly better prediction of nature conservation value by the treatments than by production suggests that the variation in species composition is not purely one-dimensional as would be expected if it were only determined by production. This is confirmed by the PCA results (Fig. 4) that showed a clearly two-dimensional variation with three distinct clusters of plots associated with the O, K and [N]PK treatments respectively. The first axis, which represents the most important direction of variation (its eigenvalue is far higher than the eigenvalue of the second axis) coincides with the +P vs. -P contrast, while the second axis coincides with the O vs. K contrast. A third, less distinct dimension coincides with the +N vs. -N contrast (Table 3, Fig. 5). The K treatment leads to a significantly higher production compared to the O treatment, but not to a significantly lower species number or a lower conservation value (Figs. 2 and 3). This means that in the K treated plots the number of [endangered] species is higher than expected on the basis of production. Or, in other words, the number of [endangered] species seems to depend on the +P vs. -P contrast rather than on production. The species composition of K treated plots is significantly different from the O plots (Table 3), but the difference entails a shift in species composition (mainly from cryptogams to grasses; Fig. 5) rather than a loss of species.

Fate of the added nutrients in the soil

In a balance study of the PK fertilised plots over the period 1958–1982, Okruszko et al. (1988) estimate the yearly N mineralisation to be 674 kg N ha−1 y−1, based on measurements of peat subsidence and N content of the peat (Table 4). Further, they estimate the uptake by the crop (and subsequent removal by hay) to be 205 kg N ha−1 y−1. The latter figure may seem high (at a yearly production of 6.3 t DM ha−1 (Fig. 2) this would mean a N content of 3.1 %). However, the productivities of the plots have substantially decreased since the start of the experiment in 1958, when the O plots produced c. 7 t DM ha−1 and the fertilised plots c. 12 t DM ha−1 (Gotkiewicz and Gotkiewicz 1987). Also the N mineralisation reported by Okruszko et al. (1988) seems very high compared to values in the order of 100 kg N ha−1 y−1 for drained peat in the same area reported by Olde Venterink et al. (2009). Okruszko et al. (1988) estimate that of the 469 kg N ha−1 y−1 that is not taken up by the crop, 63 kg N ha−1 y−1 is lost by leaching and denitrification, while the remainder is again incorporated in organic matter in ‘muck’ (degraded, mineralised peat) (Table 4). As a result, the total N content of muck is c. 30–35 % higher than in the original peat that formed it (Okruszko 1960).

Table 4 Balance of N, P and K in the PK fertilised plots over the period 1958–1982 (after Okruszko et al. 1988). Unit: kg element ha−1.y−1. The numbers are the source and sink terms of N, P and K in the soil, a positive balance difference indicates re-incorporation into organic matter and a negative value indicates depletion

In contrast to N, there is a strong depletion of K in the course of the experiment: in the PK-fertilised plots the total K content of the soil in 1982 had been reduced to 12 % of its value at the start in 1958. Therefore the unfertilised plots can be assumed to be strongly K-limited, and even the K-fertilised plots may still be K-limited. In general, the removal of K with hay in annually mown grassland on sandy or peaty soil is high, which makes such ecosystems sensitive to K-limitation (Table 4; Van der Woude et al. 1994; Van Duren and Pegtel 2000). For P this is less obvious. In the fertilized plots P accumulates in the soil, although possibly in a form unavailable to plants, through binding by amorphous aluminium and iron oxides (Okruszko et al. 1988; Freese et al. 1992). Even in the unfertilised plots the mineralisation of P still exceeds the removal through hay (Table 4). Note that the relative effect of K fertilisation on production exceeds that of P fertilisation (response ratios sensu Harpole et al. (2011) are 3.35 and 2.78, respectively; Fig. 2).

Comparison to other N addition experiments

In a stable situation in grassland on sandy soils the N mineralisation should balance the organic matter left after hay cutting (Wamelink et al. 2009). If for example 1 t ha−1 is left, containing 2.5 % N, mineralisation would amount to a maximum of 25 kg N ha−1 y−1. If this amount is increased by some 10–20 kg N ha−1 y−1, e.g. by deposition, a conspicuous loss of species usually occurs (Bobbink and Hettelingh 2011). Stevens et al. (2004, 2010) and Maskell et al. (2010) report a decrease from c. 20 species per plot (4 m2) to less than 10 if deposition increases from c. 10 to c. 40 kg N ha−1 y−1 for acid grassland in north-western Europe. In our plots species numbers are c. 33 species per plot (36 m2) in plots that do not receive P and c. 25 in the P-fertilised plots (Fig. 3). These numbers are high compared to those given by the authors cited above, which may be partly due to the larger plot size in our study. However, our mineralisation figures exceed the deposition values used in the above-cited studies by almost an order of magnitude. Extrapolation of Stevens’ (2010: Fig. 1) regression line to our mineralisation value of c. 500 kg N ha−1 y−1 would result in far less than one species per plot (of 4 m2), and the value of 100 N ha−1 y−1 reported by Olde Venterink et al. (2009) would result in only two species, while in fact the minimum encountered in our study is 20 species per plot (of 36 m2). Our results therefore suggest that the decrease in species number with increasing N deposition is not a direct effect of N availability but rather an indirect effect that comes about through, for example, production.

Effect of N vs. P on species number and nature conservation value

Olde Venterink et al. (2003) and Wassen et al. (2005) found indications that at a given production the number of threatened species (i.e. appearing on the Dutch Red List: Van der Meijden et al. 2000) is higher under P limitation than under N limitation. We cannot use our data as a test of this hypothesis as we did not find any N limitation; however we did find a comparable phenomenon. The K-fertilised plots have a significantly higher production than the unfertilised ones, but there is no significant difference with the O treatment in relation to species richness or in nature conservation value (Figs. 2, 3 and 4). Both the unfertilised and the K-fertilised plots can be assumed to be P (co-) limited, while in the PK and NPK-fertilised plots this P-limitation is (at least partly) lifted. These plots have both a significantly higher production and a significantly lower species richness and conservation value (Fig. 2). Apparently lifting the P-limitation causes a decrease in the number of species while lifting the K-limitation does not. Note that N-limitation is absent in our data and present only at high production levels in Olde Venterink et al.'s (2003) data. An explanation may be that there are more forms of P in the soil and, consequently, more strategies to acquire P (such as mycorrhiza or exudation of phosphatase) than to acquire N or K, both of which are highly mobile in soils (Olde Venterink 2011). There is a striking resemblance between our species richness vs. production plot (Fig. 3) and Olde Venterink et al.'s (2003) plots for P-limited situations (i.e. high maximum species number at very low production and sharp decline at higher production levels, and absence of production levels above c. 6 t ha−1). This suggests that the same mechanisms are indeed active in spite of very different N availabilities.

General conclusions

In relation to whether species richness is primarily determined by production, by abiotic conditions or by a combination (see e.g. Gough et al. 1994; Grace 1999; Olde Venterink et al. 2001; Laliberté et al. 2014), our results are somewhat contradictory. We found slightly different effects for K and P additions on numbers of [endangered] species and on production; which is an indication that species richness not only depends on production, but at least partly on the relative availabilities of these nutrients. However, in our experiment addition of N does not have any effect on species richness, which is in sharp contrast to the majority of other N addition experiments where N leads to a sharp decline in species numbers. In our view this can only be explained by assuming that our case is an exceptional one where nutrients other than N are limiting production. This hypothesis is supported by the absence of a significant effect of N addition on production. However, for this hypothesis to be true, one has to assume that species composition depends entirely on production. Taking everything together the most probable hypothesis is that the number of [endangered] species is primarily determined by production, but can be modified by the nature of the limiting nutrient(s). Olde Venterink et al.’s (2003) hypothesis of a higher species number under P limitation compared to N limitation at a given production is coherent with the above hypothesis; however, this would need to be confirmed with further testing.

Implications for nature conservation

There is overwhelming evidence, both experimental and observational, that nitrogen enrichment leads to a loss of species, at least for vegetation in the terrestrial environment; this evidence is summarised by Bobbink and Hettelingh (2011) or De Schrijver et al. (2011). Species loss through N enrichment appears to be widespread although it is most prominent in the temperate zone (Sala et al. 2000). This agrees with the widespread occurrence of N limitation that has been shown by the large-scale meta-studies cited earlier. If P is important as a limiting nutrient worldwide, as put forward by Wassen et al. (2005), it will usually be co-limiting so that even in that case addition of N will lead to increased production and loss of species. It can therefore be concluded that there is little scientific evidence for Wassen et al.'s (2005) view that the occurrence of a widespread limitation by P alone, together with a strong P enrichment of the environment, makes N emission reduction an ineffective tool for preserving biodiversity.

We also observed the sharp decline in species richness as a response to the lifting of P-limitation, even at a moderate increase of production, that was previously seen by Olde Venterink et al. (2003). In practice this may mean that if P is added to a species-rich, P-(co)limited site, a strong loss of botanical diversity may occur, irrespective of effects of N. As P is highly immobile in soils this is most likely to occur in sites that are influenced by external supply of surface-water such as fens or mires. However, spreading of fertiliser dust (Ceulemans et al. 2014) or chemical changes in the soil that increase P availability (e.g. acidification; Hinsinger 2001) may also play a role. Our results also lend support to the hypothesis that unfertilised hay meadows on peaty soil may become strongly K-limited in the long run (e.g. Van Duren and Pegtel 2000; see also Fay et al. 2015), but we did not find indications that this will lead to a higher species richness. In cases with a very strong K and/or P-limitation (like ours) N emission reduction would indeed be ineffective for preserving diversity (as has been put forward by Wassen et al. 2005). However, we believe that the case reported here is extreme due to the very strong mineralisation that is in turn caused by drainage of the peat (groundwater in summer is 85–95 cm below surface level). In practice such deep drainage of peat soils is never applied in natural areas because the resulting mineralisation and concomitant soil subsidence are considered highly unwanted phenomena (Borger 1992). We conclude that N emission reduction is likely to be an effective strategy for preserving or restoring botanical diversity, even though limitation by P and/or K may be common. However, an increase in P availability due, for instance to application of fertiliser or influx of nutrient-rich surface water will most probably lead to a further loss of species. In general our data lend experimental support to the view of Ceulemans et al. (2014) that the general decrease of species numbers in the grasslands north-western Europe is not only due to enrichment by N but also by P.