Direct and plant community mediated effects of management intensity on annual nutrient leaching risk in temperate grasslands

Grassland management intensity influences nutrient cycling both directly, by changing nutrient inputs and outputs from the ecosystem, and indirectly, by altering the nutrient content, and the diversity and functional composition of plant and microbial communities. However, the relative importance of these direct and indirect processes for the leaching of multiple nutrients is poorly studied. We measured the annual leaching of nitrate, ammonium, phosphate and sulphate at a depth of 10 cm in 150 temperate managed grasslands using a resin method. Using Structural Equation Modeling, we distinguished between various direct and indirect effects of management intensity (i.e. grazing and fertilization) on nutrient leaching. We found that management intensity was positively associated with nitrate, ammonium and phosphate leaching risk both directly (i.e. via increased nutrient inputs) and indirectly, by changing the stoichiometry of soils, plants and microbes. In contrast, sulphate leaching risk was negatively associated with management intensity, presumably due to increased outputs with mowing and grazing. In addition, management intensification shifted plant communities towards an exploitative functional composition (characterized by high tissue turnover rates) and, thus, further promoted the leaching risk of inorganic nitrogen. Plant species richness was associated with lower inorganic nitrogen leaching risk, but most of its effects were mediated by stoichiometry and plant community functional traits. Maintaining and restoring diverse plant communities may therefore mitigate the increased leaching risk that management intensity imposes upon grasslands.


