Mechanical weeding enhances ecosystem multifunctionality and profit in industrial oil palm

Oil palm is the most productive oil crop, but its high productivity is associated with conventional management (that is, high fertilization rates and herbicide application), causing deleterious environmental impacts. Using a 22 factorial experiment, we assessed the effects of conventional vs reduced (equal to nutrients removed by fruit harvest) fertilization rates and herbicide vs mechanical weeding on ecosystem functions, biodiversity and profitability. Analysing across multiple ecosystem functions, mechanical weeding exhibited higher multifunctionality than herbicide treatment, although this effect was concealed when evaluating only for individual functions. Biodiversity was also enhanced, driven by 33% more plant species under mechanical weeding. Compared with conventional management, reduced fertilization and mechanical weeding increased profit by 12% and relative gross margin by 11% due to reductions in material costs, while attaining similar yields. Mechanical weeding with reduced, compensatory fertilization in mature oil palm plantations is a tenable management option for enhancing ecosystem multifunctionality and biodiversity and increasing profit, providing win–win situations. A full factorial experiment in a state-owned industrial oil palm plantation in Indonesia evaluates whether reduced management via reduced fertilization rates and mechanical weeding can decrease the negative impacts on ecosystem functions and biodiversity while maintaining profitability of the plantation.

Oil palm is the most productive oil crop, but its high productivity is associated with conventional management (that is, high fertilization rates and herbicide application), causing deleterious environmental impacts. Using a 2 2 factorial experiment, we assessed the effects of conventional vs reduced (equal to nutrients removed by fruit harvest) fertilization rates and herbicide vs mechanical weeding on ecosystem functions, biodiversity and profitability. Analysing across multiple ecosystem functions, mechanical weeding exhibited higher multifunctionality than herbicide treatment, although this effect was concealed when evaluating only for individual functions. Biodiversity was also enhanced, driven by 33% more plant species under mechanical weeding. Compared with conventional management, reduced fertilization and mechanical weeding increased profit by 12% and relative gross margin by 11% due to reductions in material costs, while attaining similar yields. Mechanical weeding with reduced, compensatory fertilization in mature oil palm plantations is a tenable management option for enhancing e co sy stem m ul tifunctionality and biodiversity and increasing profit, providing win-win situations.
Oil palm production has increasingly expanded in large areas of Southeast Asia, with Indonesia currently the world's largest producer of palm oil 1 , which also coincides with increased rates of deforestation in the country 2 . It is estimated that between 2001 and 2019, oil palm plantations were responsible for 32% of the total forest area lost in Indonesia 3 . Despite the high environmental costs, oil palm production is highly attractive because of its short-term economic returns and increasing global demand for food, fuel and cosmetics 4 .
Industrial oil palm plantations (>50 ha planted area and owned by corporations 5 ) are currently estimated to account for about 60% of the total cultivated oil palm area in Indonesia 6 . Compared with smallholder farms (typically <50 ha of land and owned by individuals 5 ), the productivity of industrial oil palm plantations is ~50% higher 7 but is largely driven by high fertilizer and herbicide applications to optimize productivity 8 at the expense of other ecosystem processes. In contrast to the forests, industrial oil palm plantations tend to be structurally simplified and highly disturbed [9][10][11] and have reduced capacity to provide several ecosystem functions simultaneously (so-called 'multifunctionality' of ecosystems) 4,12 . Although various studies have addressed the deleterious effects of oil palm expansion on forest and biodiversity losses 4,12,13 , to date, there has been no holistic multiyear and spatially replicated assessment of the effect of different oil palm Article https://doi.org/10.1038/s41893-023-01076-x rates × two weeding practices, with conventional vs reduced (that is, equal to nutrients exported with harvest) fertilization rates and herbicide vs mechanical weeding (that is, brush cutter), was established in 2016 in a mature (≥16 years old) industrial oil palm plantation. Treatments had four replicate plots of 50 m × 50 m each. During the first four years of the experiment, we measured indicators of eight ecosystem functions, seven indicators of biodiversity and six economic indicators linked to the level and stability of yield and profit (Supplementary  Tables 1 and 2). We hypothesized that compared with conventional high fertilization rate and herbicide treatment, reduced fertilization rate and mechanical weeding will enhance ecosystem functions and biodiversity while maintaining high productivity and increased profit.

