Multifunctionality of temperate alley-cropping agroforestry outperforms open cropland and grassland

Intensively managed open croplands are highly productive but often have deleterious environmental impacts. Temperate agroforestry potentially improves ecosystem functions, although comprehensive analysis is lacking. Here, we measured primary data on 47 indicators of seven ecosystem functions in croplands and 16 indicators of four ecosystem functions in grasslands to assess how alley-cropping agroforestry performs compared to open cropland and grassland. Carbon sequestration, habitat for soil biological activity, and wind erosion resistance improved for cropland agroforestry (P ≤ 0.03) whereas only carbon sequestration improved for grassland agroforestry (P < 0.01). In cropland agroforestry, soil nutrient cycling, soil greenhouse gas abatement, and water regulation did not improve, due to customary high fertilization rates. Alley-cropping agroforestry increased multifunctionality, compared to open croplands. To ameliorate the environmental benefits of agroforestry, more efficient use of nutrients is required. Financial incentives should focus on conversion of open croplands to alley-cropping agroforestry and incorporate fertilizer management. Alley-cropping agroforestry enhances carbon sequestration compared to open cropland and open grassland, and improves soil biological habitat and erosion resistance relative to open cropland, based on measured ecosystem functions at five sites in Germany.

C urrent agricultural practices in industrialized countries focus on high farm-level productivity and profit; however, external costs (e.g., soil degradation, water pollution, increased greenhouse gas (GHG) emissions, biodiversity loss) [1][2][3][4] are not included in the price and carried by the entire society. Nonetheless, intensively managed cropland monocultures undoubtedly show extraordinary achievements in agricultural production 5 . Several detrimental environmental consequences have raised awareness that modern agricultural systems should not only focus on high production but also on providing important ecosystem functions and landscape features, which stimulate biodiversity and carbon sequestration and reduce environmental pollution and soil degradation. Maintaining healthy soils and their functions are key requisites 6,7 in the pursuit of sustainable intensified agricultural systems 8 .
Agroforestry is projected as a promising form of agro-ecological management 9 . Presently, there is discussion to include financial incentives linked to agroforestry's environmental performance as compared to cropland monocultures, e.g., in the European Common Agricultural Policy 10 . Such financial incentives require comprehensive evaluation of ecosystem functions, including their capacity to provide several ecosystem functions simultaneously (also called "multifunctionality" 11 ). Our study fills this knowledge gap by comparing the multifunctionality of temperate alley-cropping agroforestry (i.e., the combined mechanized cultivation of rows of crops or grass alternated with rows of short-rotation trees [12][13][14] to open croplands and open grasslands without any trees. Whereas individual studies foster an increasing awareness of improved soil properties and ecosystem functions in temperate agroforestry, e.g., increases in soil organic carbon 12 , diversity of soil micro-and macroorganisms 13,15 , nutrient utilization efficiency 14 , wind erosion resistance 16 and reductions in nitrate leaching 17 , there is a lack of systematic comparisons of combined ecosystem functions between temperate agroforestry and open croplands or grasslands within a single multidisciplinary study that employs a replicated field-based design. This study aimed to quantify the multifunctionality of alleycropping agroforestry with short-rotation trees versus open croplands and open grasslands across different soil types and climatic conditions in Germany. Open croplands in our study were conventionally managed rotations of crop-monocultures (receiving customary applications of fertilizers and agrochemicals; Supplementary Table 1) without trees; open grasslands were permanent grassland without any trees. We used multiple distinct indicators of different ecosystem functions 9 Table 1). We hypothesized that alley-cropping agroforestry will promote beneficial ecosystem functions as compared to open croplands or open grasslands and will foster multifunctionality. Based on the ecosystem functions considered vital in assessing the benefits of agroforestry 9 , we quantified 47 indicators of seven ecosystem functions in croplands and 16 indicators of four ecosystem functions in grasslands (Supplementary Tables 2-4), which included the following: provision of food, fiber and fuel, carbon sequestration, soil nutrient cycling, habitat for soil biological activity, soil GHG abatement, water regulation, and erosion resistance.