Introduction
Grasslands are widespread and highly diverse ecosystems providing a multitude of ecosystem functions and services, including nutrient cycling and food production (Bengtsson et al. 2019). Management intensification of grassland ecosystems places them at considerable risk of nutrient losses via leaching, especially for nitrate  and references therein), which represent the loss of important resources and, at the same time, potential threats to human health and aquatic ecosystems (EEA 2018; EU-Nitrate Directive 2018; WHO 2011). Compared to nitrogen (N), management effects on leaching of other elements, like phosphorus (P) and sulphur (S), are less studied in grasslands, despite their agricultural and environmental importance (Eriksen 2009;Withers et al. 2014;Gallejones et al. 2012). Unravelling the relationships between ecosystem management, biotic drivers and nutrient leaching will provide important insights for sustainable grassland management and could help secure associated ecosystem services.
In unfertilized grasslands, inorganic nutrients that are available for leaching represent the net balance between nutrient inputs to the soil solution (e.g. by organic matter mineralization or dissolution) and outputs from the soil solution (e.g. due to uptake by plants, or immobilization in microbial biomass), as well as nutrient exchange with the soil matrix (i.e. via (de)sorption to minerals). Management, such as fertilization and grazing, changes these relations by increasing both inputs and outputs of nutrients from the ecosystem. In the past, grasslands were considered to have low N leaching fluxes, and introducing temporary (mowed) grasslands in arable crop rotations can even reduce nitrate leaching (Kunrath et al. 2015). However, studies investigating the effects of grassland management on soil nitrate show increased leaching risk in pastures compared to meadows (Di and Cameron 2002;Ryden et al. 1984), and when broader management gradients are included, both fertilization and grazing equally increase nitrate leaching risk . Grassland intensification can affect nutrient cycling directly by changing inputs and outputs (i.e. fertilization, grazing and harvest, Rumpel et al. 2015) and indirectly by changing plant communities (e.g. loss of diversity or functional change) and altering interactions between soils, plants and soil microorganisms (de Vries et al. 2012;Klaus et al. 2018). However, at present, it is unclear if management intensity predominantly drives leaching risk directly and/or indirectly by changing plant and/or microbial communities.
Grassland management, mainly fertilization, mowing and grazing, influences soil, plant and microbial properties (Herold et al. 2014;Kleinebecker et al. 2014;Boeddinghaus et al. 2019) with potential effects on nutrient leaching. For instance, management intensification reduces the C:N ratio of grassland soils and increases enzymatic activities (Herold et al. 2014) and, potentially, the decomposition of organic matter. In addition, fertilization increases plant biomass productivity and its nutrient concentrations, reducing its C:N stoichiometry, and it reduces plant diversity and root litter decomposition (Gilhaus et al. 2017; Vol.: (0123456789) Solly et al. 2014). Many studies have reported that a high plant species diversity reduces N and P availability and leaching in grasslands (de Vries et al. 2016;Leimer et al. 2015Leimer et al. , 2016Oelmann et al. 2011;Tilman et al. 1996), especially when compared to monocultures (Scherer-Lorenzen et al. 2003), possibly due to resource complementarity (Kahmen et al. 2006). However, plant diversity effects on nitrate leaching are still debated and have not yet been comprehensively assessed up to now, as most of the respective studies were conducted in experimental grasslands or mesocosm experiments. Such studies are useful to gain a mechanistic understanding of the underlying processes controlling soil and plant interactions. However, they often do not consider realistic management and/or environmental gradients. Thus, they might fail to capture the complexity of real-world ecosystems (Klaus et al. 2020a) since, for instance, management is an important determinant of plant community diversity and composition Midolo et al. 2019). To overcome these constraints, and gain knowledge transferable to sustainable grassland management, we need to investigate grassland ecosystems along broad gradients of management intensity and environmental conditions .
Plant community composition, e.g. the abundance of functional groups or specific traits, can also play an important role in grassland nutrient cycling (Boeddinghaus et al. 2019;de Vries and Bardgett 2016). Plant functional groups are often used as predictors of nitrogen leaching. The group of legumes is, for example, associated with higher nitrate leaching, while the group of grasses indicates lower nitrate leaching (de Vries and Bardgett., 2016;Leimer et al. 2015Leimer et al. , 2016. In addition, plant traits can be used to categorize plant communities as fast-growing, and exploitative or slow-growing and conservative (Bruelheide et al. 2018;Reich 2014;Díaz et al. 2016). Exploitative plant communities have higher specific leaf area and leaf N (and P) and lower leaf dry matter content than conservative communities (Díaz et al. 2016;Wright et al. 2004). Regarding belowground traits, exploitative communities have higher specific root length and lower root tissue density and fine root diameter than the conservative communities (Bergmann et al. 2020). Plant, and especially root, traits are thus broadly indicative of nutrient acquisition strategies of plant species and these are relevant for nutrient cycling in soils (de Vries and Bardgett, 2016). For instance, a high specific leaf area indicates high photosynthetic capacity, which relates to fast growth rates, short tissue lifespan and high decomposability, and high needs for nutrients (de Vries and Bardgett, 2016). Moreover, a high specific root length indicates a better screening of soil for nutrients, and a low fine root diameter indicates high nutrient uptake efficiency (Bardgett et al. 2014;Freschet et al. 2021). In managed grasslands, plant community traits are influenced by management-induced shifts in plant community composition (Allan et al. 2015;Busch et al. 2018) and might explain the effects of vegetation change on nutrient leaching.
Recent evidence highlights the importance of soil microbial communities and their properties for N cycling (Grigulis et al. 2013;Pommier et al. 2018). Understanding, and even controlling, plant-microbe interactions could provide opportunities for enhanced N retention in grasslands, with desired effects both on biomass production and N leaching mitigation. Studies on element cycling that include soil microbial properties are, in general, scarce compared to those addressing abiotic and/or plant effects (Pommier et al. 2018), and so are studies that simultaneously study the cycling of multiple nutrients in managed grasslands. However, management strongly changes microbial properties since, for instance, fertilization shifts the microbial community structure towards bacteriadominated communities (de Vries et al. 2006;Donnison et al. 2000;Herold et al. 2014), which are characterized by lower microbial C:N ratios and higher mineralization rates (Six et al. 2006;Strickland and Rousk, 2010). To elucidate the drivers of ecosystem (dis)services like nutrient leaching, a holistic approach that includes management, but also soil, plant and microbial effects and their interactions is needed (Koncz et al. 2015;Pommier et al. 2018).
Reliable assessments of the drivers of nutrient leaching need to cover long time spans, for example due to the large variability in weather conditions and the timing of land management activities (de Vries et al. 2011;Klaus et al. 2018). This is even more important since recent evidence suggests that plant diversity effects on N leaching follow seasonal patterns with more pronounced effects in fall or winter (Compton et al. 2021;Leimer et al. 2015), seasons that are often ignored as most studies focus on the growth period only. However, accounting for both high spatial and temporal variability imposes considerable logistical constraints. Resin methods for the evaluation of nutrient availability and leaching have received increasing attention in the literature (Grahmann et al. 2018;Grunwald et al. 2020;Predotova et al. 2011;Willich and Buerkert, 2016), and are often successfully used to integrate the temporal effects of different drivers on nutrient leaching (Klaus et al. , 2020b. In this study, we measured the annual downward flux of nitrate, ammonium, phosphate and sulphate with the percolating water, i.e. leaching, at a depth of 10 cm in 150 grassland sites from spring 2018 to spring 2019. The grassland sites were located in three regions in Germany as part of the Biodiversity Exploratories project (Fischer et al. 2010), and covered broad management intensity gradients within each region. We used a resin method, namely the Self-Integrating Accumulators (SIAs), to determine the cumulative annual leaching, and not just the availability of nutrients in the soil solution (Bischoff, 2009). Considering that grassland management, here fertilization and grazing, has direct (via increased inputs and outputs) and indirect (via soils, plants and microbial biomass) effects on nitrate, ammonium, phosphate and sulphate leaching, we addressed the following hypotheses: 1. Fertilization and grazing intensity increase nutrient leaching risk from topsoils due to increased nutrient inputs (direct effects).
2. Fertilization and grazing intensity increase nutrient leaching risk by lowering the C:N stoichiometry of soils, plants and soil microbial biomass (indirect effects). 3. Fertilization and grazing intensity increase nutrient leaching risk (i) by shifting the functional composition of plant communities from conservative (i.e. slow-growing) communities at extensively managed grasslands to exploitative (i.e. fast-growing) communities at intensively managed grasslands, and (ii) by reducing plant species richness (indirect effects).

Study regions
This study was conducted in the Biodiversity Exploratories project, which investigates land use effects on biodiversity and ecosystem functions and services in real-world ecosystems (Fisher et al. 2010). The project includes three study regions in Germany; the Schwäbische-Alb (ALB), the Hainich-Dün (HAI) and the Schorfheide-Chorin (SCH). The study regions differed in their geology, topology, dominant soil types and/or climate (Table 1). Soils in ALB developed mainly on Jurassic limestone and were clay-rich Leptosols or Cambisols (IUSS Working Group WRB, 2014). In HAI, soils had a loamy or clayey texture due to the dominant geological substrate of loess over limestone and the main soil types were Cambisols, 6.0-7.0 6.5-8.0 8.0-8.5 Annual precipitation (mm) 700-1,000 500-800 500-600 Stagnosols and Vertisols. Soils in SCH were drained Histosols with a loamy texture, Gleysols, Cambisols, Luvisols and Albeluvisols. In each study region, 50 plots of 50 m × 50 m were established within agriculturally managed permanent grasslands (Fisher et al. 2010). Information regarding the regional management practices as well as soil, plant and microbial properties is given in Table 2 (and Tables S1 and S2). Management intensity ranged from extensively managed grasslands to intensively grazed and/or fertilized pastures and meadows, covering a wide and realistic management intensity gradient. Fertilization, expressed as N additions, was applied as organic or mineral (chemical) fertilizers. Grazing, expressed as equivalent-livestock unit per area multiplied

Soil properties
In May 2017, we collected 14 soil samples of the upper 10 cm along two intersecting 20 m transects in each of the 150 grasslands. Soil samples were used to prepare one composite sample per grassland. All soil samples were air-dried and sieved to < 2 mm, while a portion was ground for elemental analysis. Total carbon (C) and total N were determined by dry combustion at 1,100 °C with an elemental analyser VarioMax (Elemental, Hanau, Germany). Soil inorganic C was determined through the same process after removing organic carbon (OC) by exposing 250 mg of soil to 450 °C for 16 h. Soil OC was calculated as the difference between total and inorganic C. Soil pH was measured in the extraction solution of 10 g of soil with 25 mL of CaCl 2 (0.01 M) with a glass-electrode and a pH meter (WTW pH meter 538, Weilheim, Germany). Total soil S concentration was determined with ICP-OES (Acros, SPECTRO Analytical Instruments GmbH, Kleve, Germany) in the extracted solution of an aliquot of 1.5 g of soil with 3.5 mL HNO 3 and 10.5 mL HCl. Soil bulk density was determined by drying a known volume of soil collected in May 2014 at 105 °C to constant weight. Soil texture analysis of composite samples collected in May 2011, was performed by destruction of the organic matter with hydrogen peroxide, dispersion of soil and clay aggregates into discrete units, and separation of particles of different sizes by sieving and sedimentation (DIN ISO 11277 2002;Solly et al. 2014). Soil oxalate extractable Al (in the composite samples collected in May 2011) was determined with ICP-OES (Optima 3300 DV, PerkinElmer) as a measure for non and poorly crystalline oxi-hydroxides. In soil samples collected in May 2014, we extracted NaHCO 3 -P and NaOH-P with subsequent extractions and we determined inorganic P with a continuous flow analyser (CFA, AA3, XY2, Seal-Analytic, Norderstedt, Germany). Labile and moderately labile inorganic P (P labile,inorg ) was calculated as their sum (Sorkau et al. 2018).

Microbial soil properties
Aliquots of sieved soil (from the composite samples collected in May 2017) were used for microbial analysis. We determined the microbial C and N concentration with the chloroform-fumigation-extraction method (CFE), according to Vance et al. (1987) and Keil et al. (2011). Fumigated sample replicates were incubated with CHCl 3 for 24 h. C and N were extracted from each fumigated and non-fumigated replicate (5 g) with 40 ml 0.5 M K 2 SO 4 . The suspension was horizontally shaken (30 min, 150 rpm) and centrifuged (30 min, 4400 × g). C and N concentrations in diluted extracts (1:4, extract:deionized H 2 O) were measured with a TOC/TN analyser (Multi N/C 2100S, Analytik Jena AG, Jena, Germany). No correction factors accounting for the extractable fraction of microbial C and N were used to calculate microbial biomass.

Plant properties
Plant communities were surveyed in May 2018. Standing biomass was measured in eight 0.5 m × 0.5 m subplots in each grassland, which were fenced to prevent livestock grazing or mowing prior to harvest. Aboveground vegetation was clipped at 2 cm height and dried at 80 °C for 48 h. Dry mass was then weighed at one decimal of the gram, and ground for chemical analysis. In a nearby 4 m × 4 m subplot, we recorded all vascular plant species and estimated their percentage cover. Afterwards, we calculated the plant species richness, and the cover of two main functional groups, legumes, which typically have nitrogen rich tissues, and grasses. We determined the concentrations of C, N and P in aboveground plant biomass with a Near-Infrared Reflectance spectrometer (Spec-traStar 2400, Unity Scientific, Columbia, MD, USA) Kleinebecker et al. 2011). The reflectance spectrum of each sample, averaged over 24 scans, was recorded from 1250 to 2350 nm with 1 nm intervals.
We used plant species identity, and their abundances in our grasslands, to calculate the Community Weighted Mean (CWM) traits, for four above-(specific leaf area, SLA; the leaf dry matter content, LDMC; leaf N and P) and three belowground (fine root diameter, fineRD; root tissue density, RTD; specific root length, SRL) traits (Table 2). Aboveground traits were extracted from the TRY database (Kattge et al. 2020), and belowground traits were measured in a greenhouse experiment on plant species from the Biodiversity Exploratories grasslands (Lachaise et al. 2021). First, we matched the species identities found in our grassland sites to the TRY database (including synonyms) and excluded (i) duplicates, (ii) non-mature plants, (iii) non-healthy or dead plants, (iv) measurement contributions from experimental manipulations, and (v) measurement contributions where mean data were not reported. Then, for each species-trait combination, we averaged the trait value for each contributing dataset and excluded outlying datasets. We calculated average trait values for each species over all contributing datasets, weighted by the number of replicates in each dataset. Finally, the community weighted mean was calculated as the average trait values of all species, weighted by their relative abundance.
In our grasslands, the cover of legume species was negatively correlated with the CWM leaf dry matter content and specific root length, and, to some extent, it was positively correlated with leaf N and fine root diameter ( Fig. 1). Moreover, the cover of grass species was negatively correlated with fine root diameter and root tissue density, and positively correlated with specific root length. Specific leaf area was not related to the plant functional group of legumes or grasses. Grazing intensity was negatively correlated with specific root length and leaf dry matter content, and was positively related to leaf N (and P). Fertilization was positively correlated with specific leaf area, and, to some extent, negatively related to leaf dry matter content and root tissue density.

Annual nutrient leaching
Annual, downward fluxes of nitrate, ammonium, phosphate and sulphate with the percolating water, i.e. leaching, were measured with the Self-Integrating Accumulators (SIAs) provided by the TerrAquat GmbH (Nürtingen, Germany, www. terra quat-gmbh. com; Bischoff, 2009). SIAs consisted of a PVC cylindrical cartridge (diameter and height: 10 cm) with a net at the lower end and were filled with a mixture of cation and anion exchange resins with a silica sand of specified grain size distribution.
SIAs were installed below the upper 10 cm of soil, and we measured the leaching from the upper 10 cm of soil for three reasons. Firstly, we wanted our leaching data to match the soil and microbial properties data, which were also measured at the upper 10 cm of the soil (Table 2). Secondly, we wanted to follow a standardized methodology in every grassland, because the total soil depth varied among and within the three study regions (Table 2), and it was often ≤ 20 cm deep (in 58 out of the 150 grasslands). Finally, we were interested in biotic drivers of nutrient leaching, which are probably easier identified below the biologically most active topsoil layer than in the deeper subsoil, where geochemical drivers might become more important. Since we measured nutrient leaching from the upper 10 cm of soil, and not below the rooting zone, we cannot conclude on nutrient losses from the ecosystem, but our data indicate their loss risk.
We installed four SIAs in each of the 150 grasslands (total of 600 SIAs) in spring 2018, a year in which a summer drought occurred in central Europe (Apostolakis 2022). We selected four 1 m × 1 m subplots at a 2 m extension of the cardinal transects used for soil sampling. In each subplot, we opened a 40 cm × 40 cm pit down to 30 cm. We selected the wall that was against the slope and, on this wall, we  Table 2), two plant functional groups (i.e. cumulative cover of Legumes and Grasses) and management intensities (Grazing and Fertilization) for the three study regions; Schwäbische-Alb (ALB), Hainich-Dün (HAI) and Schorfheide-Chorin (SCH). Abbreviations: fineRD = fine root diameter, LDMC = leaf dry matter content, Leaf N (or P) = leaf nitrogen (or phosphorus) concentration, RTD = root tissue density and SLA = specific leaf area, SRL = specific root length opened a lateral cubic tunnel (side of 15 cm) at a depth of 10 to 25 cm. We lined the floor of the tunnel with silica sand. SIAs were saturated with deionized water from below and then placed in the lateral tunnels, so that SIAs were installed under undisturbed soil. The empty space around and above the SIAs was filled with silica sand, and the pit was re-filled with the removed soil. The silica sand ensured that SIAs were not in direct contact with soil impeding, thus, direct exchange reactions between soils and resins and, at the same time, it allow the undisturbed penetration of water.
SIAs stayed in the soil for one year, until spring 2019. For the installation and removal, we visited the study regions in the same order. Installation and removal dates were recorded together with field observations (e.g. presence of roots and soil fauna) for a quality control of the samples. After removal, SIAs were stored at 4 °C, and each SIA was separated to three resin samples by separating the resin column into three resin layers. The first resin layer accounted for nutrient leaching, and had a depth of 5 cm. The second resin layer was a control, and the third layer served as a buffer for upward nutrient movement, and they had a length of 1 and 4 cm, respectively. Each resin sample was mixed to ensure homogeneity, weighed, and stored at 4 °C. To calculate leaching, we analysed the first two layers for their nutrient concentrations. If the second layer was contaminated with NO 3 (> 10% of first layer NO 3 -N concentration), then the SIA was excluded from our dataset. Out of the 600 SIAs, five SIAs were excluded due to disturbances by roots and/or soil fauna, and 34 SIAs due to a contaminated second layer.
Resin samples of 15 g were weighed. After adding 100 mL of 1 M NaCl (Natriumchlorid zur Analyse Acs, Merck KGaA, Darmstadt, Germany) the resin was shaken for 45 min. The extracted solution was filtered by gravity using N-free filters (MN 619 DE 1/4, Macherey-Nagel Gmbh & Co, Germany), and the filtrate was transferred to a 15-mL falcon tube, which was stored at − 20 °C. Ammonium (NH 4 -N), nitrite (NO 2 -N) and nitrite and nitrate (NO 2 -N + NO 3 -N) and phosphate (PO 4 -P) concentrations in the solution were determined by a flow-injection analyser (Quickchem QC85S5, Lachat Instruments, Hach Company, Dreieich, Germany) after applying appropriate dilutions. Nitrate concentration was calculated as the difference of [NO 2 -N + NO 3 -N] from NO 2 -N concentration. Sulphate (SO 4 -S) concentration was determined by an ion chromatography system (Dionex-DX 500, Thermo Fisher Scientific GmbH, Dreieich, Germany), after applying appropriate dilutions. Nutrient concentrations of the first and second layer were used to calculate leaching fluxes.
SIAs might underestimate the effects of grazing intensity on nutrient leaching compared to the effects of fertilization intensity due to the higher spatial heterogeneity of nutrient inputs through dung and urine patches compared to the more homogeneous nutrient inputs through fertilization. However, i) the high number of grazed sites included in our study (118 out of 150 permanent grasslands), ii) the use of average management history (fertilization and grazing intensities averaged over 2015-2018) instead of management applied only in 2018, and iii) the determination of leaching as an integrating measurement over a year, support the overall representation of management effects on nutrient leaching risk.

