Assessing multifunctionality of agricultural soils: Reducing the biodiversity trade‐off

Soils are indispensable for the provision of several functions. Agricultural intensification and its focus on increasing primary productivity (PP) poses a threat to soil quality, due to increases in nutrient loads, greenhouse gas emissions and declining biodiversity. The EU Horizon 2020 Landmark project has developed multi‐criteria decision models to assess five soil functions: PP, nutrient cycling (NC), soil biodiversity and habitat provision (B‐HP), climate mitigation and water regulation, simultaneously in agricultural fields. Using these algorithms, we evaluated the supply of PP, NC and B‐HP of 31 grasslands and 21 croplands as low, medium or high. The multi‐criteria decision models showed that 38% of the farms had a medium to high supply of all three soil functions, whereas only one cropland had a high supply for all three. Forty‐eight per cent of the farms were characterized by a high supply of PP and NC. We observed a clear trade‐off between these two functions and B‐HP. Multivariate statistical analyses indicated that higher organic inputs combined with a lower mineral fertilization concur with higher biodiversity scores while maintaining a medium delivery of PP and NC. Additionally, we compared the outputs of the model predictions to independent variables that served as proxies for the soil functions and found: (a) croplands (but not grasslands) with high PP had a higher standardized yield than those with medium PP; (b) grasslands (but not croplands) with high NC had a significantly lower fungal to bacterial biomass ratio, suggesting faster decomposition channels; and (c) a positive though non‐significant trend between B‐HP score and rank according to soil invertebrate biodiversity. These comparisons suggest a successful upscaling of the models from field to farm level. Our study highlights the need for systematic collection of management‐related data for the assessment of soil functions. Multifunctionality can be achieved in agricultural soils; however, without specifically managing for it, biodiversity might come at a loss.


Highlights:
• We study how well soils can provide primary productivity, nutrient cycling and biodiversity.
• We study trade-offs and synergies among soil functions, as well as the drivers of these relationships.
• Soil biodiversity is largely sacrificed for primary productivity and nutrient cycling • Changes in pesticide and fertilizer management can increase soil multifunctionality.

K E Y W O R D S
cropland, grassland, land management, nutrient cycling, primary productivity, soil biodiversity, soil functions, yield