Management effects on ecosystem functions and biodiversity
Ecosystem multifunctionality, an aggregated measure of multiple ecosystem functions, was, on average, higher in mechanical weeding than herbicide treatment ( Fig. 1 and Table 1). Similarly, threshold multifunctionality, calculated as the number of ecosystem functions above a specific threshold (Methods and Supplementary Fig. 3), showed a higher number of functions in mechanical weeding than herbicide treatment that were above 70% and 90% thresholds and marginally significant at thresholds of 50% and 80% (Table 1 and Supplementary  Fig. 3). Although neither fertilization nor weeding effect was detected when analysing only for each individual ecosystem function (Table 1 and Supplementary Fig. 1), the higher multifunctionality of mechanical weeding over herbicide treatment was brought about by the higher mean z values of measured indicators of litter decomposition, soil fertility, water filtration and plant refugium ( Supplementary Fig. 2). management strategies on ecosystem functions, biodiversity and economic productivity. Primary, experimental research studies that assess the effect of management on multiple outcomes, although presently rare, are pertinent to providing clearer oil palm management recommendations 14 because oil palm plantations are inherently complex agricultural systems. Application of fertilizer and herbicide at high frequencies are arguably the most important management activities in industrial oil palm plantations due to their effect on oil palm yield. Thus, a major leverage for a more environmentally sustainable production of palm oil is identifying optimal fertilizer application rates and weeding method to maintain sufficiently high productivity and economic profitability while minimizing associated losses of ecosystem functions and biodiversity.
Fertilization contributes to maintaining high productivity levels but also represents a substantial share of management costs in industrial plantations. High fertilization rates are associated with high nutrient-leaching losses 15 and greenhouse gas (GHG) emissions 11,16 . Additionally, over-fertilization can have multiple deleterious effects on the soil, including reductions in soil microbial biomass 17 , pH, base cation availability 18 and notable changes in the composition of bacterial and fungal communities 19 . Economically, excess fertilization of oil palm may reduce profitability due to saturation of the yield curve despite increasing fertilizer application so that the linear increase in fertilizer costs could substantially diminish net returns 20 . This is vital to consider given global reduced availability and consequent price increase of fertilizers 21 . The effects of over-fertilization could be reduced by adjusting fertilization rates to levels that compensate only for the quantity of nutrients exported through harvest and through organic fertilization by putting back palm litter 22 and other processing by-products (for example, palm oil mill effluents and empty fruit bunches) 23 . Such a reduction in fertilization intensity may consequently promote efficient retention and recycling of soil nutrients, decrease soil acidification and the associated need for lime application and generally increase profitability due to fertilizer cost savings and hence present a more ecologically and economically attractive option.
Another important management practice in industrial oil palm plantations is weed control, needed to balance the positive and potential negative effects of the understory vegetation and to facilitate access for management operations. Industrial plantations conventionally employ the use of herbicides for easy and quicker removal of the understory vegetation, but this is often complemented by mechanical weeding, for example, during rainy periods and for removing woody plants. Herbicide weed control results in lower understory vegetation diversity 24 while promoting the invasion of herbicide-resistant weeds 25 . It also reduces the habitat complexity, strongly impacting vegetation-dwelling biodiversity [26][27][28] . A more sustainable alternative could be mechanical weeding, which allows for fast regeneration of the understory vegetation due to preservation of their roots. This consequently leads to increased understory plants cover and diversity which can promote nutrient cycling 29,30 and enhance habitat complexity, supporting animal abundance and diversity [26][27][28] . Nevertheless, these ecologically desirable effects may have adverse economic effects as mechanical weeding requires higher labour input compared with herbicide application, which may reduce profits. Additionally, competition between oil palms and understory vegetation may put yields at risk, particularly during drought periods, given the sensitivity of fructification and fruit bunch weight to water stress 31,32 .
In this study, we report primary data from a multidisciplinary management experiment in a state-owned industrial oil palm plantation in Jambi, Indonesia, that evaluates whether reduced management (that is, reduced fertilization rates and mechanical weeding) as an alternative to conventional management (that is, high fertilization rates and herbicide use) can reduce the negative impacts on ecosystem functions and biodiversity while maintaining high and stable production levels 8 . A full factorial field experiment of two fertilization   Biodiversity, estimated as the average of indicators of taxonomic richness across seven trophic groups, was also higher in the mechanical weeding than herbicide treatment ( Fig. 2 and Table 1). This effect was mainly driven by the strong increase in understory vegetation species richness in the mechanical weeding during four-year measurements (2017-2020; Fig. 2 Table 2). Nevertheless, when removing understory plant species richness from the biodiversity index, there was still a positive marginal effect of mechanical weeding on biodiversity (P = 0.09). Of the 126 understory plant species in the plantation, 33% more plant species occurred in the reduced management than in conventional management. Across the years, the most abundant plant species were the herbicide-resistant invasive shrub Clidemia hirta (L.) D.Don, the invasive herb Asystasia gangetica subsp. micrantha (Nees) Ensermu, the native grass Centotheca lappacea (L.) Desv. and the native fern Christella dentata (Forssk.) Brownsey & Jermy. Compared with herbicide treatment, mechanical weeding increased the ground cover of all these plant species (P < 0.04) except for C. hirta for which ground cover was reduced by 55% (P = 0.01).

Management effects on yield and profit
The four-year cumulative yield did not differ among treatments ( Fig. 3a and Table 1). To assess the notion that reduced fertilization (with or without mechanical weeding) will have less yield stability than conventional fertilization (with or without herbicide treatment), we looked at the left end of the yield distribution. We found that cumulative yield of palms with the lowest performance (that is, the lower fifth quantile of yield per palm per plot) did not differ among treatments (Fig. 3b). Additionally, the probability of the palm yield to fall below 75% of the average yield (that is, yield shortfall probability; Fig. 3c) was similar among treatments, indicating that replacing herbicide with mechanical weeding did not have a negative effect on worst-case yield. When analysing separately for 2017-2018 and 2019-2020, we did not find any effect of treatments on yield and profit indicators ( Supplementary  Fig. 4), signifying that there was no evidence for possible delayed yield reductions in response to reduced management treatment during the four-year measurements. The main management costs considered in the plantation were material (chemical products) costs and labour costs for harvest and weeding operations. Material costs in the reduced fertilization and mechanical weeding treatment were 41% lower (Fig. 3d) than in the conventional fertilization and herbicide treatment (Table 1 and Supplementary Table 2). Conversely, the labour costs were 10% higher in the mechanical weeding than in the herbicide treatment ( Fig. 3d and Supplementary Table 2). The most intensive labour activity was harvesting (accounting for 39-45% of total labour costs) and cleaning of the palm circle (12-14% of total labour costs), which were the same in both weeding treatments. Consequently, profit was 12% higher ( Fig. 3e) in mechanical weeding than herbicide treatment (Table 1 and  Supplementary Table 2). There was an increase in relative gross margin in the reduced fertilization and mechanical weeding by 11% (Fig. 3f) compared with the conventional fertilization and herbicide treatment (Table 1 and Supplementary Table 2).

Discussion
Reduced management (reduced fertilization rates and mechanical weeding) promoted ecosystem functioning and maintained high levels of yield and profit, providing a win-win situation (Fig. 4). In particular, the mechanical removal of weeds instead of using herbicide was found as a tenable management practice that enhanced ecosystem multifunctionality and understory plant diversity, which together with reduced fertilization, increased profit.