Data analysis
We tested for differences of annual nutrient leaching in the three study regions using analysis of variance (ANOVA) models and the Tukey's (honest significant difference) test. We use a backward stepwise function to assess the impact of (i) study regions, (ii) plant properties (aboveground biomass, species richness, legume cover and grass cover), (iii) soil OC, (iv) microbial biomass C and (v) grassland management (fertilization and grazing intensities) on nutrient leaching based on the Akaike Information Criterion. Using the variables selected by the stepwise function, we modelled nutrient leaching with analysis of covariance (ANCOVA). ANCOVA models were performed using a type II sum of squares, which is not influenced by the order in which the explanatory variables are introduced in the model (Zuur et al. 2009). For every model, we used diagnostic plots to test the assumptions of linearity, normality and homoscedacity of variance in the residuals, and to check for outliers. Heteroscedasticity was tackled with logarithmic or square-root transformations. Collinearities between explanatory variables were assessed with the variance inflation factor and values < 3.0 were accepted. We followed the same approach to separately investigate the importance of several physico-chemical soil properties (Table S3) and CWM traits on nutrient leaching (Table S4) listed in Table 2.
We used Structural Equation Modeling (Shipley, 2016) to investigate how fertilization and grazing, and their impacts on soils, plants and soil microorganisms, affect nutrient leaching. Grazing and fertilization intensities were included as exogenous variables. Plant species richness is an important predictor of soil functions (Chen et al. 2019;Leimer et al. 2016) and sensitive to management intensification (Midolo et al. 2019), and it was included as an endogenous variable. Substrate stoichiometry controls the cycling of C and nutrients in soils (Manzoni et al. 2008;Schleuss et al. 2021), and so we used stoichiometric ratios (C:N, C:P or C:S) of soils, plants and soil microbial biomass as endogenous variables in our models. We removed the effects of study regions (with one-way ANOVA models) from all endogenous variables, and we used their residuals in our Structural Equation Models. Proposed causal structures with p-values ≥ 0.05 were accepted. In addition, we calculated the root mean square errors of approximation (RMSEAs) and the standardized root mean square residuals (SRMRs), which both should be ≤ 0.05. The high number of observations compared with the number of estimated parameters in our models provided safety against multivariate non-normality issues. As an additional measure, we re-ran our models using bootstrapping. The p-values of bootstrapped models (not shown) were always similar to those of the classic models. We followed the same approach to separately investigate the importance element pools (C, N, P or S, Fig. S1) and plant CWM traits (selected based on Table S4) on nutrient leaching.
Structural Equation Modeling is among the few (semi)statistical methods claiming to prove, or disprove, causal relationships in a proposed structure. Natural ecosystems are characterized by complex interactions between processes involved in element cycling. This makes the development of a unique causal structure challenging. We recognise that results from our Structural Equation Models represent associations between variables, and not necessarily causal relationships. However, we interpret these associations as evidence of management effects and, for simplicity, we often use terms such as 'effects' and 'drivers' hereafter.
Statistical analysis was performed with the R statistical software (R Core Team 2019). Figures were designed with the ggplot2 package. ANOVA and ANCOVA models were performed with the lm function from stats package, and the backward stepwise elimination with the step function. The variance inflation factor was calculated with the vif function of the faraway package (Faraway, 2016). Tukey's test and PCA were performed with the TukeyHSD and the prcomp functions, respectively, of the stats package. The Structural Equation Modeling was performed with the sem function from the lavaan package (Rosseel, 2012).