Results
Indicators of provision of food, fiber, and fuel. Conversion of croplands and grasslands to alley-cropping agroforestry did not affect crop yield or grass biomass per cropped area or grassallocated area (Figs. 1 and 2). Although crop yield declined close to the tree rows, this was compensated by increased yield towards the center of the crop row ( Supplementary Fig. 2a). It should be noted that in our alley-cropping agroforestry systems, 80% of the area (consisting of 12 m tree row and 48 m crop or grass row) was planted to crop or grass; we also provided the crop yield or grass biomass per system's area (Supplementary Table 4) as these values may be useful for comparing with other temperate agroforestry designs. The gross margin analysis (see below) includes both crops and tree biomass that is harvested for biofuel feedstock (Supplementary Table 1). Fiber and protein contents remained unaffected by grassland agroforestry (Fig. 2), whereas crop quality partly improved for cropland agroforestry (Fig. 1) as shown by greater wheat (Triticum aestivum) crude starch and canola (Brassica napus) crude protein contents as well as greater canola 1000-grain weight compared to open cropland (P ≤ 0.03, Supplementary Table 4). In most cases, abundances of phytopathogenic fungi in crops were not affected following conversion to agroforestry, with the exception of Verticillium longisporum in canola, which was reduced in agroforestry compared to the open cropland (P = 0.03, Supplementary Table 4).
Indicators of carbon sequestration, soil nutrient cycling, habitat for soil biological activity, soil GHG abatement, water regulation, and erosion resistance. The large net primary production of trees strongly contributed to the increase in carbon sequestration (P < 0.01, Figs. 1 and 2) following agroforestry establishment in both cropland and grassland agroforestry. In cropland agroforestry, fine root density increased compared to open cropland (P < 0.01, Supplementary Table 4); however, soil organic C (SOC) stocks did not differ between agroforestry and open cropland or grassland (P > 0.05, Supplementary Table 4). Conversion of open cropland to agroforestry improved the soil biological habitat (P = 0.03, Fig. 1) as shown by greater earthworm biomass (P < 0.01, Supplementary Fig. 2b) and soil bacterial and fungal population sizes (P ≤ 0.03, Supplementary Table 4). In cropland agroforestry, the tree rows greatly reduced wind speed and number of days with wind erosion risk compared to open cropland ( Supplementary Fig. 2c), which resulted in a substantial increase in wind erosion resistance (P < 0.01, Fig. 1). Although soil nutrient cycling, soil GHG abatement and water regulation functions did not change following conversion of open cropland or grassland to agroforestry (P > 0.05, Supplementary Table 5), two indicators improved: larger plant-available P in grassland agroforestry than in open grassland (P = 0.02) and greater gross rates of nitrous oxide (N 2 O) uptake in the soil under cropland agroforestry compared to open cropland (P = 0.01, Supplementary Table 4). In summary, conversion to agroforestry enhanced carbon sequestration, habitat for soil biological activity and erosion resistance functions in croplands ( Fig. 1) and improved carbon sequestration in grasslands (Fig. 2). Implementation of agroforestry did not lead to a decrease in any of the measured indicators of ecosystem functions (Figs. 1 and 2 and  Supplementary Tables 4 and 5).
Correlations between ecosystem functions. This analysis was performed to identify whether there were potential trade-off or win-win situation between ecosystem functions-suitable to carry out when multiple ecosystem functions are quantified. We did not detect any correlations between ecosystem functions in cropland agroforestry; however, carbon sequestration and nutrient cycling were positively correlated in grassland agroforestry ( Supplementary Figs. 3a and 4a). There was no correlation between ecosystem functions in open grassland although carbon sequestration was positively correlated with provision of food, fiber and fuel in open cropland. In contrast, water regulation was