Improved multifunctionality and biodiversity
Several ecosystem functions in oil palm plantations tend to be interrelated 23 , which makes it difficult to evaluate the effect of management practices on single functions. Analysing across multiple ecosystem functions, mechanical weeding clearly exhibited higher multifunctionality than herbicide treatment ( Fig. 1 and Supplementary Fig. 3), which was concealed when analysing only for a single ecosystem function ( Supplementary Fig. 1). To note, the indicators of GHG regulation function revealed that this mature plantation located on mineral soils was a slight C source 16 (that is, as indicated by the net ecosystem productivity (NEP) of which fruit harvest was subtracted; Supplementary Table 2). Additionally, it had also large soil N 2 O emissions due to the high N fertilization rate, which is ~50-75% higher than the N fertilizers applied in smallholder oil palm plantations 11 . Nevertheless, the lack of management effects on GHG regulation function or on any single ecosystem function may be due to the relatively short-term span of this management experiment, which was 2016-2020 against at least 16 years of conventional management when this plantation was established (1998)(1999)(2000). The legacy effect of the long-term conventional management practices (Methods) may have dampened the effects of the reduced management practices. Additionally, some of the ecosystem indicators were measured as early as one year after the start of the management experiment (Supplementary Table 1). Indeed, studies that investigated the effects of reduced fertilization and mechanical weeding on root-associated soil biota one year after the start of the experiment found no significant treatment effects, which they attributed to legacy effect of conventional fertilization and herbicide use 33 . Similarly, different weeding treatments did not affect vegetation cover (our indicator for erosion prevention), litter decomposition rates (one of our indicators for organic matter decomposition) and soil physical and biochemical characteristics during the first one to two years after treatment 15,26 . When considering across the eight ecosystem functions, within four years of a shift from herbicide to mechanical weeding, an increased ecosystem multifunctionality was achieved ( Fig. 1 and Supplementary Fig. 3). Mechanical removal of weeds promoted fast vegetation regrowth due to the preservation of the root biomass, which may have resulted in positive feedback effects on plant diversity 24 and ultimately on multifunctionality (Figs. 1 and 2). Experimental studies in mature oil palm plantations that compared herbicide use with manual removal of weeds have reported substantially higher levels of plant biomass and cover in manual weeding 24 , with greater potential to support ecosystem functions 26 . We expect the fast recovery of the understory vegetation under mechanical weeding to provide high organic matter input 34 and a more suitable microclimate for soil microbial and faunal activity 26 , which together promote decomposition 26 and soil nutrient retention 35,36 . Indeed, the positive effect of mechanical weeding on multifunctionality in our study was largely related to its effect on soil functions, such as decomposition, soil fertility and water filtration and on plant refugium functions (Supplementary Figs. 2 and 3).  Table 1).
Article https://doi.org/10.1038/s41893-023-01076-x The positive effect of mechanical weeding on biodiversity was mainly driven by plant diversity effects with a small positive effect on other trophic groups (soil microorganisms, soil invertebrates and aboveground insects; Fig. 2). We found a strong increase in understory vegetation species richness in the mechanical weeding ( Fig. 2), which was similar to other findings from oil palm 24 and Eucalyptus 29 plantations when changing from herbicides to manual weeding. An increase in plant diversity enhances nutrient use efficiency and availability 37 through increased decomposition of diverse plant litter and has been linked to enhanced functional diversity, promoting multifunctionality 38 . High plant diversity is likely to create several complex food webs capable of supporting high macrofaunal abundance, which can, in part, stimulate decomposer and pollinator communities 26,39 . Thus, an increase in plant diversity may promote biodiversity at higher trophic levels due to bottom-up effects, with direct benefits for primary consumers and soil microbes due to increased diversity of substrate 40 .
The method of weeding can also have considerable effects on plant species composition due to the encroachment of non-native, invasive plant species 8,24 . Particularly, the suppressing effect of herbicides on vegetation diversity is well documented as it causes selective pressure on weeds and promotes herbicide-resistant species 41 . Accordingly, the plant community in the herbicide treatment plots was dominated by the herbicide-resistant weed C. hirta, but mechanical weeding was effective in slowing down the spread of this species. This has important implications for sustainable management of oil palm plantations as this highly invasive alien species, in addition to being problematic for weed management in the plantation, can also invade tropical forests and threaten their biodiversity 42 . The reduction of C. hirta in the mechanical weeding was also concomitant with a higher cover of the invasive plant A. gangetica subsp. micrantha, an attractive plant for pollinators 43 , and of some native species (such as the grass C. lappacea and the fern C. dentata), which resulted in a higher plant refugium function ( Supplementary Fig. 2).
Given the four-year span of this management experiment, the positive effects of mechanical weeding on biodiversity and multifunctionality signified that such easily adoptable field practices can reap substantial benefits within a short period. Significant plant diversity effects take some time to propagate through trophic levels, which tend to increase considerably with experimental duration 44 . Therefore, it is possible that the positive effects of mechanical weeding on biodiversity and ecosystem functions could get stronger over time, especially on soil invertebrates and aboveground insects 45 . Employing mechanical weeding during the early stage of oil palm establishment may generate timely benefits on ecosystem multifunctionality and biodiversity, although there is no study so far investigating this in young plantations.