Annual nutrient leaching from the topsoil in grassland
The median nitrate (annual) leaching in a depth of 10 cm was 7.4 kg NO 3 -N ha −1 (interquartile range: 2.2-22.2 kg NO 3 -N ha −1 , Table S1) over the three study regions and increased in the order ALB < SCH = HAI (Fig. 2a). The median leaching of the less mobile ammonium was 4.4 kg NH 4 -N ha −1 (interquartile range: 3.8-6.0 kg NH 4 -N ha −1 ) and increased in the order HAI < ALB = SCH (Fig. 2b). The cumulative inorganic N (NH 4 -N + NO 3 -N) leaching followed the patterns of nitrate leaching and increased in the order ALB < SCH = HAI, underlining the importance of nitrate for total N leaching losses. Interestingly, nitrate and ammonium leaching show comparable fluxes despite that nitrate is more mobile than ammonium. This might be due to the relatively shallow depth of our measurements since, within the rooting zone, nutrient inputs (e.g. from fertilization or decomposition) might be more important for leaching than nutrient mobility. The three study sites differ in atmospheric N inputs, which can be relevant for inorganic N leaching since, for instance, the annual bulk deposition of total N in SCH (11.8 ± 1.7 kg N ha −1 , Schwarz et al. 2014) is nearly four times higher than the annual fertilization inputs ( Table 2).
The median phosphate leaching flux was 0.8 kg PO 4 -P ha −1 (interquartile range: 0.6-1.7 kg PO 4 -P ha −1 , Table S1) over the three study regions and there were no detectable differences between the regions (Fig. 2c). Lastly, the median sulphate leaching flux was 21.8 kg SO 4 -S ha −1 (interquartile range: 13.7-29.5 kg SO 4 -S ha −1 ), and leaching in SCH was significantly higher than ALB (Fig. 2d), but not than HAI. The regional patterns of nutrient leaching differed from the patterns of soil N, P labile,inorg and total S concentrations for nitrate, phosphate and sulphate leaching, while ammonium leaching differences between the study regions matched those of soil N concentrations (Table 2).