Discussion
Cropland agroforestry substantially improved ecosystem functions without affecting yield of the cropped area. Long-term yield stability was supported by an eight-year investigation at one of our sites, in which crop yields are similar between agroforestry crop row and open cropland 18 . This observed yield stability was probably related to the combined effect of limited tree height (a result of short rotation management) and relatively wide (48 m) crop rows, which may limit competition for light and other resources. In two cropland sites on Vertic and Gleyic Cambisol soils (Supplementary Table 1), where there was a decrease in crop yield at 1-m distance from the tree row, these were compensated by the increase in yield in the central part of the crop rows 14 ( Supplementary Fig. 2a). However, it should be noted that for other temperate alley-cropping agroforestry systems, where trees are harvested less frequently (e.g., timber trees) and/or the crop rows are not wide, reduced crop yields with increasing tree heights have been observed 19,20 .
The three improved ecosystem functions under cropland agroforestry, i.e., carbon sequestration, habitat for soil biological activity and wind erosion resistance (Figs. 1-3), all address major global change problems. First, the increases in the trees' net primary production and combined root density from crops and trees signified that agroforestry performed better in carbon sequestration than open cropland or grassland (Fig. 3). This substantiates projections that European agricultural systems currently sequester less carbon than they potentially could under improved land-use strategies 21 including agroforestry 22 . Although SOC stocks in our ≤9-year-old agroforestry systems did not differ from the open cropland or grassland, other studies in ≥15-year-old alley-cropping agroforestry systems show increases in SOC stocks 23 . As claimed for temperate agroforestry systems based on a meta-analysis 12 , we expect that SOC stocks at our sites will also increase as agroforestry ages. Second, the increased population sizes of earthworms, bacteria, and fungi in our cropland agroforestry systems (Fig. 1) support positive effects of agroforestry for biological habitat, as was reported earlier for macro-detritivorous arthropods 13 , earthworms 24 , and plant diversity 25 . This improved habitat for soil biological activity will possibly have positive bottom-up effects on higher trophic levels 26 . Hence, agroforestry may counter the effect of intensive management in croplands on loss in species abundance 27 . In grassland agroforestry, however, habitat for soil biological activity remained comparable to that in open grassland (P > 0.05, Supplementary Table 5), which we attribute to the fact that both systems have perennial plants, manifesting permanent roots 28  Table 5). and high plant species diversity 29 combined with low fertilizer application at the studied sites 30 . Third, susceptibility of soil to wind erosion is not only common in Central Europe 31 but is a serious global problem 32 that reduces SOC, nutrient stocks, and agricultural productivity 33,34 . The strong reduction of wind speeds to levels below risk of wind erosion ( Supplementary  Fig. 2c) was a direct result of the introduction of tree rows (Fig. 1), a positive effect that is well known from shelterbelts 35 . This is a major motivation mentioned by farmers for their acceptance of agroforestry 36 .
Alley-cropping agroforestry did not perform better than open cropland or grassland in ecosystem functions related to fertilization, i.e., soil nutrient cycling, soil GHG abatement and water regulation ( Fig. 1), although such functions were reported to improve in other agroforestry systems 37,38 . Soil N 2 O and CH 4 fluxes are inherently variable in space and time, which can affect the results if the spatial and temporal sampling schemes do not take these into account. We dealt with this challenge by including in our spatial design potential effects of distance to the tree rows ( Supplementary Fig. 1). Our temporal resolution was monthly sampling from March 2018 to February 2019, which included periods following fertilizer application when soil mineral N, water content and GHG fluxes are strongly elevated 39 . Thus, we are confident that our sampling design was able to capture a substantial part of the spatial and temporal patterns of soil GHG fluxes; the annual soil GHG fluxes used in our analysis (Supplementary Table 4) were interpolations across months (including the winter period) and encompassed any spatial patterns, as the values for agroforestry were weighted by the areal coverages of the tree row and the various sampling distances within the crop row 39 . Hence, the lack of difference in soil GHG abatement function between the two management systems was probably real and not caused by the inadequate sampling design. Management practices, including fertilization rates, at our sites followed customary management operations in Germany (Supplementary Table 1) and were the same for agroforestry crop row and open cropland. Both agroforestry crop rows and open croplands show low nutrient response efficiency (yield to plantavailable nutrient ratio) and both systems are nutrient saturated 14 . Thus, potential improvements of nutrient cycling following the introduction of tree rows, such as nutrient recycling through litterfall 40 and permanent tree roots acting as a safety net for excess nutrients beyond uptake of crop roots 41 , did not recuperate the ecosystem functions related to fertilization and were probably masked by the current high fertilization rates of the agroforestry crop rows. This suggests opportunities to reduce fertilizer application rates. In our cropland agroforestry systems, a reduction of 50% in available nutrients could increase the nutrient response efficiency by 67% for nitrogen (N) and 83% for phosphorus (P) with only relatively small (17% for N; 8% for P) reductions in yield 14 . We anticipate that in the agroforestry systems an improved nutrient cycling through tree leaf-litter input 14 and nutrient uptake by deep tree roots with associated mycorrhizal fungi 42 , a safety net against leaching losses, may even moderate yield losses. However, reductions in fertilization rates still lack testing under field conditions in both systems. Additional benefits of reduced fertilization would be reduced leaching of excess nutrients and soil N 2 O emissions 39,43 and the air pollutant nitric oxide 44 . The increased gross N 2 O uptake in the soil of cropland agroforestry was brought about by the unfertilized tree rows, which have low ratio of soil mineral N to decomposable carbon (e.g., due to absence of fertilization in the tree row combined with litter input from trees) and favors for microbial reduction of N 2 O to N 2 45 . This attests to the benefit of agroforestry in increasing uptake of N 2 O in soils 45 .
Although conversion of open cropland to alley-cropping agroforestry with short-rotation trees has been shown to increase the gross margin (including crops and tree biomass used for biofuel feedstock) from €489 ± 5 ha −1 to €518 ± 5 ha −1 46 , a recent farmers' survey revealed that a major obstacle for agroforestry implementation is establishment cost, which was estimated at €1800 ha −1 36 . Additional perceived hindrances were irregular income from harvest of aboveground tree biomass (commonly carried out every four to seven years for biofuel feedstock) and price fluctuations of wood 46 . Yet, there appeared to be a widespread willingness among the interviewed farmers to establish agroforestry if these financial obstacles would be ameliorated. The financial compensation that the farmers perceived as necessary to encourage agroforestry establishment was much higher than the above-mentioned costs. Farmers estimated an initial financial support for establishing agroforestry on their fields at €2718 ± 345 ha −1 , followed by an annual support of €511 ± 54 ha −1 36 . Interviewed taxpayers, outside the farming community, showed a willingness to pay a tax with a median value of €20 year −1 per person for environmental benefits and landscape amenities associated with agroforestry 47 .
The full ecological and economic potential of agroforestry will only be achieved if nutrient management practices are optimized. The following suggestions can help to reach this goal. First, current fertilization rates need to be adapted to the local yield levels (not on generalized maximum yield potential) as overfertilization is the major cause of inefficient nutrient use, triggering high external environmental costs 1,7 . Commonly, crops use less than half of applied N fertilizer 48 . As discussed earlier, there is potential to reduce applied nutrient rates which will improve the functions of soil nutrient cycling, soil GHG abatement and water regulation with only marginal yield losses 14 . Second, the introduction of tree rows in agroforestry causes spatial and temporal variations of not only nutrient inputs from leaf litter and root turnover but also of light and evapotranspiration 49 , which cause spatial variability in yield 18 . Using this knowledge, precision farming techniques have the potential to further optimize nutrient inputs both in space and time. At present, fertilizer is applied using standard farm equipment in our studied cropland agroforestry and open cropland systems, which leads to uniformly high nutrient input even in areas closest to the tree rows with low yields 14,18,50 . Third, agroforestry may be further optimized by diversifying crops: (a) inclusion of shade-tolerant crops in the rotations as tree height increases with age, and (b) growing shade-tolerant crops near the tree rows in combination with light-demanding crops in the common rotations. We anticipate that these management improvements will reduce the current trade-offs between soil nutrient cycling, habitat soil biological activity and water regulation in croplands, which possibly facilitate more positive functional interactions than are currently observed (Supplementary Figs. 3 and 4). With the exception of carbon sequestration, ecosystem functions in grasslands did not change following conversion to agroforestry. Furthermore, permanent open grasslands in Europe are protected from conversion into other land uses 51 , which calls for prioritizing financial incentives for the conversion of open cropland to agroforestry.
Our findings support our main hypothesis that alley-cropping agroforestry with short-rotation trees enhances multifunctionality and hence represents a more sustainable management compared to open cropland or grassland, even under current high fertilization rates. Nevertheless, the full potential of cropland agroforestry as a sustainable and multifunctional management alternative can only unfold if policy changes are implemented that provide incentives for farmers to adjust their management practices. Therefore, appropriate financial support should not only address the financial risks in establishing agroforestry systems 52 but also focus on measures that further reduce external environmental costs 9,53 , which are mainly associated with current high fertilizer applications. Although there are several promising strategies that can be implemented in agroforestry to improve nutrient management 14 , these concepts still need long-term evaluation under actual field conditions. The current initiative of the European Commission to couple environmentally beneficial agricultural practices with a substantial fraction of the general subsidy for agricultural land 54 provides an opportunity to create such incentives to promote and optimize agroforestry in our pursuit towards a more sustainable yet productive agriculture.