Management effect on yield and profit
Reduced management (reduced fertilization rates and mechanical weeding) was as effective in attaining similar yield as the conventional management but at lower costs (Fig. 3, Table 1 and Supplementary  Fig. 5). The average annual yield of 29.6 Mg ha −1 yr −1 under the reduced management system was very close to the recently published 'attainable yield' of 30.6 Mg ha −1 yr −1 for large-scale oil palm plantations, defined as the yield attained with the adoption of economically optimal inputs and management practices 1 . The similar yields among treatments suggest a more efficient use of applied fertilizers in the reduced management system, particularly on N, as shown by a 56% decrease in dissolved N losses in the reduced fertilization and mechanical weeding compared with conventional management (Supplementary Table 2) 15 . Given the time for sex ratio differentiation and fruit development of oil palm 20 , it is estimated that the effects of management practices on the yield may have a delayed response of approximately two years 46 . Our yield measurements covered four years of this management experiment, and the lack of treatment effect on four-year cumulative yield and similar yields between two separate periods (that is, 2017-2018 and 2019-2020; Supplementary Fig. 4a) suggest that reduced fertilization was a more sustainable alternative to existing conventional fertilization regimes that can lead to higher profitability. This was further reinforced by our findings that the reduced management did not lead to an increased risk of only obtaining very low yields, for example, during a dry period in late 2018 to early 2019 caused by El Niño/Southern Oscillation, as shown by our indicators of yield stability (Fig. 3b,c). The reduced management even tended to show a higher worst-case yield and lower probability to fall under the threshold of 75% of average yield compared with conventional management. It is important to mention that our notion of production risks refers only to yields whereas we disregarded the fluctuations of producer prices for oil palm fruit bunches on the one hand and trends in input prices on the other hand. While declines in producer prices would not lead to any changes in the ranking of the management alternatives, increasing scarcity and consequent price increase of mineral fertilizers might even increase the economic superiority of the reduced fertilization treatments. We found that the reduced fertilization rates would outperform the conventional fertilization rates in terms of profit if fertilizer costs increased by 100% ( Supplementary  Fig. 6). Global fertilizer prices are currently at near-record levels 21 , exceeding those seen during the food and energy crisis in 2008. The high fertilizer prices may last long, especially given the scarcity in P 47 , vastly over-priced K in agriculture 48 and rising consensus on the social costs of excessive N fertilization 49 , which were not monetized here. The high yield under reduced fertilization suggests that this plantation may be over fertilized. However, this statement holds for mature oil palm plantations that have been heavily fertilized since planting, such as the plantation of this study, but may not be true for younger plantations in which the soil has not accumulated nutrients from fertilization and from the decomposition of senesced fronds 22 . Therefore, fertilization rates may need to be adjusted during the life cycle of the plantation, considering the different nutrient accumulations and nutrient requirements over the course of the oil palm life cycle 50 .
Lastly, the higher profit and relative gross margin in reduced compared with conventional management treatments were mainly due to a reduction in material cost and the maintenance of high yield ( Fig. 3 and Supplementary Fig. 5). Due to the low minimum wage in Indonesia (US$184 per month per person), labour cost constituted only a small share of the total management cost (Fig. 3d) and did not translate into low profit in the reduced management treatment. Yet, an increase in minimum wage in Indonesia or higher wages in other production areas would most likely not translate into a sharp difference in labour costs among treatments. This is because the relative change in the labour cost difference among treatments would still be small because the main time-consuming activities, such as harvesting and cleaning of the palm circle, required the same time in all treatments. Our sensitivity analysis demonstrates that the relative differences in profitability appear robust against changes in assumptions on per unit labour costs. Only at an increase of labour costs by 1,250% would the herbicide treatment (with reduced fertilization rates) become more profitable than the respective mechanical weeding ( Supplementary  Fig. 6). Nevertheless, under labour shortage, it needs to be considered that even if the difference in labour hours was small (approximately 12 extra man-hours ha −1 yr −1 ), mechanical weeding increases total labour hours of plantation management.
This management experiment was carried out in an industrial plantation with all works being done by staff members of the plantation within the course of the common plantation management. While we cannot fully exclude any bias, for example, labour time as recorded by researchers observing the operations, we are still convinced that this inherent variation will be smaller than the differences in total management cost, profit and relative gross margin we found between reduced and conventional management treatments (Table 1).
Overall, the results of our four-year management experiment provided early indications that mechanical weeding, together with reduced, compensatory fertilization rates in mature, industrial oil palm plantations, can help in minimizing soil nutrient leaching, decrease water pollution risk, eliminate the effect of herbicide on native vegetation and other non-target soil trophic groups 33 , reduce risks to health of plantation workers and thus contribute to sustainability guidelines of certification bodies such as Roundtable of Sustainable Oil Palm 51 . We acknowledge that when analysing across multiple ecosystem functions, the positive effects of the reduced management on multifunctionality may be still small during the first four years of this experiment as there was no significant effect when analysing for individual ecosystem functions. However, the improved ecosystem multifunctionality under mechanical weeding exemplifies a win-win situation, given its high yield and profit on one hand and the environmental costs associated  Fig. 1), seven multitrophic richness for biodiversity (green; Fig. 2) and six indicators for yield and profit (orange; Fig. 3). 2 2 factorial treatments: + indicates conventional fertilization and herbicide treatment; − denotes reduced fertilization and mechanical weeding.
Article https://doi.org/10.1038/s41893-023-01076-x with herbicide use on the other. Reduced management is therefore a viable management option to maintain optimal yield, lower material cost and improve biodiversity and ecosystem multifunctionality in industrial oil palm plantations that are located on mineral soils.

Ethics
No ethics approval was required for this study. Our study was conducted in a state-owned industrial oil palm plantation where we established a cooperation with the estate owner to access the site and collect data. No endangered or protected species were sampled. Research permits were obtained from the Ministry of Research, Technology and Higher Education, and sample collection and sample export permits were obtained from the Ministry of Environment and Forestry of the Republic of Indonesia.