Management effects on nutrient leaching
Including management intensities and selected soil, plant and microbial properties in our ANCOVA models for nitrate, ammonium and phosphate annual leaching (Table 3) considerably increased the explained variance compared to models including only soil properties (Table S3). Our models explained Fig. 2 Annual leaching of (a) nitrate (NO 3 -N), (b) ammonium (NH 4 -N), (c) phosphate (PO 4 -P) and (d) sulphate (SO 4 -S) in a depth of 10 cm for the three study regions; Schwäbische-Alb (ALB), Hainich-Dün (HAI) and Schorfheide-Chorin (SCH). Bars show average values and error bars represent standard deviations. Numbers at the base of the bars show the number of observations. Lower-case letters indicate significant differences (p-value < 0.05) among the study regions based on ANOVA and Tukey's HSD models 39, 27, 23 and 38% of the variance in nitrate, ammonium, phosphate and sulphate leaching, respectively (Table 3). Study regions were important predictors of leaching for every nutrient. Fertilization intensity increased nitrate, ammonium and phosphate annual leaching, while grazing had an increasing effect only on nitrate and phosphate leaching. Management was not important for sulphate leaching. Plant aboveground biomass was positively related to ammonium and phosphate leaching, but it was negatively linked to sulphate leaching. Plant species richness was negatively associated with nitrate leaching, and it was not important for the other nutrients. The relative cover of grasses was negatively related to nitrate leaching, and that of legumes was marginally (and positively) related to nitrate leaching. Soil OC was positively related to ammonium leaching. Finally, microbial biomass C was negatively associated with ammonium and phosphate leaching, but it was positively associated with sulphate leaching.

Management and stoichiometry effects on nutrient leaching
According to our structural equation models, fertilization intensity had (direct) positive effects on nitrate, ammonium and phosphate leaching risk, but not on sulphate, and these effects were not explained by plant species richness or soil, plant and microbial stoichiometry (Fig. 3). Similarly, grazing intensity had direct positive effects on nitrate and phosphate leaching, but not on ammonium or sulphate. Plant species richness was directly associated with ammonium and sulphate leaching, and these relationships were negative and positive, respectively. In contrast to our ANCOVA model (Table 3), plant species richness was not directly related to nitrate leaching (Fig. 3a). Interestingly, soil, plant and microbial C:N ratios were not associated with fertilization intensity directly, but rather indirectly through plant species losses caused by fertilization. In contrast, grazing intensity was not related to C:N ratios directly or indirectly. A higher plant species richness was associated with higher C:N ratios of soils, plants and microbial biomass. Both plant and microbial C:N ratios had a negative relationship with nitrate leaching, while soil and plant C:N ratios had a positive relationship with ammonium leaching. In a similar model with N pools instead of C:N ratios, microbial and plant N were Grey arrows indicate paths of marginal significance. Standardized path coefficients are given together with their significance level (*p < 0.050; **p < 0.010; ***p < 0.001). Coefficients of determination are given for the response variables positively linked to nitrate leaching, and plant N was negatively linked to ammonium (Fig. S1a, b).
Grazing and fertilization intensities were negatively associated with C:P ratios of plants and soils and with microbial biomass C (no microbial P data available, Fig. 3c). Apart from management, soil C:P ratio was the only parameter related to phosphate leaching, and this relationship was negative. In a similar model with P pools instead of C:P ratios, soil P labile,inorg was positively related to phosphate leaching, while plant P and microbial C were not important (Fig. S1c).
For S, no information on plants or microbes was available. Microbial C and plant species richness were (directly) positively related to sulphate leaching (Fig. 3d). However, plant species richness was negatively linked to microbial C, and, therefore, its effect on sulphate leaching was offset. Soil C:S ratio was not associated with fertilization or grazing intensities, and it was only marginally associated with sulphate leaching. In a similar model with soil S concentration instead of C:S ratio, soil S was positively related to sulphate leaching, and there were no other effects from management, plants or microbial biomass on sulphate leaching (Fig. S1d). The net effects of fertilization and grazing on sulphate leaching were negative and rather small (Table S5).