Methods
Study area, experimental design and management. Our study was conducted on three cropland systems and two grassland systems (Supplementary Fig. 1). The croplands were located in Thuringia (Calcaric Phaeozem 55 soil), Lower Saxony (Vertic Cambisol 55 soil) and Brandenburg (Gleyic Cambisol 55 soil), and the grassland sites were located in Lower Saxony (Histosol and Anthrosol 55 soils), Germany. At each site, we studied alley-cropping agroforestry (12 m wide rows of trees alternated by 48 m wide rows of crops or grassland) paired with open cropland or open grassland. Each of these management systems at each site was represented by four replicate plots in cropland systems and three replicate plots in grassland systems (total of 3 sites × 2 management systems × 4 replicate plots = 24 plots in croplands, and 2 sites × 2 management systems × 3 replicate plots = 12 plots in grasslands). To take into account any possible spatial variations of the measured indicators of ecosystem functions (see details below) with distance from the tree row, measurements at each agroforestry replicate plot were conducted in the tree row and at various distances from the tree row within the crop or grass row (Supplementary Figs. 1 and 2). This sampling scheme was employed for all indicators, unless described otherwise (see ecosystem function indicators below). Our sites are managed following the farmers' customary crop rotation (Supplementary Table 1). Agroforestry crop rows and open croplands are conventionally managed, which included cultivation once a year and application of recommended mineral fertilizers and agrochemicals. Tree rows (see Supplementary Table 1 for tree species, distances between rows and between trees within a row, planting density, and first cutting of wood biomass) were not fertilized and were harvested after four to seven years for bioenergy, i.e., the wood biomass was exported from the field. Details on site and management practices can be found in Supplementary Table 1.
Multifunctionality and ecosystem functions. Several approaches have been applied to assess multifunctionality of land use or management systems. So far, there is no consensus which approach is preferable; however, a recent review made several recommendations 56 , which we largely followed in our study. We wanted to test whether the arguments that are typically used in favor of agroforestry 9 can be quantified in an ecosystem function framework for our studied management systems, and whether this resulted in differences in multifunctionality. We started with a list of ecosystem functions (categorized under provisioning, regulating and supporting services) that are considered vital in assessing benefits of temperate agroforestry systems 9 : provision of food, fiber and fuel, air quality regulation, climate regulation (often separated into carbon sequestration and GHG abatement), water quality and regulation, erosion regulation, pest and disease regulation, biodiversity regulation, and regulation of natural hazards and extreme events. Finally, soil formation and nutrient cycling are often included as an important ecosystem function that support other functions 9 . The ecosystem functions and their indicators (Supplementary Tables 2-4), which we measured in order to compare alley-cropping agroforestry with open cropland or grassland, were adjusted for various reasons. Regulation of air quality and of natural hazards and extreme events were not available, because these effects can only be measured for larger areas and not at the scale of the experimental plots in our study. We also excluded soil formation (categorized under supporting service) because this function is often indicated by SOC and soil microbial communities, which we measured but used as indicators for the more fitting ecosystem functions (see below). We separated climate regulation into carbon sequestration and soil GHG abatement because stakeholders often express interest in carbon sequestration (e.g., to calculate carbon credits) 9 . Finally, we measured phytopathogen abundances and these were included as indicators under provision of food, fiber and fuel (and not for pest and disease regulation), because phytopathogen abundances directly influence yield quantity and quality. Below are the details of the ecosystem functions and their corresponding indicators that were measured to evaluate the contrasting management systems. For each function, we explain why we used those indicators and how they were measured. Although some indicators can be used in multiple ecosystem functions, we used each indicator in only one ecosystem function so as not to out-weigh its effect on aggregated ecosystem functions to depict multifunctionality (Fig. 3).
Provision of food, fiber, and fuel. This function was assessed in all the cropland and grassland sites. We used indicators not only of grain yield but also of yield quality (affecting the marketability and suitability for food processing) and of food safety (abundances of different phytopathogens). Grain yield was determined using a combine harvester to sample an area of 17.5 m 2 at specific sampling distances from the trees (Supplementary Fig. 1 Carbon sequestration. This function was assessed in all the cropland and grassland sites. The indicators included tree aboveground net primary production, as this short rotation coppice is used for carbon-neutral energy source; root density, which through its turnover may be a timely indicator of SOC source; and SOC stock that serve as a surrogate for its slow-changing process 50 . Tree aboveground net primary production was calculated as the sum of woody biomass production and leaf litterfall during 2016 and 2017. Woody biomass production was measured for 50% of all trees in a plot by measuring stem diameter at breast height (DBH) and applying a site-and year-specific allometric equation 59 , derived from the stem diameters and measured biomass of 25 trees that represented the DBH range of the trees present. Tree leaf litterfall was collected bi-weekly in 2016 and 2017 using 0.14 m 2 litter traps. Combined fine root density of trees and crops was determined within 0-100 cm depth with soil cores of 6.5 cm diameter in the growing season of 2017, by washing out the soils and collecting all roots of ≥1 cm length and drying them at 70°C for three days. SOC was measured for the top 30 cm depth in 2016 from air-dried, 2-mm-sieved and acid-fumigated 60 (in case of pH ≥6) soils, using a CN analyzer (Elementar Vario EL; Elementar Analysis Systems GmbH, Hanau, Germany). Soil bulk density was measured in 2016 during the growing season for the top 30 cm using stainless steel soil cores with a volume of 250 cm 3 61 . As there was no statistical difference in soil bulk density among agroforestry tree row, crop row and open cropland or grassland at each site, we used the average soil bulk density for each site in order to have the same soil mass for the calculations of SOC stocks in the two management systems 14 .
Soil nutrient cycling. This function was assessed in all the cropland and grassland sites. The indicators included gross rate of N mineralization in the soil (a process providing extant mineral N for plant and microbial use), plant-available P and K (which reflect the P and K stocks in the soil accessible for root uptake), base saturation (signifying the proportion of base nutrient elements adsorbed on the soil exchange sites), effective cation exchange capacity (indicating the capacity of the soil to adsorb elements and thus an index of its slow changing process), and nifH gene abundance (denoting the genetic N 2 fixation potential, responsible for encoding nitrogenase enzyme in soils). Soil gross N mineralization was measured in the top 5 cm between April and June 2017 using 15 N pool dilution techniques 62 on intact in situ incubated soil cores. Plant-available phosphorus (P) was measured monthly during the growing seasons of 2016 and 2017 in the top 5 cm depth, and the average value over the growing season represented each replicate plot. Plantavailable P was determined as the sum of bicarbonate-extractable and resinexchangeable P 63 , analyzed using an inductively coupled plasma-atomic emission spectrometer (ICP-AES; iCAP 6300 Duo View ICP Spectrometer, Thermo Fischer Scientific GmbH, Dreieich, Germany). Plant-available potassium (K) was determined for the top 30 cm depth in 2016 by percolating the soils with unbuffered 1 mol l −1 NH 4 Cl and the percolates were analyzed for exchangeable K (as well as for cations mentioned below) using ICP-AES. Base saturation was measured as the percentage of base cations on the effective cation exchange capacity, which was determined as the sum of exchangeable Ca, Mg, K, Na, Mn, Al, Fe, and H. Soil N 2 -fixation nifH gene abundance in the top 5 cm was determined in April 2017 for grasslands and March 2019 for croplands using real-time PCR. Sampled soils were frozen immediately in the field for DNA extraction 30,64 and nifH genes were amplified using the primer pair IGK3-DVV 65 .
Habitat for soil biological activity. This function was assessed in all the cropland and grassland sites. The indicators were earthworm biomass (denoting bioturbation which influences soil aggregation), microbial biomass C and N (as expressions of microbial biomass size which influences nutrient cycling), bacterial and fungal abundances (as expressions of prokaryotic and eukaryotic microbial population size), and β-glucosidase activity (implying the potential for enzymatic break down of complex carbohydrates and thus reflects the acquisition for energy by the microbial biomass). Earthworm biomass was determined in 2018 by extracting individuals from 0.25 × 0.25 × 0.25 m soil blocks using hand-sorting 66 . Soil microbial biomass C and N were determined in the top 5 cm in spring 2017 by CHCl 3 fumigation-extraction method 67 . One pair of soil samples was immediately extracted with 0.5 mol/l K 2 SO 4 , and the other pair was fumigated for 5 days and then extracted. Extractable organic C was analyzed by UV-enhanced persulfate oxidation using a carbon analyzer with an infrared detector (TOC-VWP, Shimadzu Europa GmbH, Duisburg, Germany). Extractable N was determined by ultravioletpersulfate digestion followed by hydrazine sulfate reduction using continuous flow injection colorimetry (Autoanalyzer Method G-157-96; SEAL Analytical AA3, SEAL Analytical GmbH, Norderstedt, Germany). Microbial biomass C and N were calculated respectively as the differences in extractable organic C and N between fumigated and unfumigated soils divided by k C = 0.45 and k N = 0.68 for 5-day fumigated samples 67,68 . Bacterial and fungal abundance in soils were measured in spring 2017 at the grasslands and in spring 2019 at the croplands. Soil samples collected in the top 5 cm were frozen immediately in the field for DNA extraction, and the abundance of soil bacteria and fungi was determined using real-time PCR 30,64 . The soil β-glucosidase activity in the top 5 cm, sampled in October-November 2015, were analyzed using a fluorescently labeled substrate (4methylumbelliferone-β-D-glucopyranoside; Sigma-Aldrich Chemie GmbH, Taufkirchen, Germany) 69 . Fluorescence intensity was measured at 30-minute intervals over 180 minutes at 355 nm excitation and 460 nm emission wavelengths using a FLUOstar Omega microplate reader (BMG, Offenburg, Germany).
Soil greenhouse gas abatement. This function was assessed in the cropland sites, and the indicators were net fluxes of the GHG, N 2 O and CH 4 , from the soil surface, and the gross rate of N 2 O uptake in the soil (signifying the capacity of the soil to convert the GHG N 2 O into unreactive N 2 via microbial denitrification) 45 . Soil CH 4 and N 2 O fluxes were measured monthly from March 2018 to February 2019 using vented static chambers 70 . Gas samples were analyzed using a gas chromatograph (SRI 8610C, SRI Instruments Europe GmbH, Bad Honnef, Germany) with a flame ionization detector (for CH 4 concentration) and an electron capture detector (for N 2 O concentration). Gross rate of N 2 O uptake in the soil was measured monthly from March 2018 to February 2019 using the 15 N 2 O pool dilution technique 71 . We injected 7 ml of 15 N 2 O label gas into the chamber headspace, which was composed of 100 ppm v of 98% single labeled 15 N-N 2 O, 275 ppbv sulfur hexafluoride (SF 6 , as a tracer for possible physical loss of gases from the chamber headspace) and balanced with synthetic air (Westfalen AG, Münster, Germany). At 0.5, 1, 2, and 3 h in situ incubation, gas samples of 100 ml and 23 ml were taken from the chamber headspace and injected respectively into a pre-evacuated 100 ml glass bottle and a 12 ml glass vial (Exetainer; Labco Limited, Lampeter, UK) equipped with rubber septa. The 100 ml gas samples were analyzed for isotopic composition using an isotope ratio mass spectrometer (Finnigan Deltaplus XP, Thermo Electron Corporation, Bremen, Germany). The 23 ml gas samples were analyzed for N 2 O and SF 6 concentrations, using the same gas chromatograph with electron capture detector. Calculation of gross N 2 O uptake in the soil is given in our earlier study 71 .
Water regulation. This function was assessed in the cropland sites, and the indicators were: actual evapotranspiration (the amount of water taken up by the plants), soil saturated conductivity (Ks) (indicating the flow of water in the soil), change in water storage in the soil (signifying the amount of water available for plant uptake), and leaching fluxes of N, P, S, and base cations (implying decrease in quality of drainage water in the soil). Actual evapotranspiration was determined in 2016 and 2017 with eddy covariance masts (3.5 m height in open cropland, 10 m height in agroforestry), using the eddy covariance energy balance (ECEB) method 49 . The ECEB measurements integrate over the entire footprint area, and the entire cropland agroforestry and open cropland at each of the three cropland sites were represented each by an annual value of actual evapotranspiration. The annual value was the sum of the half-hourly, ECEB-derived actual evapotranspiration rates. The Ks was measured in July-August 2019 by taking intact soil cores (250 cm 3 each) that were slowly saturated with water in the laboratory (after securing a filter paper to one end of the cores and placing them on a basin filled with 1 cm of degassed water). The water level was raised over the course of 8 h to avoid trapping air in the soil samples. Once saturated, the intact soil cores was installed on the UMS-KSAT device (UMS, München, Germany) using either falling head or constant head method depending on the permeability of the soil 72 . Measurements of Ks were based on the Darcy equation and were normalized to 20°C. Change in water storage was estimated using the water sub-model of the Expert-N 5.0 model 73 for the top 60 cm of the soil for the period of April 2016 to March 2017. The same water model was used to estimate monthly drainage fluxes for calculation of nutrient leaching fluxes 74 . Soil water was sampled monthly from 2016 to 2018 at 60 cm depth, using suction-cup lysimeters (P80 ceramic, maximum pore size 1 μm; CeramTec AG, Marktredwitz, Germany), and was analyzed for dissolved total N (using the same continuous flow injection colorimetry mentioned above), total P, total S, Na, K, Ca, and Mg concentrations (using the same ICP-AES mentioned above). Soil texture was measured in 2016 during the growing season for the top 30 cm, using the pipette method with pre-treatments to remove organic matter and carbonate (for soil pH >6) 75 . The water sub-model inputs of Expert-N were climate (global radiation, temperature, precipitation, relative humidity and wind speed), soil (texture, bulk density, and Ks), all measured at our study sites, and vegetation characteristics (root biomass and leaf area index), specific to the crops or trees at the sites. Modeled soil water contents were validated with measured soil moisture contents, conducted monthly by gravimetric measurements. Nutrient leaching fluxes were calculated by multiplying the nutrient concentrations with the water drainage fluxes during the sampling period and summed for the entire year 74,76 . open cropland act as a reference wind field, from which the wind speed measurements at each of the replicate plots of cropland agroforestry (influenced by the tree rows) were compared to. We used a LES model to simulate the spatially varying wind field for three wind directions (north, north-west and west) in the agroforestry tree row (with tree height set to 5 m), agroforestry crop row and open cropland, which would be impossible to achieve with only wind speed measurements. We extracted the wind speed at the lowest model domain height of 0.5 m above the ground from the 3D model output for the three wind directions in the agroforestry and open cropland 77 . We then derived a ratio of horizontal wind speed of these three wind directions in the agroforestry to horizontal wind speed of the same directions in the open cropland, as indicators of wind speed reduction. This wind speed ratio from the model output was then multiplied with half-hourly mean wind speed measurements from the open cropland to derive a temporally and spatially varying wind speed. The spatially varying potential for wind erosion was obtained by the number of days with wind speed exceeding a threshold wind speed of 5 m s −1 at 0.5 m above the ground for each management system and was used as indicator for wind erosion resistance 78,79 . We only calculate potential wind erosion for this one cropland site because (i) its paired agroforestry offered a large extended area with different planting densities of trees, and (ii) this cropland site had independent measurements of wind speed in a horizontal transect 80 , both allowed validation of our simulated values at wide distances between tree rows and under different planting densities. Furthermore, from this cropland site we were able to understand effects of distance between tree rows and tree density on wind speed reduction. As the agroforestry plots at our other two cropland sites had comparable tree and crop row widths, we expect that the findings from this site are also relevant to those other sites in our study.
Data compilation and statistical analysis. To derive the value for a whole replicate plot of agroforestry, the values from the tree row and various sampling distances within the crop row were weighted by their areal coverages. As the agroforestry systems were represented by a 12 m tree row and a 48 m crop row, the weighting factors were: 12/60 for the tree row, 5/60 for the 1 m, 6/60 for the 4 m, 33/60 for the 7 m that represented a wide area of the crop row based on comparison tests among these sampling points using suitable indicators, e.g., for provision of food/fiber/fuel, habitat for microbial activity and GHG abatement 14,30,45,64 ; and 4/60 for the 24 m sampling distances within the crop row that accounted the 4-m overlap of fertilizer broadcaster as it turned around in the middle of the crop row 45 . When indicators were not measured at all sampling distances within the crop row (e.g., 1 and 4 m distances were statistically comparable in some indicators, and further measurements excluded the 4 m but instead increased the temporal sampling), we adjusted the weighting factors for plot-level agroforestry values based on the areal coverages of the sampled distances within the crop and tree rows. For each indicator, area-weighted values for agroforestry plots and open cropland or grassland plots were z-standardized (z = (individual value − mean value across plots and sites)/standard deviation). This prevents the dominance of one or few indicators over the others (e.g., brought by the units of their values) within a specific ecosystem function, and z-standardization allows several distinct indicators to best characterize an ecosystem function. The z-standardized, plot-level values were used in assessing differences between management systems (i.e., agroforestry versus open cropland or open grassland) for a specific ecosystem function with multiple distinct indicators 81 . For indicators in which high values signify an increased undesirable effect (i.e., grain pathogen abundance, nutrient leaching losses, soil GHG fluxes, wind speed as erosion risk indicator), we used the additive inverse (multiplied by −1) of their values in order to have uniform intuitive meaning as the rest of the indicators (i.e., large value indicates high desirable effect). Differences in specific ecosystem functions between agroforestry and open cropland or open grassland were assessed using linear mixed effects (LME) model with management system (cropland or grassland agroforestry vs. open cropland or grassland) as fixed effect and indicators and replicate plots across sites (4 plots × 2 management systems × 3 sites for cropland systems = 24 plots, and 3 plots × 2 management systems × 2 sites for grassland systems = 12 plots) as random effects (Supplementary Fig. 1 and Supplementary Tables 2-5) 81 . As indicator variables may systematically differ in their responses to management systems, we also tested the interaction between indicator and management system (Supplementary Table 5) 81 . Diagnostic plots of LME model residuals were visually inspected to check whether these meet the criteria for normal distribution and equal variance. To assess whether there were potential trade-off or win-win situations, we conducted Spearman rank correlation test among ecosystem functions within each management system on the z-standardized values of replicate plots (Supplementary Figs. 3a, b and 4a, b); such analysis is commonly perform when multiple ecosystem functions are quantified 81 . Additionally, we conducted Spearman rank correlation test among ecosystem functions attributable to conversion from open cropland or grassland to agroforestry using the relative change of z-standardized values (i.e., agroforestry − open cropland or grassland/open cropland or grassland) ( Supplementary Figs. 3c and 4c). For testing differences in specific indicators between management systems, we used either an Independent t test (when the data showed equal variance and normal distribution), a Wilcoxon test (when the data showed equal variance but non-normal distribution) or a Welsh test (when normal distribution and equality of variance were not met). Data were analyzed using R (version 4.0.4).
Reporting summary. Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.