Study area and experimental design
Our study was conducted in a state-owned industrial oil palm plantation (PTPN VI) located in Jambi, Indonesia (1.719° S, 103.398° E, 73 m above sea level). Initial planting of oil palms within the 2,025 ha plantation area started in 1998 and ended in 2002; planting density was 142 palms ha −1 , spaced 8 m apart in each row and between rows, and palms were ≥16 years old during our study period of 2016-2020. The study sites have a mean annual temperature of 27.0 ± 0.2 °C and a mean annual precipitation of 2,103 ± 445 mm (2008-2017, Sultan Thaha Airport, Jambi). The management practices in large-scale oil palm plantations typically result in three contrasting management zones: (1) a 2 m radius around the base of the palm that was weeded (four times a year) and raked before fertilizer application, hereafter called the 'palm circle'; (2) an area occurring every second inter-row, where pruned senesced palm fronds were piled up, hereafter called 'frond piles'; and (3) the remaining area of the plantation where less weeding (two times a year) and no fertilizer were applied, hereafter called 'inter-rows'. Within this oil palm plantation, we established a management experiment in November 2016 with full factorial treatments of two fertilization rates × two weeding practices: conventional fertilization rates at PTPN VI and other large-scale plantations (260 kg N-50 kg P-220 kg K ha −1 yr −1 ), reduced fertilization rates based on quantified nutrient export by harvest (136 kg N-17 kg P ha −1 yr −1 -187 kg K ha −1 yr −1 ), herbicide and mechanical weeding 15 . The reduced fertilization treatment was based on quantified nutrient export from fruit harvest, calculated by multiplying the nutrient content of fruit bunches with the long-term yield data of the plantation. Fertilizers were applied yearly in April and October following weeding and raking of the palm circle. The common practice at PTPN VI and other large-scale plantations on acidic Acrisol soils is to apply lime and micronutrients, and these were unchanged in our management experiment. Before each N-P-K fertilizer application, dolomite and micronutrients were applied to the palm circle in all treatment plots using the common rates) 52 : 426 kg ha −1 yr −1 dolomite and 142 kg Micro-Mag ha −1 yr −1 (containing 0.5% B 2 O 3 , 0.5% CuO, 0.25% Fe 2 O 3 , 0.15% ZnO, 0.1% MnO and 18% MgO). Herbicide treatment was carried out using glyphosate in the palm circle (1.50 l ha −1 yr −1 , split into four applications per year) and in the inter-rows (0.75 l ha −1 yr −1 , split into two applications per year). Mechanical weeding was done using a brush cutter in the same management zones and at the same frequency as the herbicide treatment.
The 2 2 factorial design resulted in four treatment combinations: conventional fertilization with herbicide treatment, reduced fertilization with herbicide treatment, conventional fertilization with mechanical weeding and reduced fertilization with mechanical weeding. The four treatments were randomly assigned on 50 m ×\ 50 m plots replicated in four blocks, totalling 16 plots. The effective measurement area was the inner 30 m × 30 m area within each replicate plot to avoid any possible edge effects. For indicators (below) that were measured within subplots, these subplots were distributed randomly within the inner 30 m × 30 m of a plot. All replicate plots were located on flat terrain and on an Acrisol soil with a sandy clay loam texture.

Ecosystem functions and multifunctionality
Our study included multiple indicators for each of the eight ecosystem functions 23 , described in details below (Supplementary Tables 1 and 2). All the parameters were expressed at the plot level by taking the means of the subplots (that is, biological parameters) or the area-weighted average of the three management zones per plot (that is, soil parameters). (1) Greenhouse gas (GHG) regulation was indicated by NEP, soil organic C (SOC) and soil GHG fluxes. (2) Erosion prevention was signified by the understory vegetation cover during the four-year measurements. (3) Organic matter decomposition was indicated by leaf litter decomposition and soil animal decomposer activity. (4) Soil fertility was signified by gross N mineralization rate, effective cation exchange capacity (ECEC), base saturation and microbial biomass N. (5) Pollination potential was designated by pan-trapped arthropod abundance and nectar-feeding bird activity. As such, it does not quantify the pollination potential for oil palm, which is mainly pollinated by a single weevil species, but rather as a proxy for a general pollination potential for other co-occurring plants. (6) Water filtration (the capacity to provide clean water) was indicated by leaching losses of the major elements. (7) Plant refugium (the capacity to provide a suitable habitat for plants) as signified by the percentage ground cover of invasive plants to the total ground cover of understory vegetation during the four-year measurements. (8) Biological control (the regulation of herbivores via predation) was indicated by insectivorous bird and bat activities and the soil arthropod predator activity.
All the ecosystem functions were merged into a multifunctionality index using the established average and threshold approaches 12 . For average multifunctionality, we first averaged the z-standardized values (Statistics) of indicators for each ecosystem function and calculated the mean of the eight ecosystem functions for each plot. For threshold multifunctionality, this was calculated from the number of functions that exceeds a set threshold, which is a percentage of the maximum performance level of each function 12 ; we investigated the range of thresholds from 10% to 90% to have a complete overview. The maximum performance was taken as the average of the three highest values for each indicator per ecosystem function across all plots to reduce effect of potential outliers. For each plot, we counted the number of indicators that exceeded a given threshold for each function and divided by the number of indicators for each function 12 .