Effects of plant community traits on nutrient leaching
Based on backwards stepwise analyses and ANCOVA models, we identified CWM plant traits that were important for nutrient leaching (Table S4). Each nutrient had a different set of important CWM plant traits, while for sulphate leaching none of the traits was important. Our Structural Equation Models show that fertilization and grazing intensities exhibited multiple indirect effects on nutrient leaching via CWM plant traits, in addition to their direct effects (Fig. 4). Interestingly, while plant species richness was only affected by fertilization, CWM plant traits were, typically, influenced by both grazing and fertilization, though grazing acted more strongly on aboveground traits (with the exception of specific root length) than fertilization. Fertilization and grazing intensification seem to have shifted plant communities towards more resource acquisitive, fast-growing species with higher leaf N and P, and lower leaf dry matter content and root tissue density (Fig. 1). However, grazing was also negatively associated with specific root length.
Fine root diameter directly increased, and leaf dry matter content directly decreased, nitrate leaching (Fig. 4a). Microbial N mediated positive indirect effects from fine root diameter and leaf dry matter content on nitrate leaching. Interestingly, plant species richness directly decreased nitrate leaching, even after including selected CWM plant traits and microbial N in our model. In contrast, plant and microbial C:N ratios (from field measurements) fully mediated the effects of plant species richness on nitrate leaching (Fig. 3a), but its net effect was negative in both models (Table S5).
Root tissue density was marginally (and positively) linked to ammonium leaching, and leaf dry matter content was marginally (and negatively) linked to it (Fig. 4b, Table S4). Plant species richness was again negatively associated with ammonium leaching (as in Fig. 3b). Leaf P increased phosphate leaching, and it mediated indirect positive effects from both fertilization and grazing intensities on phosphate leaching (Fig. 4c). Specific root length was positively linked to sulphate leaching, and it mediated an indirect negative effect from grazing intensity on sulphate leaching (Fig. 4d). Finally, leaf N marginally (and negatively) related to sulphate leaching, and it mediated indirect negative effects by both fertilization and grazing intensities.