| INTRODUCTION
Humans derive multiple benefits from the soil system. Soils provide us with feed, fibre and fuel, cycle and mobilize nutrients, control diseases, regulate the climate, maintain water quality and provide a habitat to support biodiversity, amongst other things (Dominati, Patterson, & Mackay, 2010). The continued delivery of these soil functions is at risk, due to the effects of economic growth, population growth and climate change (Montanarella et al., 2016). These threats result not only in the loss of functional soils, but also reduce the capacity of the soils to function optimally. To face this challenge we need to change land management and extend our focus away from the delivery of only primary productivity (PP; Clay, Garnett, & Lorimer, 2019;Rutgers et al., 2012;Schulte et al., 2014;Tilman, Balzer, Hill, & Befort, 2011).
Agricultural soils are particularly complex to manage. Farmers face a continuous challenge to maintain high productivity whilst considering the resilience of their soil to ensure that the land can support continued productivity for future generations. Additionally, maintaining or even increasing agricultural yield is imperative to sustain a growing population, but to do so, farmers have relied heavily on external inputs and intensive crop rotations, which result in increased greenhouse gas emissions, eutrophication, pollution, loss of organic carbon and loss of above-and belowground biodiversity (Foley et al., 2005;Foley et al., 2011;Matson, Parton, Power, & Swift, 1997;Tilman et al., 2001;Tsiafouli et al., 2015). These effects clash with societal goals regarding nature conservation and climate change (e.g., the Common European Agricultural Policy, which intends to contribute to sustainability and climate change mitigation; https://ec.europa.eu/info/food-farming-fisheries/key-policies/ common-agricultural-policy).
Since the introduction of the ecosystem service concept (Costanza et al., 1997), many studies have attempted to quantify ecosystem services, but only a few have studied several simultaneously, and even less have explored the synergies and trade-offs therein (Seppelt, Dormann, Eppink, Lautenbach, & Schmidt, 2011). Understanding the trade-offs and synergies between soil functions as well as recognizing which variables might be important in swaying the system, is a step forward in providing management solutions to farmers and policymakers alike (O'Sullivan et al., 2015). Schulte et al. (2014) proposed such a framework (the functional land management framework) to assess the supply of soil functions in order to simultaneously meet agronomic and environmental objectives. It is based upon the principle that the suitability of soils to provide a given function, for example PP, depends on soil composition and environmental conditions, as well as the choices made by the land manager, and it has steered the development of tools to aid land management, from farmers' decision making (Soil Navigator; http://soilnavigator.eu) to European policy assessment (Schulte et al., 2019) and monitoring (Zwetsloot et al., this issue).
The functional land management framework allows for the simultaneous assessment of five predominant soil functions related to agricultural soils . These are: (a) PP, defined as the capacity of a soil to supply nutrients and water and to produce plant biomass for human use, providing food, feed, fibre and fuel ; (b) water purification and regulation, as the capacity of a soil to remove harmful compounds from the water that it holds and to receive, store and conduct water for subsequent use and the prevention of both prolonged droughts and flooding and erosion; (c) climate mitigation, as the capacity of a soil to reduce the negative impact of increased greenhouse gas emissions on the climate; (d) nutrient cycling (NC), as the capacity of a soil to receive nutrients in the form of by-products, to provide nutrients from intrinsic resources or to support the acquisition of nutrients from air or water, and to effectively carry over these nutrients into harvested crops (Schröder et al., 2016;Schulte et al., 2015); and (e) soil biodiversity and habitat provision (B-HP), as the multitude of soil organisms and processes, interacting in an ecosystem, providing society with a rich biodiversity source and contributing to a habitat for aboveground organisms .
In this study, we evaluated the performance of three out of the five soil functions (PP, NC and B-HP) in 52 Dutch farms using the multi-criteria decision models (MCDMs) developed by Debeljak et al. (2019) after adjustment for use at the farm level. We explore the synergies and trade-offs between these three functions as well as the drivers of high functionality. The following research questions were answered: 1. Were the soils under study capable of delivering a medium to high supply of more than one function simultaneously? 2. Do synergies and trade-offs exist between the soil functions? 3. Which input variables in the model lead to high supply of soil functions and soil multifunctionality?

| Assessing soil functions
To answer all the research questions, we first calculated soil function scores. In order to do this, we first (a) gathered the necessary data to calculate soil functions using (b) a decision support system. Because the data we obtained were mostly at the farm level, but the decision support system was designed to be used at the field level, we (c) had to adjust the model from field to farm level. We checked the success of this adjustment to the model by selecting proxy indicators of each soil function and comparing them to the scores obtained using the adjusted decision support system. All statistical analyses were carried out using R version 3.6.1 (R Core Team, 2017).

| Data acquisition
The first step to answer all research questions was gathering farm data on soil quality and management. Throughout the years 1993 and 2014, biotic and abiotic measurements of soil quality were assessed as part of the Netherlands Soil Monitoring Network (Rutgers et al., 2009a). The assessment included new sites over the years, covering a range of agricultural and natural ecosystems. Simultaneously, data were collected by Wageningen Economic Research from commercial farms through the Minerals Policy Monitoring Programme (RIVM 2019(RIVM -0026, 2019, aimed at monitoring the effectiveness of the Nitrates Directive (EU, 1991) and the National Minerals Policy (MANFQ & MIWM, 2017). These data provide data on farm nutrient balances, including information on the type and amount of manure and mineral fertilization as well as the yield and area cover for each crop (Tables S1 and S2). We used these sources of information to find farms that had been sampled by both institutions in the same year, leading to a total of 52 farms: 31 dairy farms (referred to as grasslands from here onwards, although they often reserve some land for maize) and 21 arable farms (or croplands) sampled between 2006 and 2014. Weather variables (Table  S3) were extracted from the weather station closest to each farm and are specific to the year that samples were taken from said farm using data from the Royal Dutch Meteorological Agency. Bulk density, cation exchange capacity, groundwater depth and salinity were not measured in the field but were calculated using published transfer functions or extracted from local maps (Tables  S2 and S3 ;De Vries, 2007;Dufour, 2000;Helling, Chesters, & Corey, 1964;Rawls, 1983;RIZA, 2004).

| Decision support system
In order to calculate the performance of soil functions in the selected farms, we used MCDMs designed by the Landmark project as part of the EU Horizon 2020 . The initial models were developed by experts in each of the soil functions and were further evaluated and validated (Sanden et al., 2019;Van de Broek et al., 2019;. The experts established a hierarchical model that builds up to a final function score (low, medium or high supply) starting from qualitative input variables related to the management, environment and soil properties. Discretization of these variables is carried out prior to running the model in accordance with a set of thresholds suited to different land uses and climatic zones. The input variables aggregate to the upper level of the model structure through if-then functions and aggregation continues until the top node is reached and an assessment of the performance of a specific function is made. For example, in order to assess B-HP (top node), the model is fitted with input variables such as "Soil pH" and "Liming" at the lowest nodes. These two variables aggregate through if-then functions: if they are both low, then "pH condition" is also low. "pH condition" when aggregated with "Nutrient content" and "Nutrient input" (which each are assessed by their own input variables) leads to a score in the super-attribute "Nutrients". Lastly four super-attributes (nutrients, biology, structure and hydrology) are aggregated to a final B-HP score. The Decision Expert integrative methodology was used to build these decision models (Bohanec, 2008;Bohanec & Rajkovič, 1990;Bohanec, Žnidaršič, Rajkovič, Bratko, & Zupan, 2013). The combination of the two datasets did not include all of the variables necessary to compute the five main soil functions, even after we complemented the dataset by analysing dried soil samples taken by the Netherlands Soil Monitoring Network (Rutgers et al., 2009b) for total nitrogen, carbon to nitrogen ratio, pH, organic carbon (%), organic matter content (%), clay content (%), phosphorus and magnesium (Tables S2 and S3). We therefore chose to focus on PP, NC and B-HP. Additionally, we had to make some assumptions about the farms. Slope was considered low (<3 ) for all farms. Because organic arable farms did not use chemical pesticides, we assumed they used some form of biological pest management. Crop failure risk, defined as the number of times in the past 20 years that a field provided no yield in the MCDMs (http://landmark2020.eu/landmark-glossary/), was assumed to be low in all farms. In the Netherlands, drought can have a strong effect on crop yield, but it rarely leads to complete crop failure. A study into the most extreme droughts since 1970 found that the most affected crop (onions) can experience up to a 50% decrease in yield when compared to years with average precipitation. Such a drop in yield was documented only once since 1976 (Prins, Jager, Stokkers, & van Asseldonk, 2018). Soil phosphorus status was assessed according to the Dutch Nitrates Directive (2014-2017; Table S3). Due to a number of missing variables while calculating PP and NC (biological-pest management, rooting depth, grass cutting frequency and share of legumes), some of the farms did not receive a full score, but rather were classified between two performance scores. In such cases, we downgraded the farm to the lower of the two suggested scores. In total, this measure affected three grasslands, which scored medium to high in NC and were downgraded to medium NC; one cropland that was reclassified as having medium PP; and all the grasslands performing at medium PP, which were initially performing at medium to high PP.

| From field to farm level
The MCDMs were intended to be used at the field scale; however, the data within the Netherlands Soil Monitoring Network on soil quality as well as some of the information collected in the Minerals Policy Monitoring Programme, were collected at the farm level. Consequently, parts of the models had to be adjusted to account for farm-level variability. The input variable "number of crops in rotation" was changed to "crop diversity" to reflect the spatial variability in crops within a farm using the Shannon-Weaver index, calculated using the "diversity" function of the "vegan" package for R (Oksanen et al., 2019). The threshold values for low or high crop diversity were the lower and higher quartiles of the calculated indices: 1.10 and 1.52, respectively, with an average crop diversity of 1.35. The input variable "loss of ammonia" from manure per farm was calculated within the Minerals Policy Monitoring Programme (RIVM 2019-0026, 2019; Appendix 2), and as with crop diversity, the thresholds for low, medium and high were derived from the distribution of the data, such that losses below 11 kg of nitrogen per hectare were considered low, between 11 and 30.3 kg they were considered medium, and above 30.3 kg they were considered high.
To assess the success of upscaling the MCDMs from field to farm level, we compared the results obtained from the farm-level MCDMs with selected proxies for soil functionality. Choosing appropriate proxies was complicated because the chosen variables had to be independent from those required by the MCDMs (Tables S2 and S3) and the MCDMs encompass many aspects of each soil function.
We compared the results from the PP MCDM with standardized yield. The Minerals Policy Monitoring Programme contains information on the crop, the area devoted to said crop, and the crop's yield for each farm. Using this information, we calculated standardized yield using the following formula: where n is the number of crops in farm i, MY i is the measured yield for crop i, PY i is the countrywide potential yield for that crop (Silva et al., 2020) and PA i is the proportional area that that crop occupies in the farm. This approach had three limitations: (i) potential yield was only known for the most common crops grown in the Netherlands, which means fields with uncommon crops were not included in the calculation; (ii) the grass yield was not measured, but rather calculated from information on the farm's overall performance using the amount of external feed bought, the number of cows in each dairy farm, the energetic requirements of cows, and the total milk production of the farm (RIVM 2019-0026, 2019; Appendix 2); and (iii) there was no potential yield estimate for grass, therefore, we used the distribution of grass yield to establish the low, medium and high categories, using the first and third quartiles as the threshold values. Although potential yield is likely to be different for different soil types, we chose to use a national level potential yield to be able to compare the standardized yield with the result from the PP MCDM, which includes soil type in the calculation of the PP score.
Due to the fact that there are two distinct decomposition paths in the soil (fungal and bacterial) we compared the NC MCDM scores with two biological indicators: the fungal to bacterial (F:B) biomass ratio and the ratio between the maturity of the plant parasitic nematode community (PPI; Bongers, 1990) and the maturity of the free-living nematode community (MI; Bongers, 1990). Sites with readily decomposable materials tend to be bacterial dominated (Bloem & Breure, 2003). Nematode composition can also provide insight into the efficiency of the nutrient cycle. Bongers, van der Meulen, and Korthals (1997), based on empirical observations, proposed that the PPI:MI ratio can indicate a nutrient enrichment effect. With increasing fertilization, an increase in opportunistic species causes the MI to decrease and the PPI to increase due to an increase in carrying capacity for plant parasitic nematodes. The microbial and nematode communities were sampled as part of the Netherlands Soil Monitoring Network. The MI and PPI were calculated using the online tool Ninja (Sieriebriennikov, Ferris, and de Goede (2014); accessed on May 2019). Neither of the two selected proxy indicators (PPI:MI ratio and F:B biomass ratio) can fully compare to the full NC model because both proxy indicators can only represent variations in the capacity of a soil to mineralize materials ("Mineralization" in the NC MCDM) and, to a certain degree, the ability of a soil to make nutrients available to plants ("Nutrient recovery" in the NC MCDM), but neither can capture the variation within the Harvest index, the third upper node in the NC model.
Lastly, we compared the results from the B-HP model with the farms' biodiversity rank, calculated by using principal component analyses (PCAs; Legendre & Legendre, 2012;Ter Braak, 1986). The PCAs explained the variation in the invertebrate diversity, including earthworm, micro-arthropod and enchytraeid species and nematode genera. As input variables for the PCA we used the Shannon-Weaver diversity, H 0 = P n i = 1 p i × lnðp i ), where p i is the proportional abundance of species i, and Pielou's evenness indices, where H 0 max is the Shannon-Weaver diversity if all species present in a site were present with equal proportional abundance (Pielou, 1966;Shannon, 1948) of the aforementioned faunal groups, scaling the input variables to unit variance.
We performed either one-way ANOVAs or a Kruskal-Wallis test (when the assumptions of homogeneity of variance and normality of the residuals were not met) using each of the selected indicators as the dependent variable and the MCDM scores for each soil function as the grouping variable (Chambers & Hastie, 1992;Hollander & Wolfe, 1999). PCA was calculated with the "rda" function in the "vegan" package (Oksanen et al., 2019).

| Soil multifunctionality
To answer our first research question, we created frequency plots for the low, medium and high scores in each soil function (PP, NC and B-HP) and assessed the number of farms with a high and low supply of all three, two, one or none of the studied functions.

| Synergies and trade-offs
To understand the synergies and trade-offs that occurred between PP, NC and B-HP, we used frequency plots and observed the overlap between sites with low, medium or high supply scores. We also studied how often functions tended to be supplied at a high or low level together.

| Variables associated with soil multifunctionality
Studying the variables that associate with multifunctionality was a two-step process. Firstly, we looked into the scores of the highest tier of each of the MCDMs. Each of these super-attributes is a result of several aggregations of discretized input variables into several layers of tiers. The upper tiers represent distinct parts of each soil function. They are "Soil", "Management", "Land userelated management", and "Environment" for PP; "Nutrient recovery", "Mineralization" and "Harvest index" for NC; and "Structure", "Nutrients", "Hydrology" and "Biology" for B-HP. The distribution of high, medium and low scores for these super-attributes can help us understand the drivers behind the MCDM scores. Secondly, we performed redundancy analyses (RDA; Legendre & Legendre, 2012;Ter Braak, 1986) to understand which variables better explained the capacity of grasslands and croplands to provide high PP, NC and B-HP. We calculated six separate RDAs, specifically three for grasslands and three for croplands, using either the soil, the environmental or the management-related variables that served as inputs for the MCDMs. To reduce the number of explanatory variables and avoid over-parameterization, we eliminated redundant variables, variables with very little variation and those that were co-linear when they had a variation inflation factor higher than 10 ( Greenacre, 1984;Gross, 2003). Cation exchange capacity and bulk density were not independent from soil organic matter and clay content; hence, these parameters were not included in the RDA models. Yearly average temperature and precipitation correlated strongly with the number of days with an average temperature above 5 C and precipitation during the first month of the growing season, respectively. Only the latter two were included into the RDA models. Soil pH, phosphorus and potassium content had high variation inflation factors and were excluded from the RDA model. We used the "ordistep" function for variable selection. We calculated the adjusted R 2 of each model using the Peres-Neto, Legendre, Dray, and Borcard (2006) permutation approach and an ANOVAlike permutation test for redundancy analysis of principal coordinates for each of the covariates, axes and for the whole RDA model (Legendre, Oksanen, & ter Braak, 2011). All calculations were performed using the "vegan" package for R (Oksanen et al., 2019).

| From field to farm level
There was a significantly lower standardized yield in croplands with medium PP than in those with high PP (Kruskal-Wallis χ 2 = 4.9, p-value = .03). For the grassland farms we found no significant differences between the standardized yield for farms categorized into the low, medium or high supply of PP (ANOVA, F = 2.13, p-value = .15).
We observed no significant differences in the PPI:MI ratio of grasslands (Kruskal-Wallis χ 2 = 0.39, pvalue = .22, with a mean PPI:MI of 1.54) or croplands (Kruskal Wallis χ 2 = 1.79, p-value = .23, with a mean PPI:MI ratio of 1.57) with low, medium or high NC. The F:B ratio was significantly lower in grasslands with high NC compared to those farms with medium NC, but neither differed significantly from the low NC supply (Kruskal Wallis χ 2 = 8.31, p-value = .02). There was no significant relationship in the croplands, where overall the F:B ratio was very low (Kruskal Wallis χ 2 = 2.95, pvalue = .23, with an average value of .56).
The first two axes of the PCAs explaining the variation in diversity indicators in grasslands and croplands explained 54% and 56% of the total variation for grassland and cropland farms, respectively ( Figure S1). In both cases, enchytraeid evenness was strongly positively correlated with the first axis and nematode diversity with the second, as can be observed from the direction of increase of both indices along the axes of the biplot ( Figure S1).
General positive trends were apparent between the scores of the first PCA axes and the B-HP scores, but none of the differences were statistically significant (ANOVA, F = 2.34, p-value = .12 for grasslands and ANOVA, F = 1.57, p-value = .23 for the croplands). The rank of the sites along the second RDA was also not significantly related to the score received by the MCDM (ANOVA, F = 0.04, p-value = .96 for grasslands and ANOVA, F = 0.95, p-value = .41 for the croplands).

| Soil multifunctionality
Only two grasslands performed poorly in two functions simultaneously (Figure 1c), in other words from the 31 grasslands included in our study, 29 of them delivered two or more functions at a medium or high level simultaneously, whereas all the croplands delivered two or more functions (Figure 1d). Additionally, none of the farms had a low supply of all three functions (Figure 1c,d). The MCDM, however, predicted that only one cropland soil could score high in all three functions (Figure 1b). PP was high in 18 grasslands and 18 croplands (58% and 86%, respectively; Figure 1a,b), and only three out of the 52 farms (6%) had a low supply of PP. Seven out of the 52 sites (13%) showed a low NC, whereas 46% of the farms had a low B-HP (Figure 1c,d). The proportion of farms with low, medium and high NC was similar in both grasslands and croplands (Figure 2a,b).

| Synergies and trade-offs
As the MCDM model provided the supply of the three functions for any given farm, it enabled further assessment of potential synergies and trade-offs between the soil functions. High PP scores co-occurred with high NC in 12 grasslands (39%) and 13 croplands (62%), but only one grassland and three croplands had both high B-HP and high PP simultaneously (Figure 1a,b). Twenty farms did not receive a low score for any of the functions and sites with medium PP did not necessarily get a medium B-HP score, which implies that the relationship between the provision of PP and B-HP is not strictly linear (Figure 3).

| Variables associated with soil multifunctionality
The analysis of the top tier of the MCDM (or super-attributes) that adds up to the PP score gave us insight into which parts of the soil system might be driving the score.
Out of the four super-attributes related to PP, "Environment" always scored high (Figure 4a,b), because the input variables related to "Environment" (slope, annual precipitation and annual temperature) were considered suitable for high PP by the MCDM. "Soil", which is an aggregate of input variables such as pH, soil organic matter, carbon to nitrogen ratio, abundance of macroelements, cation exchange capacity or salinity, was often either medium or low (Figure 4a,b). "Grassland management" (which results from the evaluation of variables such as stocking rate and legume presence; Figure 4a) showed a lower score than "cropland management" (crop diversity, legume cover and type of crops; Figure 4b).
The super-attributes determining the NC supply all indicate mostly medium to high supply capacity (Figure 4c,d). The "Harvest index" was high on all farms, which is the result of a combination of an assumption of low residues left in the field and low crop failure risk (both input variables into the MCDM). "Mineralization", which in the MCDM is a result of the aggregation of temperature, precipitation and pH, was medium on only 13% of the farms, and for the remainder was either high or low. No croplands scored low in terms of nutrient recovery (Figure 4d).
According to the scores of the B-HP super-attributes, the soil "Structure" is more suited to B-HP in croplands than grasslands. "Nutrient status" and "Hydrology" of the F I G U R E 1 Number of grasslands (left panels) and croplands (right panels) with a high (a,b) or low (c,d) supply of none, one, two or three soil functions. Scores were obtained from multi-criteria decision models for primary productivity, nutrient cycling and biodiversity and habitat provision. Colours denote which soil functions/s had a high or low score. B-HP, biodiversity and habitat provision; NC, nutrient cycling; PP, primary productivity F I G U R E 2 Frequency of functionality scores derived from multicriteria decision models for primary productivity (grey), nutrient cycling (pink) and biodiversity and habitat provision (green) on grasslands (a) or croplands (b) in the Netherlands studied soils, however, show a higher potential for high B-HP in grasslands than croplands (Figure 4e,f). The "Biology" super-attribute, however, received either a medium or a low score in all sites, except one cropland (Figure 4f). The variables that determine the "Biology" super-attribute were the use of pesticides, the presence of grassland in the arable rotation, crop diversity and legume presence, variables that were low in most of the studied farms.
According to the results from the RDAs, management variables ( Figure 5a) and environmental variables (Figure 5c) explained an important part of the variation in the distribution of high supply scores for different soil functions in grasslands (with an adjusted R 2 of 0.52 and 0.28, respectively; Table 1). Grasslands with a high capacity to support NC and PP tended to have a low groundwater table and high clay content (above 25%; Figure 5c) and were measured in years where the number of days with an average temperature above 5 C was medium (i.e., 180-230 days). In grasslands, clay content was positively correlated with magnesium content (Pearson ρ = 0.7). According to the RDA, these sites also show a positive correlation with high mineral nitrogen fertilization (more than 150 kg N/ha per year) combined with cows spending a high amount of time in the fields. High mineral nitrogen fertilization showed collinearity with the use of pesticides and only the first was included in the models. In contrast, grasslands with a high capacity to provide B-HP tended to have a high groundwater table and were sampled in years with a higher number of days with an average temperature above 5 C (above 230 days; Figure 5c) as well as a mineral nitrogen fertilization below 50 kg N/ha per year (Figure 5a).
For croplands, important parts of the variation in the high supply of PP, NC and B-HP were explained by soil and management variables (Table 1). Croplands with high PP and NC scores were characterized by very low soil organic matter (below 2%; Figure 5d), and high levels of mineral nitrogen fertilization (>150 kg N/ha per year) combined with low levels of organic nitrogen fertilization (Figure 5b; <50 kg N/ha per year). Soil organic matter was strongly and negatively correlated with bulk density (Pearson, ρ = −0.75), but high fertilization rates were not strongly correlated with any of the variables not included in the RDAs, that is with a correlation coefficient equal to or larger than 0.7. High B-HP scores were correlated with a medium to high nitrogen to phosphorus ratio (above 20) and a high carbon to nitrogen ratio, and more than 8% soil organic matter, but correlated negatively with the use of cover crops and crop diversity. The nitrogen to phosphorus ratio was negatively correlated with pH (Pearson, ρ = −0.71), meaning that this observation could in fact be a result of lower pH values in sites with high B-HP than in those with a low score. Similarly, the carbon to nitrogen ratio of croplands was negatively correlated with clay content, bulk density and pH (Pearson, ρ = −0.73, ρ = −0.74 and ρ = −0.97, respectively) and positively with the phosphorus and potassium load (Pearson, ρ = 0.82 and ρ = 0.73, respectively).

| DISCUSSION
We upscaled the MCDMs developed by the Landmark project from field to farm scale to assess the soil's ability to provide PP, NC and B-HP in 31 grasslands and 21 croplands in the Netherlands. The soils often delivered more than one function at a time with a medium or high score, but it was often a combination of PP and NC. B-HP received a low score on 24 out of 53 farms (46%), whereas NC and PP were low on seven (13%) and three (6%) of the farms, respectively. Management variables, namely the source of nitrogen fertilization and the type of crops in rotation, play an important role in explaining the trade-offs and synergies between B-HP, NC and PP.

| Multi-criteria decision models: The issue of data collection
The MCDMs were originally created as a tool for farmers or advisors to assess the multifunctionality of their soils.

F I G U R E 3 Frequency of grasslands (left panels) and
croplands (right panels) with low (L), medium (M) or high (H) capacity to provide primary productivity and nutrient cycling (a,b) and primary productivity and biodiversity and habitat provision (c, d). Soil function scores were obtained using multi-criteria decision models As such, they were designed to utilize soil management information alongside basic soil property and environment data to derive the multifunctionality of soils at field scale. These models have provided a framework for the assessment of soil functions across Europe and they are being utilized to support scientific assessment of multifunctionality, but in this context the models often suffer from the lack of management data, which are essential for each model (Schröder et al., 2016;. This proved a disadvantage in our research. The soil data and environment data were readily available at a farm level, and when not directly measured we could access soil maps that helped fill in the gaps. But overall, few existing research projects or organizations have field or farm-level information that includes details about soil/land management. This resulted in 20% of the input variables, particularly management-related variable, being based on either proxy values derived from literature research (such as the depth of the organic matter layer) or on expert assumptions of common agricultural practices in the Netherlands (such as a low amount of crop residues left in the fields). We accept that these assumptions could introduce a potential error into the calculation of the final function score within the model and reduce the precision of the predictions, but the final effect of this on our results is difficult to assess without comparing them to results gathered with a full set of input variables. More importantly, because not one single variable is solely responsible for the score obtained but rather a result of a tiered aggregation process, we do not expect the resulting score to be inaccurate. We suggest that future studies that use these MCDMs carry out a sensitivity analysis. The importance of soil management variables within this study emphasizes the need for future F I G U R E 4 Frequency of scores for the top tier of multi-criteria decision models for primary productivity (grey), nutrient cycling (pink) and biodiversity and habitat provision (green) on grasslands (a) or croplands (b) in the Netherlands. Each tier represents the suitability of part of the system to sustain a high, medium or low supply of each function. ''Soil'' stands for soil properties in primary productivity. In biodiversity and habitat provision "structure" stands for soil structure and "nutrients" for nutrient status [Color figure can be viewed at wileyonlinelibrary.com] research projects to collect data additional to soil-based measurements to really be able to define the functional capacity of a soil. We suggest, specifically, that studies exploring soil quality in European systems collect at least the information on farm management collected by the Landmark project, or a variation of this that includes more categories for each variable .

| From field to farm level: Challenges and opportunities
Comparison of the MCDM with proxy variables was performed to see if the upscaling from the field to farm level had been successful. We did not, however, expect a perfect relationship between the proxy variables and the functioning scores according to the MCDM. The proxy variables serve as a representation of one part of a soil function, whereas the MCDMs attempt to encompass all the main processes that affect the same soil function.
We observed a clear difference between the average standardized yield observed in croplands with medium and high PP scores. These results are in line with those obtained during the model's validation process . However, comparison of the standardized yield with the functional capacity of PP in the grassland sites did not result in a significant relationship. The grassland PP model had not previously been validated (unlike F I G U R E 5 Redundancy analyses of the effect of management-(a and b), environmental-(c) and soil-related variables (d) on the high (H) capacity of the soil to provide primary productivity (PP), nutrient cycling (NC) and biodiversity and habitat provision (B-HP) in grasslands (left panel) and croplands (right panels). Circles denote the position of the farms and the colour the number of farms in that position. Only RDAs where variables (indicated by crosses) explained a significant variation in distribution of the farms' high scores have been included in this figure, leaving out soil variables in grasslands and environmental variables in croplands. SOM stands for soil organic matter and GWT for ground water table [Color figure can be viewed at wileyonlinelibrary.com] the PP crop model) and might require adjustment to better represent potential PP. The processes underlying the PP MCDM, however, are very similar in grasslands and croplands and should represent the relationship between soil properties and PP even after being adjusted to the farm level. In our dataset, an average 18% of a dairy farm was dedicated to maize. The upscaling to the farm scale led us to include these areas in the calculation of standardized yield but the PP MCDM, even after adjustment for farm level, did not include these maize fields. Although this could have meant a small overestimation of standardized yield in farms with a successful crop maize it should affect farms with a medium and high PP equally. More importantly, all the grasslands with a medium score were in fact scored between a high and a medium PP by the MCDM. We believe that the difference in grass yield between farms with medium to high PP and farms with high PP is simply not large enough to be detected within our dataset. Either a larger sample size or a more diverse set of grasslands (with a broader PP spectrum) could shed further light on the relationship between standardized yield and the PP MCDM. Another approach would be to increase the number of final scores from three to five and increase the overall precision of the model without compromising interpretability.
Selecting proxy variables for NC was difficult due to the holistic nature of the NC MDCM model, which encompasses the process of NC and recovery from the application of residues/manure to the off-take of nutrients in the plant yield. The F:B ratio provided some assessment of the NC MDCM potentially representing the mineralization process within the model, which also contributes to the nutrient recovery (Schröder et al., 2016). We found a lower F:B biomass ratio in grasslands with high NC when compared to those with medium NC (although neither group was significantly different than grasslands with low NC capacity, which was probably because one of the four grasslands with low NC had an F:B ratio of 0.46, a value close to the mean F:B ratio of sites with high NC). Lower F:B ratios are typical in more intensively managed systems, with higher nitrogen inputs, and readily decomposable materials (Ruess & Ferris, 2004). We observed similarly low F:B values in all the cropland farms, regardless of their capacity to supply NC according to the MCDM, indicating a bacterial-dominated decomposition pathway for all the croplands.
The comparison of the B-HP function with the diversity of soil invertebrates resulted in promising trends that were, however, not significant, indicating a high level of overlap and variability. A significant but weak correlation T A B L E 1 Statistics from six redundancy analyses of the distribution of the high supply of primary productivity, nutrient cycling and biodiversity and habitat provision in grasslands and croplands, as explained by either environmental (Env), soil or management (Mng) variables. In bold are values with an adjusted R 2 avobe 0.3 or P-values below 0.05. Df stands for degrees of freedom. was observed when the B-HP MCDM scores were compared with those assigned by experts in the soil biodiversity field . Due to the discrete nature of the MCDMs, much of the variation in the data is discarded. The use of pesticides, irrigation or the use of drainage are all variables that can be either positive or negative as input data in the MCDMs, yet we know that pesticide load or the duration and timing of droughts (and therefore the resilience to water shortages) can both have an effect on the belowground biodiversity (Chelinho et al., 2011;Cyco n & Piotrowska-Seget, 2009;Korthals et al., 1996). We recommend either increasing the number of categories or using a continuous scale for these variables to further refine the B-HP MCDM, and although it can be costly, we recommend the inclusion of at least some of the soil fauna in the assessment of this soil function.

| Drivers of soil multifunctionality
The MCDMs indicated that up to 38% of the farms in our study had either a medium or a high supply of all three soil functions, which answers our first research question: a proportion of the farms under study were delivering three functions at a medium or high level simultaneously. Similar patterns have been observed recently in Europe, where more than a third of the visited fields achieved a supply of five functions at either medium or high levels (Zwetsloot et al., this issue). Only one farm in our study (a cropland) had high levels of all three functions. What is also significant is that only two farms (both grasslands) had a medium to high supply of only one soil function, which was PP in both cases, indicating that a vast majority of sites could supply at least two functions at a higher level simultaneously.
In answer to our second research question, the MCDM predicted clear synergies between PP and NC, but a trade-off between the two functions and B-HP was clearly evident (Figure 3). Only 27% of the farms in the study were expected to have a high level for B-HP. These results are in line with many previous studies that have found strong negative relationships between agricultural intensification and biodiversity losses (Geisen, Wall, & van der Putten, 2019;Tilman et al., 2001;Tsiafouli et al., 2015). Soil fauna are paramount for the delivery of soil functions and insuring the delivery of adequate B-HP should be a priority (Brussaard, 2012;Kibblewhite, Ritz, & Swift, 2008). Increases in organic fertilization, over time, increase the organic matter content in the soil and provide a more suitable habitat for soil biota by providing food for the microbial community and more stable soil moisture conditions (Birkhofer et al., 2008;Zsolnay & Görlitz, 1994). In our study, high soil organic matter was a common attribute amongst the croplands with a high B-HP. Additionally, lack of chemical pest management and high input of organic fertilizers were highlighted in our results as common input variables in farms with high B-HP. In fact, 80% of the farms with a high B-HP supply were organic farms. The use of insecticides and pesticides has been linked to a widespread loss of biodiversity in European farmland (Geiger et al., 2010). Conversion to organic farming can lead to more erratic yields, but measures such as the use of green manure or increased fertilization can be taken to reduce the variation (Knapp & van der Heijden, 2018). In fact, a meta-analysis showed that biodiversity increases on average by 30% under organic farming management, and this effect is larger in intensively managed landscapes (Tuck et al., 2014).
The type and amount of fertilizer were important variables in explaining multifunctionality: moderate values of each could lead to multifunctionality of grasslands and croplands. In grasslands, weather variables explained a significant proportion of the prevalence of high functioning. In croplands, however, moderate values of soil variables, such as soil organic matter, nitrogen to phosphorus ratio, carbon to nitrogen ratio and soil salinity, could lead to improvements in multifunctionality. Our results suggest that implementing changes in fertilization and pest management practices could have a big impact on the multifunctionality of soils even without a complete switch to organic management. Specifically, we found the type and amount of fertilizer to be an important variable explaining soil multifunctionality. Farms with a high capacity for B-HP were characterized by low mineral nitrogen fertilization and high organic fertilization, and in the case of croplands a soil organic matter content above 8%. The effects of mineral nitrogen fertilization on soil biota are often not direct. Their application can lead to shifts in the microbial composition by inducing changes in the PP, crop residue and soil organic matter content, but can also lead to soil acidification, which in turn has a negative effect on earthworms and nematodes (Bünemann, Schwenke, & Van Zwieten, 2006;Chen, Lan, Hu, & Bai, 2015;Eisenhauer, Cesarz, Koller, Worm, & Reich, 2012;Treseder, 2008); in fact, in our study, sites with the highest carbon to nitrogen ratio (which correlated strongly with low pH values) had the highest B-HP scores. There are measures to improve biodiversity that we have not tackled through the redundancy analysis ( Figure 5). The amount of residues left in the field, which we assumed to be low in our study sites, can positively impact the abundance and increase the functional diversity of the earthworm community (Frazão et al., 2019). Reduced tillage practices can lead to an increase in the richness and abundance of soil organisms (Kladivko, 2001;Sapkota, 2012). Moreover, these two practices combined can have a larger positive impact on the soil biodiversity than organic management (Henneron et al., 2015). The number and identity of crops in rotation can also have a positive impact on total organic carbon and microbial biomass, particularly when combined with the use of cover crops (McDaniel, Tiemann, & Grandy, 2014). Although our results have uncovered important relationships between PP, NC and B-HP, we could not include water or climate regulation functions, thus, some of the synergies and trade-offs that occur in these soils remain unknown. High nitrogen fertilization rates are associated with increases in productivity, but also to increased residual soil nitrate and increasing the risk of nitrous oxide (N 2 O) emissions (Snyder, Bruulsema, Jensen, & Fixen, 2009). In fact, Zwetsloot et al. (this issue) found a trade-off between PP and climate regulation in their study into European farms. The use of organic fertilizers has been linked with increased soil carbon stocks (Kukal, Rehana, & Benbi, 2009) and might therefore indicate a positive link between B-HP and climate mitigation. Although such a relationship was observed by Zwetsloot et al. (this issue) in Pannonian climatic soils, it became negative in the Atlantic climate zone, showing that there is still much to learn regarding how these soil functions interact with one another.
The functional land management framework was designed to promote soil multifunctionality by addressing soil function needs and delivery at different scales, such that local, regional and national objectives for multifunctionality can co-exist (Schulte et al., 2014;Schulte et al., 2015). Our results reflect a need to include the improvement of B-HP as a regional and national goal, to partially offset the effects of agricultural land management, but at the local scale sustainable agricultural practices should also be introduced to safeguard soil biodiversity. It is a step in the right direction that the Common Agricultural Policy in the European Union now includes "Effective Soil Management" as an objective, with specific recommendations regarding soil biodiversity, its role in soil functions, and management practices that are less impactful or even beneficial to the soil system. However, each member state can implement these measures at a different level and there is evidence to suggest that including environmental concerns in the Common Agricultural Policy is not enough to protect the environment (Alons, 2017;Pe'er et al., 2019). Better incentives for sustainable farming practices could lead to decreases in biodiversity loss, particularly because some of these practices are not without risk and farmers might not want to make changes that put their income at risk (Lefebvre et al.; Rotches-Ribalta et al.).

| CONCLUSIONS
Our study shows the potential for multifunctionality in soils in farms of the Netherlands, with most farms having a moderate or high capacity in two of three functions. The models predict trade-offs between biodiversity and habitat provision and both primary productivity and nutrient cycling. Management-related variables regarding the use of pesticides and the amounts of mineral and organic nitrogen fertilization were considered very important within the model predictions and explained part of the synergies and trade-offs in both croplands and grasslands, and soil parameters that can be changed via improved soil management (such as soil organic matter) are also important to increase the delivery of B-HP. Providing incentives that make the adoption of sustainable practices less risky to farmers might be a step forward in providing multiple soil functions simultaneously while maintaining the provision of B-HP.

ACKNOWLEDGMENTS
Special thanks to Jaap Schröder and Taru Sanden for their support in running the primary productivity and nutrient cycling multi-criteria decision models, and to Vladimir Kuzmanovski for his technical help in running all the models. We would also like to thank Auke Greijdanus and Ruud van der Meer for their help in navigating the Minerals Policy Monitoring Programme dataset. Lastly, we would like to thank Sylvana Harmsen for her aid during data analysis.

AUTHOR CONTRIBUTIONS
C.V., R.E.C. and R.G.M.G. conceived the presented idea and received important feedback from all co-authors. M. R. and T.K. provided the data necessary to perform the analyses. C.V. led the writing of the manuscript and the analysis of the data. All authors discussed the methods and results and contributed to the final manuscript.

DATA AVAILABILITY STATEMENT
The data that support this work were collected by the Netherlands National Institute for Public Health and the Environment and by Wageningen Economic Research and are protected by confidentiality agreements.