Indicators of GHG regulation
We calculated annual NEP for each plot as: net ecosystem C exchange - harvested fruit biomass C (ref. 16 ), whereby net ecosystem C exchange = C out (or heterotrophic respiration) - C in (or net primary productivity) 53 . The net primary productivity of oil palms in each plot was the sum of aboveground biomass production (aboveground biomass C + frond litter biomass C input + fruit biomass C) and belowground biomass production. Aboveground biomass production was estimated using allometric equations developed for oil palm plantations in Indonesia 54 , using the height of palms measured yearly from 2019 to 2020. Annual frond litter biomass input was calculated from the number and dry mass of fronds pruned during harvesting events of an entire year in each plot and was averaged for 2019 and 2020. Aboveground biomass production was converted to C based on C concentrations in wood and leaf litter 55 . Annual fruit biomass C production (which is also the harvest export) was calculated from the average annual yield in 2019 and 2020 and the measured C concentrations of fruit bunches. Belowground root biomass and litter C production were taken from previous work in our study area 55 , and it was assumed constant for each plot. Heterotrophic respiration was estimated for each plot as: annual soil CO 2 C emission (below) × 0.7 (based on 30% root respiration contribution to soil respiration from a tropical forest in Sulawesi, Article https://doi.org/10.1038/s41893-023-01076-x Indonesia 56 ) + annual frond litter biomass C input × 0.8 (~80% of frond litter is decomposed within a year in this oil palm plantation 8 ). SOC was measured in March 2018 from composite samples collected from two subplots in each of the three management zones per plot down to 50 cm depth. Soil samples were air dried, finely ground and analysed for SOC using a CN analyser (Vario EL Cube, Elementar Analysis Systems). SOC stocks were calculated using the measured bulk density in each management zone, and values for each plot were the area-weighted average of the three management zones (18% for palm circle, 15% for frond piles and 67% for inter-rows) 15,22 .
From July 2019 to June 2020, we conducted monthly measurements of soil CO 2 , CH 4 and N 2 O fluxes using vented, static chambers permanently installed in the three management zones within two subplots per plot 11,57 . Annual soil CO 2 , CH 4 and N 2 O fluxes were trapezoidal interpolations between measurement periods for the whole year, and values for each plot were the area-weighted average of the three management zones (above).

Indicators of erosion prevention
Diversity and abundance of vascular plants were assessed once a year from 2016 to 2020 before weeding in September-November. In five subplots per plot, we recorded the occurrence of all vascular plant species and estimated the percent cover of the understory vegetation. The percentage cover and plant species richness of each measurement year were expressed in ratio to that of 2016 to account for initial differences among the plots before the start of the experiment. For example, percentage cover in 2017 was: The values from five subplots were averaged to represent each plot.

Indicators of organic matter decomposition
Leaf litter decomposition was determined using litter bags (20 cm × 20 cm with 4 mm mesh size) containing 10 g of dry oil palm leaf litter 8 . Three litter bags per plot were placed on the edge of the frond piles in December 2016. After eight months of incubation in the field, we calculated leaf litter decomposition as the difference between initial litter dry mass and litter dry mass following incubation. Soil animal decomposer activity is described below (Soil arthropods).

Indicators of soil fertility
All these indicators were measured in February-March 2018 in the three management zones within two subplots per plot 22 . Gross N mineralization rate in the soil was measured in the top 5 cm depth on intact soil cores incubated in situ using the 15 N pool dilution technique 58 . ECEC and base saturation were measured in the top 5 cm depth as this is the depth that reacts fast to changes in management 22 . The exchangeable cation concentrations (Ca, Mg, K, Na, Al, Fe, Mn) were determined by percolating the soil with 1 mol l −1 of unbuffered NH 4 Cl, followed by analysis of the percolates using an inductively coupled plasma-atomic emission spectrometer (ICP-AES; iCAP 6300 Duo view ICP Spectrometer, Thermo Fisher Scientific). Base saturation was calculated as the percentage exchangeable bases (Mg, Ca, K and Na) on ECEC. Microbial biomass N was measured from fresh soil samples using the fumigation-extraction method 59 . The values for each plot were the mean of the two subplots that were the area-weighted average of the three management zones (above) 15,22 .

Indicators of general pollination potential
Fluorescent yellow pan traps were used to sample aboveground arthropods (to determine pollinator communities 60 ) in November 2016, September 2017 and June 2018. The traps were attached to a platform at the height of the surrounding vegetation within a 2 × 3 grid centred in the inter-rows of each plot in six clusters of three traps, totalling 18 traps per plot. Traps were exposed in the field for 48 h. We stored all trapped arthropods in 70% ethanol and later counted and identified to order and species level. The abundance of trapped arthropods in 2017 and 2018 were calculated as the ratio to the abundance in 2016 to account for initial differences among the plots before the start of the experiment. The activity of nectar-feeding birds is described below (Birds and bats).

Indicators of water filtration
Element leaching losses were determined from analyses of soil-pore water sampled monthly at 1.5 m depth using suction cup lysimeters (P80 ceramic, maximum pore size 1 μm; CeramTec) over the course of one year (2017-2018) 15 . Lysimeters were installed in the three management zones within two subplots per plot. Dissolved N was analysed using continuous flow injection colorimetry (SEAL Analytical AA3, SEAL Analytical), whereas these other elements were determined using ICP-AES. The values for each plot were the mean of the two subplots that were the area-weighted average of the three management zones 15,22 .

Indicators of plant refugium
In five subplots per plot, the percentage cover and species richness of invasive understory plant species were assessed once a year from 2016 to 2020 before weeding in September-November. We defined invasive species as those plants non-native to Sumatra 61 and among the ten dominant species (excluding oil palm) in the plantation for each year. The percentage cover of invasive understory plant species of each measurement year was expressed in a ratio to that of 2016 to account for initial differences among the plots before the start of the experiment. The values for each plot were represented by the average of five subplots.

Indicators of biological control
The activities of insectivorous birds and bats are described below (Birds and bats). In five subplots per plot, soil invertebrates were collected (Soil arthropods), counted, identified to taxonomic order level and subsequently classified according to their trophic groups that include predators 60 . The values from five subplots were average to represent each plot.

Biodiversity
Biodiversity was measured by the taxonomic richness of seven multitrophic groups, described in details below (Supplementary Tables  1 and 2).

Understory plant species richness
The method is described above (Indicators of erosion prevention), using the number of species as an indicator (Supplementary Table 2).

Soil microorganism richness
This was determined in May 2017 by co-extracting RNA and DNA from three soil cores (5 cm diameter, 7 cm depth) in five subplots per plot 62 . While DNA extraction describes the entire microbial community, RNA represents the active community. The v3-v4 region of the 16S rRNA gene was amplified and sequenced with a MiSeq sequencer (Illumina). Taxonomic classification was done by mapping curated sequences against the SILVA small subunit (SSU) 138 non-redundant (NR) database 63 with the Basic Local Alignment Search Tool (BLASTN) 64 .

Soil arthropod order richness
For determination of soil arthropods, we collected soil samples (16 cm × 16 cm, 5 cm depth) in five subplots per plot in October-November 2017. We extracted the animals from the soil using a heat-gradient extractor 65 , collected them in dimethyleneglycol-water solution (1:1) Article https://doi.org/10.1038/s41893-023-01076-x and stored in 80% ethanol. The extracted animals were counted and identified to taxonomic order level 61 . They were also assigned to the trophic groups decomposers, herbivores and predators based on the predominant food resources recorded in previous reviews and a local study 66,67 . Orders with diverse feeding habits were divided into several feeding groups, for example, Coleoptera were divided into mostly predatory families (Staphylinidae, Carabidae), herbivorous families (for example, Curculionidae) and decomposer families (for example, Tenebrionidae). The total number of individuals per taxonomic group in each subplot was multiplied by the group-specific metabolic rate, which were summed to calculate soil animal decomposer activity. The values from five subplots were average to represent each plot.

Aboveground arthropod order and insect family richness
In addition to the fluorescent yellow pan traps described above (Indicators of general pollination potential), sweep net and Malaise trap samplings were conducted in June 2018, which targeted the general flying and understory dwelling arthropod communities. Sweep net sampling was conducted within the understory vegetation along two 10 m long transects per plot, with ten sweeping strokes performed per transect. In each plot, we installed a single Malaise trap between two randomly chosen palms and exposed it for 24 h. Arthropods were counted, identified to taxonomic order level and the insects to taxonomic family level and values from the three methods were summed to represent each plot.

Birds and bats
Birds and bats passing at each replicate plot were sampled in September 2017 using SM2Bat + sound recorders (Wildlife Acoustics) with two microphones (SMX-II and SMX-US) placed at a height of 1.5 m in the middle of each plot 68 . We assigned the bird vocalization to species with Xeno-Canto 69 and the Macaulay library 70 . Insectivorous bat species richness was computed by dividing them into morphospecies based on the characteristics of their call (call frequency, duration, shape). In addition, we gathered information on proportional diet preferences of the bird species using the EltonTrait database 71 . We defined birds feeding on invertebrates (potential biocontrol agents) as the species with a diet of at least 80% invertebrates and feeding on nectar (potential pollinators), if the diet included at least 20% of nectar.

Economic indicators
We used six indicators linked to the level and stability of yield and profit: yield, lower fifth quantile of the yield per palm per plot, shortfall probability, management costs, profit and relative gross margin. We assessed fruit yield by weighing the harvested fruit bunches from each palm within the inner 30 m × 30 m area of each plot. The harvest followed the schedule and standard practices of the plantation company: each palm was harvested approximately every ten days and the lower fronds were pruned. For each plot, we calculated the average fruit yield per palm and scaled up to a hectare, considering the planting density of 142 palms per ha. Because the palms in each plot have different fruiting cycles and were harvested continuously, the calculation of an annual yield may lead to misleading differences between treatments. Therefore, we calculated the cumulative yield from the beginning of the experiment to four years (2017-2020), which should account for the inter-and intra-annual variations in fruit production of the palms in the plots and thus allowing for comparison among treatments. As effects of management practices on yield may be delayed 46 , we also calculated the cumulative yield during two consecutive years (2017-2018 and 2019-2020) and checked for treatment effects on yield and profit indicators separately for these two periods.
We computed risk indicators on the cumulative yield and on the yield between the two periods. We used the lowest fifth quantile of the yield per palm per plot (left side of the distribution) to indicate the production of the palms with lowest performance. Also, we determined the yield shortfall probability (lower partial moment 0th order), defined as the share of palms that fell below a predefined threshold of yield; the thresholds chosen were 630 kg −1 per palm for cumulative yield and 300 kg −1 per palm per year for the two-year yield, which corresponded to 75% of the average yield.
Revenues and costs were calculated as cumulative values during four years of the experiment (2017-2020) using the same prices and costs for all the years. This was because we were interested in assessing the economic consequences of different management treatments, and they might be difficult to interpret when changes in prices and costs between calendar years are included, which are driven by external market powers rather than the field-management practices. For the same reason, we abstained from discounting profits. Given the usually high discount rates applied to the study area, slight differences in harvesting activities between calendar years or months might lead to high systematic differences between the management treatments, which are associated with the variation in work schedule within the plantation rather than the actual difference among management treatments. Revenues were calculated from the yield and the average price of the fruit bunches in 2016 and 2017 61 . Material costs were the sum of the costs of fertilizers, herbicide and gasoline for the brush cutter. Labour costs were calculated from the minimum wage in Jambi and the time (in labour hours) needed for the harvesting, fertilizing and weeding operations, which were recorded in 2017 for each plot. The weeding labour included the labour for raking the palm circle before fertilization, which was equal in all treatments, and the weeding in the palm circle and inter-rows either with herbicide or brush cutter. In addition, we included the time to remove C. hirta, which must be removed mechanically from all plots once a year, calculated from the average weed-removal time in the palm circle and the percentage cover of C. hirta in each plot for each year. We then calculated the profit as the difference between revenues and the total management costs and the relative gross margin as the gross profit proportion of the revenues.

Statistics
To test for differences among management treatments for each ecosystem function and across indicators of biodiversity, the plot-level value of each indicator was first z standardized (z = (actual value − mean value across plots) / standard deviation) 4 . This prevents the dominance of one or few indicators over the others, and z standardization allows several distinct indicators to best characterize an ecosystem function or biodiversity 4 . Standardized values were inverted (multiplied by −1) for indicators of which high values signify undesirable effect (that is, NEP, soil N 2 O and CH 4 fluxes, element leaching losses, invasive plant cover, yield shortfall, management costs) for intuitive interpretations. For a specific ecosystem function (Supplementary Figs. 1 and 2) and across indicators of biodiversity (Fig. 2), linear mixed-effects (LME) models were used to assess differences among management treatments (fertilization, weeding and their interaction) as fixed effects with replicate plots and indicators (Supplementary Tables 2 and 3) as random effects. The significance of the fixed effects was evaluated using ANOVA 72 . The LME model performance was assessed using diagnostic residual plots 73 . As indicator variables may systematically differ in their responses to management treatments, we also tested the interaction between indicator and treatment (Table 1). For testing the differences among management treatments across ecosystem functions (that is, multifunctionality; Fig. 1), we used for each replicate plot the average of z-standardized indicators of each ecosystem function and ranges of thresholds (that is, number of functions that exceeds a set percentage of the maximum performance of each function 12 ; Supplementary  Fig. 3). The LME models had management treatments (fertilization, weeding and their interaction) as fixed effects and replicate plots and ecosystem functions as random effects; the interaction between ecosystem function and treatment were also tested to assess if there were systematic differences in their responses to management treatments Article https://doi.org/10.1038/s41893-023-01076-x (Table 1). As we expected that the type of weeding will influence ground vegetation, we tested for differences in ground cover of understory vegetation, measured from 2016 to 2020, using LME with management treatments as fixed effect and replicate plots and year as random effects. Differences among management treatments (fertilization, weeding and their interaction) in yield and profit indicators, which were cumulative values over four years (Fig. 3) or for two separate periods (2017-2018 and 2019-2020; Supplementary Fig. 4), were assessed using linear model ANOVA (Table 1). For clear visual comparison among management treatments across ecosystem functions, multitrophic groups for biodiversity, and yield and profit indicators, the fifth and 95th percentiles of their z-standardized values were presented in a petal diagram (Fig. 4 and Supplementary Fig. 5). Data were analysed using R (version 4.0.4), using the R packages 'nlme' and 'influence.ME' 73 .

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Corresponding author(s): Najeeb Al-Amin Iddris Last updated by author(s): Jan 24, 2023 Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Software and code
Policy information about availability of computer code Data collection not applicable. All data were measured in actual field conditions.
For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub All studies must disclose on these points even when the disclosure is negative.

Study description
We assessed the effect of management intensity on ecosystem functions, biodiversity and profitability in a state-owned industrial oil palm plantation in Jambi, Indonesia. The study consisted of a full factorial field experiment with two fertilization levels (conventional vs reduced fertilization rates) and two weeding methods (herbicide vs mechanical weeding), resulting in four treatment combinations. The 8 ecosystem functions and their indicators measured in this study were: 1) Greenhouse gas (GHG) regulation -net ecosystem productivity, soil organic carbon, soil N2O and CH4 fluxes 2) Erosion prevention -understory vegetation cover 3) Organic matter decomposition -leaf litter decomposition, soil animal decomposer activity 4) Soil fertility -gross N mineralization rate, effective cation exchange capacity, base saturation, microbial biomass N 5) Pollination potential -pan-trapped arthropod abundance, nectar-feeding bird activity 6) Water filtration -N, Al, Ca, Mg, K and Na leaching losses 7) Plant refugium -percentage ground cover of invasive plants to the total ground cover of understory vegetation 8) Biological control -insectivorous bird and bat activities, soil arthropod predator activity The 7 indicators of biodiversity were: 1) Understory plant species richness 2) Soil microorganism richness 3) Soil arthropod order richness 4) Aboveground arthropod order richness 5) Aboveground insect family richness 6) Bat species richness 7) Bird species richness The 6 economic indicators were: 1) Oil palm yield 2) Lower 5th quantile of the yield per palm per plot -indicate the production of the palms with lowest performance 3) Yield shortfall probability -indicate the share of palms which fell below a predefined threshold of yield 4) Management costs -material (chemical products) costs, labor costs 5) Profit -difference between revenues and the total management costs 6) Relative gross margin