Discussion
Intensively managed grasslands have an elevated nutrient leaching risk Previous studies stressed the importance of grazing, and the associated dung and urine patches, for nitrate leaching (Ledgard et al. 2009), and they associated (grazed) pastures with higher nitrate leaching risk than (mowed) meadows (Cameron et al. 2013;Ryden et al. 1984). In this study, nitrate leaching risk significantly increased with both fertilization and grazing intensity (Table 3, Fig. 3). We observed that fertilized grasslands also have an elevated nitrate leaching risk, and we extended previous findings from mostly shortterm observations  to an annual timescale. In addition, ammonium annual leaching risk increased with fertilization intensity, but not with grazing (Table 3). This is possibly because grazing increases N inputs as urea, which is quickly oxidized to nitrate (in various soil types, Soulides and Clark, 1958;Tzanakakis et al. 2018), especially since management intensification increases nitrification rates (Stempfhuber et al. 2014). Considering both nitrate and ammonium, we conclude that total inorganic N leaching risk increases both with grazing and fertilization intensification across an annual timescale.
Agricultural activities are responsible for excess P losses via leaching, which contribute to the degradation of (surface) aquatic ecosystems due to eutrophication (SOER 2020). We found that managed grasslands can also have a high P leaching risk, since annual phosphate leaching at 10 cm soil depth were frequently greater than 0.8 kg PO 4 -P ha −1 (overall median) (Table S1, Fig. 2). Such fluxes are higher than leaching losses from the subsoil of other agricultural land-uses, like croplands (Leinweber et al. 1999). However, P leached from the topsoil, as measured in our study, can still be fixed in the subsoil, or be taken up by plants (Kuhlmann and Baumgärtel, 1991), and so it may not reach aquatic ecosystems. Therefore, our study can only identify an increased nutrient Grey arrows indicate paths of marginal significance. Standardized path coefficients are given together with their significance level (*p < 0.050; **p < 0.010; ***p < 0.001). Coefficients of determination are given for the response variables loss risk due to management intensification. However, the P sorption capacity in subsoils depends on soil texture and pH, and it could, at some point, be exceeded with long term excess P supply (Gocke et al. 2021;Werner et al. 2006), especially in sandy soils, such as those in SCH (Table 2). Overall, we found that phosphate leaching significantly increased with management intensity (Table 3), and that fertilization and grazing equally increased phosphate loss risk.
Compared to N and P, little is known for S cycling in temperate managed grasslands and, usually, S leaching is not measured as an actual flux but it is, instead, estimated based on soil solution S concentrations and water balance modeling. In this study, we measure sulphate leaching at 10 cm soil depth and the overall median was 21.8 kg SO 4 -S ha −1 (interquartile range: 13.7-29.5 kg SO 4 -S ha −1 , Table S1), which agrees with reported values from grassland and cropland subsoils (4-45 kg SO 4 -S ha −1 ; Erisken et al. 2000; Gallejones et al. 2012). Due to ample atmospheric S depositions in the past, temperate soils are typically considered to be S-rich. However, concerns that S might be a limiting resource for plant production and quality have been raised (Eriksen, 2009;Mathot et al. 2008). In this study, sulphate was the only nutrient for which management intensity reduced its leaching, and this operated indirectly (Fig. S1d). Both grazing and fertilization reduced soil S concentrations and, thus, sulphate leaching. A possible explanation is that management intensification increases plant productivity, which creates higher S demands and plant uptake (Eriksen 2009), eventually leading to increased export of S from the ecosystem with harvest or grazing. This is supported by a negative relationship between aboveground plant biomass and sulphate leaching in our grasslands (Table 3). Another explanation for management-dependent reductions in soil S is the desorption of sulphate from soil minerals due to phosphate additions, since phosphate is adsorbed more strongly than sulphate (Bolan et al. 1988). However, this process should have led to increased sulphate leaching risk with increased fertilization, which we did not observe. The lack of information on P and S fertilization rates might obscure additional effects from fertilization on sulphate leaching.
Fertilization, but not grazing, alters stoichiometry towards increased nutrient leaching risk The positive effect of plant N on nitrate leaching (Fig. S1a) indicates faster mineralization of N-rich litter under high management intensity. Meanwhile, the negative effect of plant N on ammonium leaching (Fig. S1b) could be similarly explained by higher nitrification rates in intensively managed grasslands (Stempfhuber et al. 2014) or, alternatively, by a preference for ammonium uptake by plants. In contrast to plant N, plant biomass, which increases with management intensification , was positively related to ammonium leaching (Table 3). This could be explained by increased above-and belowground litter inputs and, thus, increased soil organic N (Leimer et al. 2016), which is then mineralized to ammonium.
In contrast to ammonium, nitrate leaching risk was also controlled by microbial properties. Microbial nitrification, following soil organic N mineralization, could explain the positive relationship between microbial biomass N and nitrate leaching. Management intensity, and especially fertilization, was shown to increase the nitrification potential of the studied soils (Meyer et al. 2013;Stempfhuber et al. 2014). Therefore, fast nitrification of ammonium at fertilized sites with low C:N ratios could be the reason for the observed reduction of ammonium leaching with increasing microbial biomass C (Table 3) and for the observed increase of ammonium leaching risk with increasing soil C:N ratio (Fig. 3b).
Interestingly, microbial C:N ratios were negatively linked with nitrate leaching risk (Fig. 3a), indicating that leaching was not only affected by changes in plant communities, but also microbial communities (Fig. S2). This may be because the less nutrientdemanding, more fungal-dominated communities at extensively managed, botanically diverse grasslands capture N more effectively and prevent leaching (de Vries et al. 2006). Alternatively, it may be that the turnover and decomposition of fungal-derived, N-poor necromass is low. In contrast to the latter, evidence from short-term observations suggests that microbial communities which are not N-limited, increase N retention in the plant-soil system (de Vries and Bardgett 2016). To our knowledge, our work is among the few studies showing such a clear relationship between microbial biomass and annual nitrate leaching in managed grasslands under field conditions.
Microbial biomass N had a positive effect on nitrate leaching risk (Fig. S1a). In experimental grasslands, the negative effect of microbial biomass on the leaching of dissolved organic N has been attributed to enhanced mineralization (Leimer et al. 2016), which could explain higher soil nitrate concentrations and, thus, leaching. Previous work has shown that N additions can increase or decrease microbial biomass (Dietrich et al. 2017;Lu et al. 2011). Soil microbial communities are also influenced by plant communities, mainly through changes in N availability (Moreau et al. 2015), and perhaps to a greater extent than by management intensity in our study sites (Boeddinghaus et al. 2019). In our study, management intensity directly reduced microbial biomass N (and C), but the net effect was still close to neutral, because the fertilization-induced reduction in species richness increased it again (Fig. S1a). This suggests that relations between long-term fertilizer additions and microbial biomass in managed grasslands are modified by plant communities and that the observed responses of microbial activity to short-term fertilizer additions to soils in the laboratory, where plant communities are absent, might not be representative of field conditions. Phosphate leaching risk was mainly associated with soil P labile,inorg and C:P ratio, which were both positively related to plant P and C:P ratios (Fig. 3c, S1c). Accordingly, fertilization-induced reductions in plant species richness were positively associated with soil P labile,inorg and, thus, with phosphate leaching risk. In contrast to inorganic N, phosphate leaching was negatively related to soil pH (Table S3), indicating high P immobilization as calcium phosphate in alkaline soils.

Exploitative plant communities increase nutrient leaching risk in intensively managed grasslands
We observed that the functional group of grasses was negatively associated with annual nitrate leaching risk, while the group of legumes was positively (but marginally) associated with it (Table 3). Neither of them was as important for ammonium, phosphate and sulphate leaching risk. A positive relationship between legumes and nitrate leaching can be explained by increased N availability caused by their N-fixing ability and N-rich tissues (Leimer et al. 2015). Grasses decrease nitrate leaching due to increased N acquisition and/or more efficient use of N and water (de Vries and Bardgett 2016; Leimer et al. 2016) and due to their, typically, lower tissue N content. Indeed, we observed that legume cover was positively correlated with CWM leaf N, and that grass cover was negatively correlated with CWM fine root diameter (Fig. 1), which could, at least partly, explain the functional group effects on nutrient leaching.
Plant traits give hints on the importance of plant nutrient acquisition for nutrient leaching. In accordance with the observation that phosphate and sulphate leaching risk were more related to soil properties than inorganic N (Table S3), plant traits were less important for them, too (Table S4). We identified the CWM leaf dry matter content as an important predictor of inorganic N leaching (Fig. 4), reducing both, nitrate and ammonium leaching risk. This effect could be explained by a slower and more efficient N cycling in conservative than exploitative plant communities (Kleinebecker et al. 2014), and/or by lower N inputs in the soil system since leaf dry matter content was negatively related to i) fertilization, ii) legume cover and iii) CWM leaf N (Fig. 1). The positive relationship between CWM fine root diameter and nitrate leaching risk could be explained by a reduced nutrient uptake efficiency with coarser roots. For ammonium, leaching risk was related to root tissue density (Fig. 4b), where a higher root tissue density would indicate coarser fine roots (Fig. 1) and a reduced nutrient uptake potential.
To date, it is not clear whether conservative, slow-growing plant communities, or exploitative, fast-growing plant communities lead to a greater N retention in the soil, plant and soil microorganisms' system. For instance, conservative plant communities can lead to decreased N retention due to low N acquisition rates (de Vries and Bardgett, 2016), while exploitative plant communities can lead to decreased N retention due to enhanced N cycling and high N availability (de Vries et al. 2012;Grigulus et al. 2013;Robson et al. 2010). We observed a differentiation between the effects of above-and belowground traits on inorganic N leaching risk, similar to findings on N-related ecosystem functions by de Vries and Bardgett (2016). Conservative aboveground plant traits reduced inorganic N leaching risk, while conservative belowground traits increased it (Fig. 4). Regarding our third hypothesis, we suggest that management intensity, and especially fertilization, shifted mostly the aboveground plant traits towards values more representative of exploitative communities ( Fig. 1; Allan et al. 2015;Busch et al. 2018;Zobel et al. 2007), which are characterized by a fast tissue turnover and decomposition (Lavorel et al. 2011) and, thus, influenced their association with annual N leaching risk.

Plant diversity reduces inorganic N leaching risk
Previous studies have reported a negative relationship between plant diversity and nitrate availability or leaching in managed and experimental grasslands (Leimer et al. 2016;Niklaus et al. 2001), but effects on ammonium leaching are less often described (Leimer et al. 2015). Here, plant species richness was negatively associated both with nitrate and ammonium leaching risk, even after considering fertilization intensity (Table 3, Fig. 3b). Interestingly, the negative effect of plant species richness on nitrate leaching was completely mediated by higher plant and microbial C:N ratios (and lower N concentration) in more botanically diverse grasslands (Fig. 3a, S1a), while aboveground plant biomass was not selected as a predictor for nitrate leaching risk. In contrast to nitrate, the negative effect of plant species richness on ammonium leaching was not completely mediated by stoichiometry, element pools or CWM plant traits (Fig. 3b, S1b, 4b). However, it could be explained by the negative relationship between plant (aboveground) biomass and species richness in managed grasslands ) since plant biomass was positively associated with ammonium leaching (Table 3). Alternatively, it could be explained by niche complementarity and ammonium uptake by plants in species-rich grasslands, as shown for plant species identity and functional groups (Kahmen et al. 2006). In addition, the negative relationship between plant species richness and plant N led to a positive indirect effect on ammonium leaching risk (Fig. S1b). The negative relationship between plant species richness and plant N could be explained by fertilization that simultaneously causes plant biomass N enrichment and plant species losses, and/or by an increased plant N use efficiency with increasing plant diversity in grasslands (Dias et al. 2010;Fargione et al. 2007). To conclude, in terms of net effects, plant species richness was negatively associated both with nitrate and ammonium leaching risk (Table S5) and most of its effects on inorganic N leaching were mediated by the stoichiometry of soils, plants and soil microorganisms, as well as several CWM plant traits.

Conclusions
By applying a systems perspective, we identified that management intensity, and especially fertilization, not only promote nutrient leaching risk directly, but also indirectly by changing the properties of soils, plants and microbes. Most of the indirect effects of management on nutrient leaching were mediated by fertilization-induced reductions in plant species richness. While ceasing fertilization will reduce its direct effect on N and P leaching over a relatively short period in time, some of these changes, such as a loss in plant diversity and functional changes in plant communities, are not immediately reversible. Thus, management intensification changes the risk of nutrient leaching on the long-term and active (restoration) measures to maintain diverse plant communities are needed, which might, in return, mitigate the increased leaching risk that management intensification imposes on grasslands. Finally, we showed that, when considering direct and indirect effects, both fertilization and grazing lead to higher nitrate and phosphate leaching risk in the topsoil, while only fertilization intensity associates with higher ammonium leaching risk. As such, shifting grassland management towards (mowed) meadows instead of (grazed) pastures should not be considered as a management strategy for N leaching mitigation, as long as meadows are fertilized.
We thank the administration of the Hainich national park, the UNESCO Biosphere Reserve Swabian Alb and the UNESCO Biosphere Reserve Schorfheide-Chorin as well as all land owners for the excellent collaboration. Field work permits were issued by the responsible state environmental offices of Baden-Württemberg, Thüringen, and Brandenburg. This work has been funded by the DFG Priority Program 1374 "Biodiversity-Exploratories'' (SCHR 1181/2-3) and the Max Planck Society.
Authors' contribution Antonios Apostolakis designed and carried out the field and lab work related to nutrient leaching, analysed and discussed the data, and was responsible for writing this manuscript. Ingo Schöning, Beate Michalzik and Marion Schrumpf designed the experiment, discussed the data, and reviewed drafts of the manuscript and IS carried out also part of the field work. Valentin H. Klaus and Peter Manning discussed the data and reviewed drafts of the manuscript. Wolf-Anno Bischoff and Andreas Schwarz provided expertise on the resin method used for the nutrient leaching measurements. Yvonne Oelmann and Elisabeth Schurig contributed soil (Olsen) P data. Runa S. Boeddinghaus, Ellen Kandeler and Sven Marhan contributed data related to microbial properties. Valentin H. Klaus, Ralph Bolliger, Norbert Hölzel, Till Kleinebecker, Daniel Prati and Markus Fischer contributed data related to plant properties. Margot Neyret, Peter Manning, and Mark van Kleunen contributed above-and belowground plant (community weighted mean) trait data. All authors contributed to the manuscript and approved the final version.
Funding Open Access funding enabled and organized by Projekt DEAL. This work has been funded by the DFG Priority Program 1374 "Biodiversity-Exploratories'' (SCHR 1181/2-3) and the Max Planck Society.

Availability of data:
This work is based on data elaborated by several projects of the Biodiversity Exploratories project. Part of the data used are available at database of the project (https:// www. bexis. uni-jena. de/ ddm/ publi csear ch/ index). However, to respect the rights of PhD students of individual projects, the remaining data will be made publicly available at the same website after an embargo period of 3 years, according to the data and publication policy of the Biodiversity Exploratories project.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.