Sampling strategy
We used a full factorial field experiment of two fertilization rates × two weeding practices, resulting in four treatment combinations that were randomly assigned on 50 m × 50 m plots. The four treatments were replicated in 4 blocks, totaling to 16 plots. Data on ecosystem functions, biodiversity and economic indicators were collected in all the replicate plots at different spatial scales, based on published standard methods. We collected spatially large scale variables such as bird and bat activity over the entire plot. Oil palm yield was restricted to the central area of the plots (inner 30 m × 30 m area) to avoid potential influence by plot edges. Spatially fine scale variables such as plants, litter decomposition, soil microbes and soil arthropods were measured in five 5 m × 5 m subplots in each replicate plot, which were all within the inner 30 m × 30 m area of each plot. Soil data were measured from two subplots located in three distinct management zones (i.e. palm circle, frond piles and inter-rows, which typically results from management activities of large-scale oil palm plantations) per replicate plot.

Data collection
All those above-mentioned indicators (see Research sample) were quantified at each replicate plot (see Sampling strategy) using published standard methods.
Timing and spatial scale The above-mentioned indicators (see Research sample) were measured once in a year or monthly, depending on the approved standard methods of those indicators, spanning from 2017 to 2020. See Sampling strategy above for spatial scale.

Data exclusions
No data were excluded from the analyses.

Reproducibility
Well-established standardized methods have been used for sample collections and analyses, as described in the Methods section. Many of the measured indicators are further measured in this ongoing research (2020-2023) beyond the years covered in the present study (2017-2020).

Randomization
Study plots were all located on relatively flat terrain and each plot had approximately the same age of palms and history of management. Plots were organized into four groups, which were treated as blocks. The four treatment combinations were randomly assigned to a block, resulting in 16 plots.

Blinding
Not applicable -our study is field-based experimental design with actual field data measurements.
Did the study involve field work?

Yes No
Field work, collection and transport

Field conditions
Tropical humid climate. The study sites had a mean annual temperature of 27.0 ± 0.2 •C and a mean annual precipitation of 2103 ± 445 mm, with two precipitation peaks in November and March.

Location
The study was conducted in a state-owned Industrial oil palm plantation in Jambi, Indonesia (1.719° S, 103.398° E, 73 m above sea level).
Access & import/export A cooperation with the estate owner was established to access the site and sampling. We only collected invertebrates that were not protected by Indonesian law, and vertebrates were passively detected without any disturbance. All research groups involved in the experiment obtained research permits from the Ministry of Research, Technology and Higher Education, and sample collection and/ or sample export permits were obtained from the Ministry of Environment and Forestry of the Republic of Indonesia.

Disturbance
There were no special disturbance as our research activities in this large-scale oil palm plantation were only common field activities.
Reporting for specific materials, systems and methods We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

Animals and other organisms
Policy information about studies involving animals; ARRIVE guidelines recommended for reporting animal research

Laboratory animals
The study did not involve laboratory animals.

Wild animals
Birds and bats were passively detected using sound recordings, and not captured. Information on species, sex, strain and age of soil invertebrates has not been collected and considered in the study, being out of the scope. Soil invertebrate data used in this study is based on field-collected soil samples. The samples were collected in the daytime. Animals were transported to the laboratory within 1-2 days after field sampling. Due to specifics of the collection method, all soil invertebrates were killed to acquire data and stored in 80% ethanol. Animals were extracted using a high-gradient Kempson extractor in water:glycol solution. Counting and identification were done under a dissecting microscope.
Field-collected samples Extracted soil athropods from the soil were collected in dimethyleneglycol-water solution (1:1) and stored in 80% ethanol. DNA and RNA extracts were protected from degradation by applying RNAprotect Bacteria Reagent (Qiagen, Hilden, Germany) to the soil samples in a ratio of 1:1 and stored at -80 degrees celcius until analyses. Samples were always kept in shadow to avoid overheating. After the analyses, the animals were kept in collection for further identification, stable isotope and molecular analyses (ongoing work).

Ethics oversight
All research groups involved in the experiment obtained research permits from the Ministry of Research, Technology and Higher