Article

Tree islands enhance biodiversity and

functioning in oil palm landscapes

https://doi.org/10.1038/s41586-023-06086-5 Delphine Clara Zemp1,2,3 ✉, Nathaly Guerrero-Ramirez2, Fabian Brambach2, Kevin Darras4,

Ingo Grass5, Anton Potapov6, Alexander Röll7, Isabelle Arimond4,8, Johannes Ballauff9,
Received: 7 April 2022
Hermann Behling10, Dirk Berkelmann11, Siria Biagioni12, Damayanti Buchori13,14,
Accepted: 14 April 2023 Dylan Craven15, Rolf Daniel11, Oliver Gailing3,16, Florian Ellsäßer7,17, Riko Fardiansah3,18,19,
Nina Hennings20, Bambang Irawan19, Watit Khokthong7,21, Valentyna Krashevska6,
Published online: xx xx xxxx
Alena Krause6, Johanna Kückes7, Kevin Li4, Hendrik Lorenz7, Mark Maraun6,
Open access Miryam Sarah Merk22, Carina C. M. Moura16, Yeni A. Mulyani23, Gustavo B. Paterno2,
Herni Dwinta Pebrianti19, Andrea Polle3,9, Di Ajeng Prameswari6, Lena Sachsenmaier2,24,25,
Check for updates
Stefan Scheu3,6, Dominik Schneider11, Fitta Setiajiati26, Christina Ani Setyaningsih10,
Leti Sundawati26, Teja Tscharntke4, Meike Wollni3,27, Dirk Hölscher3,7 & Holger Kreft2,3

In the United Nations Decade on Ecosystem Restoration1, large knowledge gaps persist
on how to increase biodiversity and ecosystem functioning in cash crop-dominated
tropical landscapes2. Here, we present findings from a large-scale, 5-year ecosystem
restoration experiment in an oil palm landscape enriched with 52 tree islands,
encompassing assessments of ten indicators of biodiversity and 19 indicators of
ecosystem functioning. Overall, indicators of biodiversity and ecosystem functioning,
as well as multidiversity and ecosystem multifunctionality, were higher in tree islands
compared to conventionally managed oil palm. Larger tree islands led to larger gains in
multidiversity through changes in vegetation structure. Furthermore, tree enrichment
did not decrease landscape-scale oil palm yield. Our results demonstrate that enriching
oil palm-dominated landscapes with tree islands is a promising ecological restoration
strategy, yet should not replace the protection of remaining forests.

The loss of megadiverse tropical lowland rainforests has accelerated alters the functioning of ecological communities and environmental
in the past decades3, with deforestation and land-use change being conditions, leading to a reduction of several ecosystem functions and
largely driven by the rapid expansion of high-yielding cash crops such services7,8.
as oil palm4. Globally, oil palm plantations occupy 21 million hectares, Many agricultural landscapes are in urgent need of ecological res-
mostly in Indonesia and Malaysia5. Although the expansion of oil palm toration to safeguard biodiversity and ecosystem functioning while
has promoted economic development and improved livelihoods of also promoting local livelihoods9–11, a central goal of the current
smallholder farmers, it has also led to dramatic negative ecological United Nations decade on Ecosystem Restoration. However, trade-offs
impacts6. Compared with tropical lowland rainforests, species diversity between biodiversity or ecosystem functioning and agricultural pro-
in oil palm-dominated landscapes is greatly reduced7, especially for ductivity may result in failed restoration efforts or lead to undesir-
forest-dependent species and species of conservation concern4. In addi- able ecological spillover effects by promoting the expansion of the
tion, the transformation of forests to oil palm-dominated landscapes agricultural frontier into natural forested areas12. One way to mitigate
1
Conservation Biology, Institute of Biology, Faculty of Sciences, University of Neuchâtel, Neuchâtel, Switzerland. 2Biodiversity, Macroecology and Biogeography, Faculty of Forest Sciences and
Forest Ecology, University of Göttingen, Göttingen, Germany. 3Centre of Biodiversity and Sustainable Land Use (CBL), University of Göttingen, Göttingen, Germany. 4Agroecology, Department
of Crop Sciences, Faculty of Agricultural Science, University of Göttingen, Göttingen, Germany. 5Ecology of Tropical Agricultural Systems, Institute of Agricultural Sciences in the Tropics,
University of Hohenheim, Stuttgart, Germany. 6Animal Ecology, J.F. Blumenbach Institute of Zoology and Anthropology, University of Göttingen, Göttingen, Germany. 7Tropical Silviculture and
Forest Ecology, Faculty of Forest Sciences and Forest Ecology, University of Göttingen, Göttingen, Germany. 8Functional Agrobiodiversity, Dept. of Crop Sciences, Faculty of Agricultural
Science, University of Göttingen, Göttingen, Germany. 9Forest Botany and Tree Physiology, Faculty of Forest Sciences and Forest Ecology, University of Göttingen, Göttingen, Germany.
10
Department of Palynology and Climate Dynamics, Albrecht-von-Haller-Institute for Plant Sciences, Göttingen, Germany. 11Department of Genomic and Applied Microbiology, Institute of
Microbiology and Genetics, University of Göttingen, Göttingen, Germany. 12Department of Palynology and Climate Dynamics, Albrecht-von-Haller-Institute for Plant Sciences, University of
Göttingen, Göttingen, Germany. 13Department of Plant Protection, Faculty of Agriculture, Institut Pertanian Bogor. Jl. Meranti, IPB Dramaga Campus, Bogor, Indonesia. 14Center for
Transdisciplinary and Sustainability Sciences, IPB University, Jalan Pajajaran, Indonesia. 15Centre for Ecosystem Modeling and Monitoring, Facultad de Ciencias, Universidad Mayor, Santiago,
Chile. 16Forest Genetics and Forest Tree Breeding, Faculty of Forest Sciences and Forest Ecology, University of Göttingen, Göttingen, Germany. 17Department of Natural Resources, University of
Twente, Enschede, Netherlands. 18Zoological Museum, Center of Natural History, Universität Hamburg, Hamburg, Germany. 19Faculty of Forestry, University of Jambi Jln Raya Jambi, Jambi,
Indonesia. 20Biogeochemistry of Agroecosystems, Faculty of Agricultural Science, University of Göttingen, Göttingen, Germany. 21Department of Biology, Faculty of Science, Chiang Mai
University, Chiang Mai, Thailand. 22Chairs of Statistics and Econometrics, Faculty of Business and Economics, University of Göttingen, Göttingen, Germany. 23Forest Resources Conservation
and Ecotourism, Faculty of Forestry and Environment, IPB University, Kampus IPB Darmaga, Bogor, Indonesia. 24German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig,
Leipzig, Germany. 25Systematic Botany and Functional Biodiversity, Institute of Biology, Leipzig University, Leipzig, Germany. 26Department of Forest Management, Faculty of Forestry and
Environment, IPB University, Kampus IPB Darmaga, Bogor, Indonesia. 27Environmental and Resource Economics, Department of Agricultural Economics and Rural Development, Faculty of
Agricultural Sciences, University of Göttingen, Göttingen, Germany. ✉e-mail: [email protected]

Nature | www.nature.com | 1
Article
Planted diversity
(number of tree species)

0 1 2 3 6 Control

Island area:
edge length (m)

5 10 20 40

0 200 400 m
5m

Standing oil palm

Felled oil palm

10 m
Planted trees
Parkia speciosa, Dyera polyphylla,
Fabacea Apocynaceae
Archidendron jiringa, Peronema canescens,
Fabaceae Lamiaceae
Durio zibethinus, Shorea leprosula,
Malvaceae Dipterocarpaceae

20 m

Fig. 1 | Experimental design that tests the ecological restoration outcomes a total of 52 tree islands established in an industrial oil palm plantation in
of tree island establishment in oil palm-dominated landscapes. Tree islands Sumatra, Indonesia. Control plots represent conventionally managed oil palm
vary in area (25–1,600 m2) and planted tree diversity (none to six species), with monocultures. Note that the islands in the map are not at scale.

trade-offs between restoration outcomes is to enrich agricultural land- overview of biodiversity and ecosystem functioning across the experi-
scapes with species-rich agroforestry systems13,14 and islands of native ment, we calculated multidiversity and multifunctionality using the
trees through planting or natural regeneration15–17. However, to be a aforementioned indicators18. The experimental design allowed us to
viable alternative for landowners, it is essential to generate empiri- test the effects of tree island area (25, 100, 400 and 1,600 m2) and of
cal evidence on whether and how these restoration strategies affect planted native tree diversity (zero, one, two, three and six species, with
biodiversity, ecosystem functioning and agricultural productivity in zero representing natural regeneration only) on restoration outcomes
cash crop-dominated landscapes2. and to compare them with conventionally managed oil palm monocul-
Here, we present the results of a large-scale, interdisciplinary ecosys- tures19 (Fig. 1). Overall, we expected tree islands to enhance biodiversity
tem restoration experiment, in which the restoration outcomes across and ecosystem functioning compared to conventionally managed oil
52 tree islands established in a landscape dominated by an industrial palm monocultures. To provide a mechanistic understanding of the
oil palm plantation (140 ha) were observed and quantified three to effects of planted tree diversity and island area on biodiversity and
five years after establishment. We assessed above- and below-ground ecosystem functioning, we also measured 12 indicators of vegetation
biodiversity across ten indicators representing a broad range of King- structure (Supplementary Table 3). On the basis of the theory of island
doms (bacteria, fungi, plants and animals; Supplementary Table 1) and biogeography20, we expected larger tree islands to have enhanced
19 indicators of ecosystem functioning associated with productivity biodiversity and ecosystem functioning compared to smaller ones.
of oil palms and planted trees, resistance to invasion, pollination, soil Larger tree islands potentially provide more habitats and sustain larger
quality, predation and herbivory, carbon and nutrient cycling and water populations, whereas smaller islands are expected to be more like the
and climate regulation (Supplementary Table 2). To provide an holistic surrounding environment, that is, the oil palm-dominated landscape.

2 | Nature | www.nature.com
 a Species richness b Simpson diversity c Ecosystem functioning

1 2
19
18 3
17 4
16 5

15 6

14 7
13 8
12 9
11 10
Tree islands
Controls

Bat Bird Arthropod Herb Pollen Seed Tree Soil Soil Soil Productivity Resistance Pollination Soil Predation Nutrient Water
fauna bacteria fungi to invasion quality and and carbon and climate
herbivory cycling regulation

d Multidiversity (species richness) e Multidiversity (Simpson diversity) f Multifunctionality

20
10.0 10.0
Tree islands Controls
Quartiles Quartiles
15
7.5 7.5 Median Median
Number of indicators

Individual plot Individual plot

5.0 10
5.0

2.5 2.5 5

0 0 0

0 25 50 75 100 0 25 50 75 100 0 25 50 75 100

Threshold (percentage of maximum observed values) Threshold (percentage of maximum observed values) Threshold (percentage of maximum observed values)

Fig. 2 | Multidimensional ecological restoration outcomes in an oil palm- island) is considered as an ecosystem functioning because of its contribution
dominated landscape. We measured 10 and 19 indicators of biodiversity and to primary productivity, as well as agricultural productivity. a–c, Indicators of
ecosystem functioning, respectively, in tree islands and compared their biodiversity calculated as species richness (a) and Simpson diversity (b), which
responses to those in plots representing conventionally managed oil palm emphasizes the contribution of abundant species and ecosystem functioning
monocultures. For ecosystem functioning, we measured: productivity as (1) oil (c) across 52 tree islands (green polygons) compared to four control plots of
palm yield and (2) above-ground biomass; resistance to invasion of (3) native conventionally managed oil palm monocultures (grey polygons). d–f, Polygon
seeds and (4) resistance to invasive plants; pollination as (5) pollinators and (6) vertices represent median values for each indicator. Multidiversity and
pollination rate; soil quality as (7) soil P, (8) soil decompaction and (9) 1/soil C:N; multifunctionality represent the number of indicators (species richness (d);
predation and herbivory as (10) predators (vertebrates), (11) predators Simpson diversity (e) and ecosystem functioning (f)) that exceed a specified
(arthropods), (12) predators (soil fauna) and (13) herbivores (soil fauna); carbon threshold, which is expressed as a percentage of the maximum observed values
and nutrient cycling as (14) decomposers, (15) litter decomposition and (16) in the oil-dominated landscape (calculated on the basis of both island and
litter input; water and climate regulation as (17) evapotranspiration, (18) water control plots combined).
infiltration and (19) microclimate buffering. Oil palm yield (calculated per

We further expected greater planted tree species diversity to favour increases in tree and bird species richness (+4.7 tree species in islands
diversity at higher trophic levels21 and enhance ecosystem functioning compared to monocultures,+2.5 bird species) and decreases in the
through complementarity among species22. Planted diversity effects diversity of the most abundant seed species (−1.2 seed species based on
on restoration outcomes are probably mediated by higher vegetation Simpson diversity; Supplementary Table 4). Overall, restoration ben-
structural complexity, that is, the three-dimensional distribution of efits of tree islands were found for ecosystem functioning (tree island:
plants within an ecosystem23. Finally, we proposed that agricultural F = 6.2, P = 0.016; Extended Data Table 1); with strongest increases for
productivity (oil palm yield) decreases at the local scale (within tree water infiltration (+174% saturated soil hydraulic conductivity), litter
islands), whereas the loss is negligible at the scale of the industrial input (+151% leaf litter biomass input), activity of insectivorous bats and
plantation or landscape16. birds (+556%) and soil fertility (+14% 1/soil C:N ratio; Supplementary
Tree islands had higher biodiversity and ecosystem functioning Table 5). Overall, multidiversity and ecosystem multifunctionality
compared to conventionally managed oil palm monocultures (Fig. 2 were higher in tree islands than in conventionally managed oil palm
and Extended Data Table 1). Yet, tree island effects on biodiversity varied monocultures, regardless of the threshold used for calculation or when
depending on the indicator (tree island × indicator: F = 2.5, P = 0.007 for considering relative species abundances (Fig. 2d–f). The calculation
species richness; F = 3.6, P  0.0002 for Shannon diversity; and F = 3.0, of multidiversity (or multifunctionality) relies on the number of biodi-
P =0.001 for Simpson diversity; Extended Data Table 1 and Extended versity (or functioning) indicators that cross a certain threshold, with
Data Fig. 1). For example, natural regeneration and colonization led to thresholds expressed as the percentage of the maximum observed

Nature | www.nature.com | 3
Article
a
Structural complexity (PC1)

Tree diversity 0.32

Open Dense and complex

Multidiversity
Species richness (R2 = 17%)
Tree dominance (- PC2)
7 Simpson diversity (R2 = 4%)
Tree island area 0.3
8 Effect of tree diversity
0.62 0.3
Direct
Indirect
Effect of tree island size
Oil palm dominance Tree dominance
Direct
Indirect

b Species richness c Simpson diversity d Ecosystem functioning

R2 Standard estimate R2 Standard estimate R2 Standard estimate

(%) –0.25 0 0.25 (%) –0.25 0 0.25 (%) –0.25 0 0.25
13 Oil palm yield
8 17
21 Aboveground biomass
4 3 Native seeds
5
2 Resistance to invasive plants
52 21 38 Pollinators
1 Pollination rate
8 13 15 Soil P
21 Soil decompaction
1 4 3 1/soil C:N
21 Predators (vertebrates)
30 7 5 Predators (arthropods)
0 Predators (soil fauna)
43 28 4 Herbivores (soil fauna)
0 Decomposers
1 19 0 Litter decomposition
48 Litter input
5 11
10 Evapotranspiration
0 21 Water infiltration
1
53 Microclimate buīering

Fig. 3 | Influence of tree island area and planted tree diversity on quantified by species richness (b), Simpson diversity (c) and ecosystem
multidimensional restoration outcomes in an oil palm-dominated functioning (d); direct effects of planted tree diversity are indicated by dark
landscape. a, Effects of planted tree diversity (directly or through structural green bars and indirect effects through structural complexity are indicated by
complexity) and tree island area (directly or through changes in tree dominance) light green bars and the direct effects of tree island area are indicated by black
on multidiversity and multifunctionality tested with SEMs. Filled arrows bars and indirect effects through tree dominance are indicated by grey bars.
(and standard coefficient estimates) indicate statistically significant effects The proportion of explained variance is shown to the left of b and c. All the bars
(P < 0.05; two-sided analysis of variance (ANOVA), n = 52 tree islands) and dashed indicate significant effects (P < 0.05; two-sided ANOVA, n = 52 tree islands).
arrows indicate non-significant effects. Percentage values indicate explained The legend for icons is presented in Fig. 2.
variance of each endogenous variable. b–d, Effects on indicators of biodiversity

values within the study system18 (here, within our landscape combin- diversity; Extended Data Fig. 1). Our structural equation models (SEMs)
ing islands and conventionally managed oil palm monocultures). For revealed that the influence of the island area acted through changes
example, at the 50% threshold, multidiversity increases by 1.5 in islands in tree dominance (Fig. 3a, Extended Data Table 2 and Supplementary
compared to conventionally managed oil palm monocultures. In other Tables 6 and 7), with higher multidiversity in larger tree islands that are
words, four and 2.5 biodiversity indicators reached at least 50% of their dominated by trees rather than oil palms. The higher tree dominance in
maximum observed species richness in tree islands and conventionally the canopy and the thicker leaf litter (Extended Data Fig. 2 and Extended
managed oil palm monocultures, respectively (Supplementary Table 4). Data Table 3) might provide habitats of sufficient quality and quantity
Similarly, six and three ecosystem functioning indicators reached at to enhance multidiversity. Large tree islands may thus act as keystone
least 50% of their maximum observed values in tree islands and conven- structures28 in oil palm-dominated landscapes that facilitate the arrival
tionally managed oil palm monocultures, respectively (Supplementary of seed rain (especially of locally rarer species, see Fig. 3b) and the colo-
Table 5). Overall, our results provide evidence of multidimensional nization, establishment and maintenance of diverse communities, such
ecological restoration benefits with tree islands in oil palm-dominated as understorey arthropods and trees (Fig. 3b,c). Although multifunc-
landscapes. Although the main priority is the protection of the remain- tionality also increased with island area and the effect was mediated
ing tropical forests24, ecological restoration with tree islands along by changes in tree dominance, the strength of the effect depended on
with other practices7,25 and riparian buffer management26,27 plays an the method used to calculate multifunctionality (Extended Data Fig. 3,
essential and complementary role in safeguarding biodiversity and Extended Data Table 2 and Supplementary Table 7). When calculated
ecosystem functioning in cash crop-dominated landscapes. for individual functions, large tree islands were pivotal for provid-
Confirming our initial hypothesis, larger tree islands resulted in ing ecosystem functions related to predation and herbivory (through
greater restoration benefits, both for ecosystem functioning (island predatory arthropods and soil herbivores) and carbon and nutrient
area: F = 12.9, P < 0.0001; Extended Data Table 1) and biodiversity. cycling (through decomposers; Fig 3d and Extended Data Table 3). By
Yet, the effects of island area on biodiversity varied across indica- using constant sampling area or rarefaction curves, we could rule out
tors (island area × indicator: F = 5.1, P < 0.0001 for species richness; that the influence of island area on biodiversity was limited to passive
F = 2.8, P = 0.003 for Shannon diversity; and F = 1.8, P = 0.06 for Simpson sampling29 (Methods). Thus, ecological mechanisms associated with

4 | Nature | www.nature.com
environmental filtering such as reduced edge effects and greater envi- planting tree islands. To enhance the establishment of tree islands in
ronmental heterogeneity probably explain the positive effects of larger oil palm- and other cash crop-dominated landscapes, they could be
tree islands on multidimensional restoration benefits. incorporated as a requirement in existing sustainability certifications
The effect of planted tree diversity on biodiversity—when consid- (for example, Roundtable for Sustainable Palm Oil), alongside other
ering abundances—depended on the biodiversity indicator (planted practices including optimized management25, ecological intensifica-
diversity × indicator: F = 2.3, P = 0.014 for Shannon diversity; F = 2.8, tion and riparian restoration26,36,37. Enhancing the status of sustainability
P = 0.004 for Simpson diversity; Extended Data Table 1). For example, certifications should, however, not come at the expense of smallholder
planted tree diversity promoted the diversity of (non-planted) trees farmers, who are often excluded from certification programmes38, nor
but decreased the diversity of arthropods (Extended Data Fig. 1). The at the expense of the protection of remaining intact forests for their
statistically non-significant effects of multidiversity are probably due to exceptional value as refugia for biodiversity and providers of ecosystem
contrasting responses of biodiversity indicators to planted tree diver- functioning39. Overall, we provide robust evidence that biodiversity and
sity, with contrasting responses mediated by vegetation structure as ecosystem functioning in cash crop-dominated tropical landscapes can
shown by the SEM (Fig. 3b–d and Extended Data Table 2). Higher planted be enhanced without compromising overall agricultural productivity
tree diversity led to structurally more complex habitats30 (Extended by planting tree islands. Although our study was conducted in a single
Data Fig. 2) that benefited some biodiversity indicators (Shannon and landscape, it adds to growing experimental37,40,41 and modelling evi-
Simpson diversity of bats and herbs), whereas others benefited from dence42 of the ecological and economic benefits in oil palm agroforestry
more open and structurally simpler habitats (for example, species rich- systems. Understanding how biodiversity and ecosystem functioning
ness of seeds and understorey arthropods; Fig. 3b,c and Supplementary change in several landscapes43,44 is urgently needed for designing and
Table 6). More open habitats also favoured native seeds, pollinators scaling-up ecological restoration of oil palm landscapes worldwide.
and predators (arthropods) and soil herbivores; whereas structural
complexity enhanced above-ground biomass, litter input and microcli-
mate buffering (Fig. 3d and Supplementary Table 7). Through changes Online content
in structural complexity, tree diversity had a negative impact on mul- Any methods, additional references, Nature Portfolio reporting summa-
tifunctionality, although the strength of the effect depended on the ries, source data, extended data, supplementary information, acknowl-
methods of calculation (Extended Data Fig. 3, Extended Data Table 2 edgements, peer review information; details of author contributions
and Supplementary Table 7). Owing to specific responses associated and competing interests; and statements of data and code availability
with the adaptability of different organisms to contrasting habitats31, are available at https://doi.org/10.1038/s41586-023-06086-5.
establishing a combination of tree islands that differ in structural
complexity may favour differences in local community composition 1. UN Decade on Ecosystem Restoration (UN, 2021); https://www.decadeonrestoration.org/
leading to increases of gamma diversity and ecosystem functioning about-un-decade.
at the landscape scale32. 2. Pashkevich, M. D. et al. Nine actions to successfully restore tropical agroecosystems.
Trends Ecol. Evol. 37, 963–975 (2022).
Our study shows that the magnitude of the effect of ecological res- 3. Lewis, S. L., Edwards, D. P. & Galbraith, D. Increasing human dominance of tropical forests.
toration on oil palm yield depends strongly on the spatial scale, with Science 349, 827–832 (2015).
4. Meijaard, E. et al. The environmental impacts of palm oil in context. Nat. Plants 6,
declines at the local scale, that is, within tree islands but no statisti-
1418–1426 (2020).
cally significant reduction at the landscape scale. At the local scale, 5. Descals, A. et al. High-resolution global map of smallholder and industrial closed-canopy
per area yields within tree islands were on average 24% lower than in oil palm plantations. Earth Syst. Sci. Data 13, 1211–1231 (2021).
6. Qaim, M., Sibhatu, K. T., Siregar, H. & Grass, I. Environmental, economic, and social
the conventionally managed oil palm monocultures (Extended Data
consequences of the oil palm boom. Annu. Rev. Res. Econ. 12, 321–344 (2020).
Fig. 4) because of the removal of oil palms, which reduced palm density 7. Grass, I. et al. Trade-offs between multifunctionality and profit in tropical smallholder
(Extended Data Fig. 5b). In contrast, no statistically significant differ- landscapes. Nat. Commun. 11, 1186 (2020).
8. Dislich, C. et al. A review of the ecosystem functions in oil palm plantations, using forests
ence was detected in per island yield when including the yield of palms
as a reference system. Biol. Rev. 92, 1539–1569 (2017).
adjacent to tree islands (Fig. 2). The yield gains per oil palm surrounding 9. The IPBES Assessment Report on Land Degradation and Restoration (IPBES, 2018).
the tree islands thus compensated for yield losses per area within the 10. Chazdon, R. & Brancalion, P. Restoring forests as a means to many ends. Science 365,
24–25 (2019).
islands, with these beneficial effects resulting from oil palm thinning in
11. Brancalion, P. H. S. et al. Global restoration opportunities in tropical rainforest landscapes.
the tree islands (Extended Data Fig. 5a,e). These beneficial effects were Sci. Adv. 5, eaav3223 (2019).
already observed a few years after establishment of the experiment33 12. Coleman, E. A. et al. Limited effects of tree planting on forest canopy cover and rural
livelihoods in Northern India. Nat. Sustain. 4, 997–1004 (2021).
and are consistent regardless of the time period considered (Supple-
13. Kremen, C. & Merenlender, A. M. Landscapes that work for biodiversity and people.
mentary Note 1). Over time, yield decreases within tree islands are Science 362, eaau6020 (2018).
expected because of competition with trees, particularly in tree islands 14. Arroyo‐Rodríguez, V. et al. Designing optimal human-modified landscapes for forest
biodiversity conservation. Ecol. Lett. 23, 1404–1420 (2020).
with higher planted diversity (Fig. 3d and Extended Data Fig. 4). Yet, 15. Corbin, J. D. & Holl, K. D. Applied nucleation as a forest restoration strategy. For. Ecol.
these effects will remain negligible on industrial large-scale plantations Manag. 265, 37–46 (2012).
because of the relatively small area covered by the tree islands. In our 16. Rey Benayas, J. M., Bullock, J. M. & Newton, A. C. Creating woodland islets to reconcile
ecological restoration, conservation, and agricultural land use. Front. Ecol. Environ. 6,
experiment, tree islands covered only 2.8 ha, less than 5% of the 140 ha 329–336 (2008).
industrial oil palm plantation. In contrast, smallholder oil palm planta- 17. Holl, K. D. et al. Applied nucleation facilitates tropical forest recovery: lessons learned
tions often only comprise a few hectares5, of which (larger) tree islands from a 15-year study. J. Appl. Ecol. 57, 2316–2328 (2020).
18. Manning, P. et al. Redefining ecosystem multifunctionality. Nat. Ecol. Evol. 2, 427–436 (2018).
would cover a more substantial proportion. In these cases, a decrease 19. Teuscher, M. et al. Experimental biodiversity enrichment in oil-palm-dominated
in yield may be compensated by extra goods from the tree islands, for landscapes in Indonesia. Front. Plant Sci. 7, 1538 (2016).
example, fruits, natural latex, timber and firewood34,35. Furthermore, 20. Wilson, E. O. & MacArthur, R. H. The Theory of Island Biogeography (Princeton Univ. Press,
1967).
smallholders could benefit from higher levels of several ecosystem 21. Scherber, C. et al. Bottom-up effects of plant diversity on multitrophic interactions in a
services, lower susceptibility to disturbance and risk diversification17. biodiversity experiment. Nature 468, 553–556 (2010).
From an economic perspective, oil palm represents a highly profit- 22. Barry, K. E. et al. The future of complementarity: disentangling causes from consequences.
Trends Ecol. Evol. 34, 167–180 (2019).
able cash crop6. Consequently, replacing oil palms with native tree 23. Coverdale, T. C. & Davies, A. B. Unravelling the relationship between plant diversity and
species typically raises concerns about high opportunity costs of lost vegetation structural complexity: a review and theoretical framework. J. Ecol. https://doi.
revenue among landholders. Our large-scale study offers unique empiri- org/10.1111/1365-2745.14068 (2023).
24. Poorter, L. et al. Multidimensional tropical forest recovery. Science 374, 1370–1376 (2021).
cal evidence on the viability of multidimensional ecological benefits 25. Iddris, N. A.-A. et al. Mechanical weeding enhances ecosystem multifunctionality and profit
without compromising yield in oil palm-dominated landscapes by in industrial oil palm. Nat. Sustain. https://doi.org/10.1038/s41893-023-01076-x (2023).

Nature | www.nature.com | 5
Article
26. Luke, S. H. et al. Riparian buffers in tropical agriculture: scientific support, effectiveness 39. Watson, J. E. M. et al. The exceptional value of intact forest ecosystems. Nat. Ecol. Evol. 2,
and directions for policy. J. Appl. Ecol. 56, 85–92 (2019). 599–610 (2018).
27. Bicknell, J. E. et al. Enhancing the ecological value of oil palm agriculture through 40. Ahirwal, J., Sahoo, U. K., Thangjam, U. & Thong, P. Oil palm agroforestry enhances crop
set-asides. Nat. Sustain. https://doi.org/10.1038/s41893-022-01049-6 (2023). yield and ecosystem carbon stock in northeast India: implications for the United Nations
28. Tews, J. et al. Animal species diversity driven by habitat heterogeneity/diversity: the sustainable development goals. Sustain. Prod. Consum. 30, 478–487 (2022).
importance of keystone structures. J. Biogeogr. 31, 79–92 (2004). 41. de Carvalho, W. R., Vasconcelos, S. S., Kato, O. R., Capela, C. J. B. & Castellani, D. C.
29. Chase, J. M., Blowes, S. A., Knight, T. M., Gerstner, K. & May, F. Ecosystem decay Short-term changes in the soil carbon stocks of young oil palm-based agroforestry
exacerbates biodiversity loss with habitat loss. Nature 584, 238–243 (2020). systems in the eastern Amazon. Agrofor. Syst. 88, 357–368 (2014).
30. Zemp, D. C. et al. Mixed-species tree plantings enhance structural complexity in oil palm 42. Khasanah, N. et al. Oil palm agroforestry can achieve economic and environmental gains
plantations. Agric. Ecosyst. Environ. 283, 106564 (2019). as indicated by multifunctional land equivalent ratios. Front. Sustain. Food Syst. https://
31. Schall, P. et al. The impact of even-aged and uneven-aged forest management on regional doi.org/10.3389/fsufs.2019.00122 (2020).
biodiversity of multiple taxa in European beech forests. J. Appl. Ecol. 55, 267–278 (2018). 43. Oehri, J., Schmid, B., Schaepman-Strub, G. & Niklaus, P. A. Terrestrial land-cover type
32. Montoya-Sánchez, V. et al. Landscape heterogeneity and soil biota are central to richness is positively linked to landscape-level functioning. Nat. Commun. 11, 154 (2020).
multi-taxa diversity for landscape restoration. Preprint at bioRxiv https://doi.org/10.1101/ 44. Fahrig, L. et al. Resolving the SLOSS dilemma for biodiversity conservation: a research
2022.10.31.514517 (2022). agenda. Biol. Rev. 97, 99–114 (2022).
33. Gérard, A. et al. Oil-palm yields in diversified plantations: initial results from a biodiversity
enrichment experiment in Sumatra, Indonesia. Agric. Ecosyst. Environ. 240, 253–260 Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
(2017). published maps and institutional affiliations.
34. Rahmani, T. A. et al. Evaluating the feasibility of oil palm agroforestry in Harapan Rainforest,
Jambi, Indonesia. For. Soc. 5, 458–477 (2021). Open Access This article is licensed under a Creative Commons Attribution
35. Zemp, D. C. et al. Tree performance in a biodiversity enrichment experiment in an oil palm 4.0 International License, which permits use, sharing, adaptation, distribution
landscape. J. Appl. Ecol. 56, 2340–2352 (2019). and reproduction in any medium or format, as long as you give appropriate
36. Darras, K. F. A. et al. Reducing fertilizer and avoiding herbicides in oil palm plantations— credit to the original author(s) and the source, provide a link to the Creative Commons licence,
ecological and economic valuations. Front. For. Glob. Change https://doi.org/10.3389/ and indicate if changes were made. The images or other third party material in this article are
ffgc.2019.00065 (2019). included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
37. Luke, S. H. et al. Managing oil palm plantations more sustainably: large-scale experiments to the material. If material is not included in the article’s Creative Commons licence and your
within the Biodiversity and Ecosystem Function in Tropical Agriculture (BEFTA) intended use is not permitted by statutory regulation or exceeds the permitted use, you will
Programme. Front. For. Glob. Change https://doi.org/10.3389/ffgc.2019.00075 (2020). need to obtain permission directly from the copyright holder. To view a copy of this licence,
38. Azhar, B. et al. Promoting landscape heterogeneity to improve the biodiversity benefits of visit http://creativecommons.org/licenses/by/4.0/.
certified palm oil production: evidence from Peninsular Malaysia. Glob. Ecol. Conserv. 3,
553–561 (2015). © The Author(s) 2023

6 | Nature | www.nature.com
Methods at sunrise, sampled at 22.05 kHz for birds. We used two 40 min mono
sound recordings from the right channel, extracted from consecu-
The biodiversity enrichment experiment tive nights, starting 20 min after sunset, sampled at 384 kHz for bats.
Our study was conducted in EFForTS-BEE, the biodiversity enrichment Twelve sound recorders were installed simultaneously in 12 randomly
experiment of the EFForTS project (Ecological and Socioeconomic chosen plots. The recordings were annotated in ecoSound-web50 to
Functions of Tropical Lowland Rainforest Transformation Systems extract the duration of each bird vocalization and bat pass and bird
(Sumatra, Indonesia))19. EFForTS-BEE is part of the global network of tree detection distances were estimated using reference sound transmis-
diversity experiments TreeDivNet45 (https://treedivnet.ugent.be/). The sion sequences51. We assigned birds to species according to Birdlife
study region is characterized by a humid tropical climate with a mean International taxonomy. Owing to the lack of standard protocols and
temperature of 26.7 ± 0.2 °C and an annual rainfall of 2,235 ± 381 mm reference collections for Southeast Asia, we could not identify bats to
and the dominant soil type is loamy Acrisol46. In December 2013, 52 species and used sonotypes instead. We appended feeding guild infor-
experimental plots (tree islands) were established in a conventionally mation to each bird species52 all detectable bats were echolocating and
managed 140 ha oil palm plantation. Following a random partition thus considered insectivorous bats. Only bird vocalizations detected
design47, we systematically varied island area (25, 100, 400 and 1,600 m2) within a 28 m radius were included, which corresponds to the diameter
and planted diversity (zero, one, two, three and six tree species). of the largest study plot (40 × 40 m2). We used the maximum number
The six planted tree species (Archidendron jiringa ( Jack) I.C.Nielsen of individuals detected simultaneously in all recordings per plot as a
(Fabaceae), Parkia speciosa Hassk (Fabaceae), Durio zibethinus conservative proxy of abundance per bird species or bat sonotype.
L. (Malvaceae), Dyera polyphylla (Miq.) Steenis (Apocynaceae), Shorea
leprosula Miq. (Dipterocarpaceae) and Peronema canescens Jack (Lami- Understorey arthropods. Arthropods were sampled in the under-
aceae)) are native to the region and widely used for their fruits, timber storey vegetation during October 2016 to January 2017. Each plot was
or latex35. Around 40% of the oil palms located inside the tree islands sampled three times with six pan traps per plot exposed for 45 h. Traps
were felled, with the number of felled oil palms differing depending on were made of white plastic soup bowls covered with yellow ultraviolet
the tree island area33. The trees were planted between the felled and spray-paint53 and were filled with water and one drop of regular soap.
standing oil palms on a 2 m triangular grid. The tree islands were fenced They were fixed in a holding system in groups of three at the height of
and the management comprised a total stop of fertilizer, herbicide and the surrounding plants and these systems then equally distributed in
pesticide application after planting. After May 2016, manual weed- distance from edge and to each other. All arthropods were preserved
ing was restricted to 1 m circles around the planted trees when these in 70% ethanol. Subsequently, all individuals were identified to higher
were shorter than the surrounding grass layer, allowing for natural taxonomic groups and morphospecies. The taxonomic groups Hyme-
regeneration. In addition to the 52 tree islands, we established four noptera, Lepidoptera and Araneae were categorized into functional
control plots in the oil palm plantation that have a fixed area (100 m2) groups (pollinators, predators and parasitoids) using different identi-
and that were managed conventionally (that is, no oil palm was felled, fication keys54–60 (Supplementary Table 10). Predators and parasitoids
no tree was planted and application of fertilizer, herbicide and pesticide were merged into the single functional group ‘predators’.
was as usual), in the main text referred to as conventionally oil palm
monocultures. In total, the experiment comprises 56 study plots19. In Soil fauna. During October–November 2016, in each plot, four soil
each study plot larger than 25 m2, one subplot of 5 × 5 m2 was estab- samples of 16 × 16 cm2 were taken randomly within the subplot with
lished in a random location at a minimum distance of 1.5 m from the a spade. Samples included litter (if present) and soil down to a depth
plot edge. of 5 cm. Animals were extracted using a gradient heat extractor61 and
collected in dimethyleneglycol–water solution (1:1) and thereafter
Field measurements transferred into 70% ethanol62. All extracted animals were counted
We conducted an interdisciplinary field campaign from October 2016 to and sorted into 28 taxonomic groups (in most cases orders) allow-
October 2018, that is, 33–57 months after establishment of the experi- ing for functional group classifications63; Extended Data Table 4. We
ment. At this early stage of the experiment, the tree islands already calculated community metabolism of all animals that were classified
differed in their structural complexity30 and the planted trees reached as detritivores, herbivores and predators in a sample by using mean
up to 16 m height35. In all the 56 study plots, several indicators related to group- and ecosystem-specific estimates derived from ref. 63. The
biodiversity, ecosystem functioning and structure were measured using estimates are based on measurements of more than 5,000 individuals
standardized procedures and constant sampling areas at the level of the of soil animals across eight different oil palm plantations in the same
plot or subplot (Supplementary Tables 1–3). Only trees were sampled at region; to estimate community metabolism, individual body masses
unequal areas (that is, all trees present in the plots were sampled) and were recalculated to metabolic rates using group-specific regressions
were therefore standardized using rarefaction curves (see ‘Trees’). Oil from ref. 64. Community metabolism was calculated by summing up
palm yield was continuously measured since the beginning of the experi- metabolic rates of all individuals; we used the mean per plot across
ment at the level of individual palm but the data were then aggregated four samples for each functional group (detritivores, herbivores and
over space and time (see ‘Per area yield’ and ‘Per island yield’). Each vari- predators) for the analysis. We also computed taxonomic diversity
able presented in the main text had one measurement per plot, such that as the number of taxonomic groups (in most cases orders) present in
blinding and randomization were not applicable. No statistical methods each plot for the analysis.
were used to predetermine sample size. The data were processed and
analysed in R v.1.2.1335 (ref. 48). Fungi. In January 2017 three soil cores (10 cm depth, 4 cm diameter)
were taken within each 5 × 5 m2 subplot. Surface leaf litter was removed
Birds and bats. We recorded audible and ultrasound in March 2017 before soil collection. The soil was sieved through a 50 × 50 mm2
using automated sound recorders (SM2Bat+ recorders, Wildlife acous- sieve and roots were separated from soil. The fungal community was
tics; with an acoustic SMX-II microphone on the left channel and one assessed using Illumina next-generation sequencing (Illumina) of the
full-spectrum Sonitor Parus49 microphone on the right channel), ITS2 marker region. The detailed protocol for amplification, amplicon
strapped to wooden poles at a height of 1.5 m in the centre of the plot. sequencing and generation of fungal operational taxonomic units
On consecutive days, we extracted sound recordings for sampling (OTU) is described in ref. 65. OTUs were classified taxonomically
birds and insectivorous bats. We used two stereo 15 min recordings using the BLAST (blastn, v.2.7.1) algorithm66 and the UNITE v.7.2 (UNITE_
starting 15 min before sunset and two 15 min stereo recordings starting public_01.12.2017.fasta) reference database67.
Article
Prokaryotes. In May 2017 three cores of topsoil (10 cm depth) were Pollen. To collect pollen/spore rain, Behling pollen traps83 were
taken in each subplot. Soil cores were then mixed, homogenized and installed from June until October 2018. Each trap consists of a plastic
freed from roots. A total of 5 ml of RNAprotect Bacteria reagent (Qiagen) tube which is placed about 30 cm above the ground and is held by a fix-
was added to 5 g of soils to prevent nucleic acid degradation. DNA and ing pole. The tube is filled with 5 ml of liquid glycerol, synthetic cotton
RNA were extracted from 1 g of soil by using the Qiagen RNeasy Power- and, on the top, it is covered by a mosquito net to reduce disturbance
Soil Total RNA Kit and the RNeasy PowerSoil DNA Elution Kit (Qiagen). from animals or litter and prevent the cotton from being removed. In
The V3–V4 region of the 16S ribosomal RNA gene was amplified and tropical regions heavy rainfalls occur, thus it is necessary to prevent the
sequenced as described in ref. 68. Paired-end sequences were quality pollen from pouring out of the pollen trap. In the Behling trap, glycerol
filtered with fastp (v.0.20.0)69 and merged with PEAR v.0.9.11 (ref. 70). is used, which has a higher density compared to water. Consequently,
Remaining primer sequences were clipped with cutadapt v.2.5 (ref. 71). the incoming rainfall can flow out of the trap without taking the pol-
Size filtering, dereplication, denoising and chimaera removal was per- len, which is trapped in the synthetic cotton and in the glycerol83. The
formed with vsearch v.2.12.0 (ref. 72). Curated sequences were then clas- Behling traps were modified to mimic the surrounding environment
sified by mapping each sequence against the SILVA database with the and maximize recovery. In total, 168 pollen traps were installed in the
BLAST73. Counts were normalized by using the GMPR normalization74. plots (3× plot). Of the total 56 plots, the pollen traps were not recovered
in three plots (P28, P34 and P47). One pollen trap from each 53 plot was
Seeds. We installed four seed traps in each of the 56 study plots for processed and analysed. Firstly, each pollen trap was washed with dis-
1 yr, that is, between 1 April 2017 and 29 March 2018. The traps were tilled water through a 2 mm mesh sieve to remove large size materials.
built using fine-mesh cloth attached to a squared structure made in Afterwards, the pollen traps were sieved through a 150 µm mesh sieve
PVC pipes of size 50 cm × 50 cm fixed at 1 m from the ground. The traps to exclude medium-sized materials from the samples. Two Lycopodium
were installed at random locations in each of the four quadrants of tablets were added as markers to each sample to estimate palynomorph
each plot, at a minimum distance of 1 m from the plot edge. The con- concentrations84 and the Erdtman acetolysis85 was applied, to remove
tents of the traps were collected twice a month, dried at approximately cellulose material. Residues were mounted in glycerol jelly for pollen
40 °C during 3–7 days. All the seeds were carefully extracted from the visualization, identification and counting. Pollen and spore analyses
samples, counted and separated by morphospecies using hand lenses were carried out using light microscopy. All identified pollen and spore
(×10 magnification) and a microscope (Leica photomicroscope with types were photographed using a Leica photomicroscope with ×400
×400 magnification) for very small seeds. Molecular identification magnification. For each trap, a total sum of at least 100 pollen grains
of the morphospecies was implemented using three universal plant were counted. Pollen and spore grains can rarely be identified to spe-
DNA barcodes (matK, rbcL and ITS2)75–78 and taxonomic assignments cies level and the level of taxonomic identification varies for different
were made using BLASTn search against the NCBI Genbank reference groups of plants. Consequently, a reduction to the family level has
sequence database79. Sequences obtained from the barcode loci were been proposed for studies involving analysis of palynological diversity
deposited in NCBI Genbank under the accession numbers OM811991– in the tropics86.
OM812021, OM837673–OM837724 and OM935782–OM935815. We clas-
sified each morphotype as native or non-native species using available Pollination. We assessed pollination rate on chilli pepper plants
literature80,81 (http://www.plantsoftheworldonline.org/). We derived (Capsicum annuum) as phytometer plants, selected for potential shade
the native seed density (number of native seeds per m2) as the total num- tolerance87, widespread home garden cultivation in this region88 and
ber of native seeds over the entire sampling duration per plot, which the potential role pollination can play in fruit quality and yield89. We
was used as an indicator of ecosystem functioning (see ‘Ecosystem raised 1,500 individuals of a locally available variety of C. annuum from
functioning’). Seed diversity, calculated on the basis of the Hill number seed outside of the study plots. During the growth period outside the
frameworks and used as indicators of biodiversity (see ‘Biodiversity’), study plots, we applied NPK fertilizer and pesticide (imidacloprid,
was derived from the pooled samples per plot over the entire sampling deltamethrin, mancozeb and abamectin) following local practices
duration for all seeds (native and non-native). to standardize growing conditions and control pest damage before
transfer to field sites. We halted fertilizer and pesticide application
Herbs. All non-woody terrestrial vascular plants (for example, angio­ 1 week before placement in the plots and only watered as conditions
sperm herbs and vines, ferns, but not epiphytes) in the subplot were required thereafter. In February 2018, we selected 224 healthy individu-
recorded from February until March 2018. They were classified to species als of comparable size to transfer to the 56 study plots (four plants per
or morphospecies and herbaceous cover (in absolute per cent ratios plot). The four chilli plants were placed, still in their pots, at the centre
from 1% to 100%) was estimated by two people. Epiphytes growing of each plot for 5 weeks for a period of open pollination and monitor-
on the stems of trees or palms were excluded, whereas vine species ing, followed by 3 weeks for fruit harvesting. We removed any flowers
that rooted in the ground and climbed up stems of trees or palms were before placement in the field, so pollinated flowers and developing
included. Herbarium specimens were collected and stored in the labo- fruits were assumed to result from pollination within the study plots.
ratory of Jambi University. All names were checked following The Plant During the period before final harvest, each plot was revisited once
List 2013, v.1.1 (http://www.theplantlist.org). per week and the number of flowers were counted per plant. The rate
of successful pollination was estimated from the fruit to flower ratio,
Trees. All planted trees were surveyed in January to February 2018 as which was the total number of harvested fruits divided by the total
part of a yearly inventory35. Furthermore, we surveyed all free-standing number of flowers observed per plot.
woody plants (trees, shrubs and bamboos) that colonized the plots
with a length of ≥130 cm from April until August 2018. For each species Per palm yield. We followed the conventional harvesting procedure
or morphospecies, one voucher specimen was collected, dried and established by the plantation manager of PT Humusindo and measured
pressed according to standard procedure. In the main text, we refer to the weight of the fresh fruit bunches directly after harvest using a port-
the colonized woody plants as ‘trees’, unless stated otherwise. Because able scale. We measured the per palm yield (kg per palm) of all palms
the number of sampled trees largely varied according to the tree island inside the 52 tree islands (N = 214) and one palm per control plot with
area, we standardized the diversity estimates using rarefaction curves conventional management (N = 4). To obtain a more solid estimate for
(R package iNEXT)82 to 24 individuals, which represent the median the conventional plantation, we measured the per palm yield of 30 more
number of individuals per plot. reference palms that were evenly distributed across the conventional
plantation at approximately an equal distance to each tree island and N
 ad, number of oil palms directly adjacent to the island (that is, adja-
whose neighbourhood is characteristic of conventionally managed cent position 1)
oil palm monocultures (Supplementary Note 2 and Supplementary Nfelled, number of removed oil palms in the island
Figs. 1–3). To examine potential changes in yield in the conventionally Yp_ref, median per palm yield of the reference palms in the convention-
managed oil palm plantation surrounding the tree islands (‘spillover ally managed oil palm plantation (kg palm−1)
effects’), we measured the per palm yield of three oil palms adjacent to Yp_in median per palm yield inside the tree island (kg palm−1)
each tree island, at increasing distance to the island’s edge (at position Yp_adj, per palm yield directly adjacent to the tree island (that is, adja-
number 1, 2 and 3)33. For direct comparison with earlier findings, we cent position 1 (kg palm−1)).
analysed the data following established methodology33. The tests were
based on linear mixed-effect models with the annual yield of individual Above-ground biomass. For all the planted trees, we measured the
oil palms as the response variable and the plot identity as the random basal diameter (at 10 cm above ground), the diameter at breast height
effect. Pairwise comparison was conducted with a post hoc Tukey test. (130 cm above ground) and the tree height in 2017 as part of a yearly
Because our results indicate that only the palm directly adjacent to the inventory35. In January and February 2017, we also measured the height
tree island (in position number 1) was affected by the experimental of the oil palms at meristem, that is, the point of attachment of the
treatment (Extended Data Fig. 5a, in agreement with ref. 33), we do not young leaves to the oil palm trunk30. We estimated above-ground bio-
consider the palms in position 2 and 3 in the yield calculation per island mass of the trees (equation (5)) and the oil palms (equation (6)) using
(Per island yield). The per palm yields in and adjacent to the 56 study the respective allometric equations of refs. 90,91:
plots have been continuously monitored since the establishment of
the experiment in December 2013. The extra 30 palms were established AGBtree = 0.0673 × (ρ × DBH2 × H )0.976 (5)
in December 2016. For consistency with other indicators (Extended
Data Tables 1–3), we reported yield data for 1 yr (November 2017 until AGBpalm = 71.797 × H − 7.0872 (6)
October 2018) in the main text and for 2 yr (November 2016 until
October 2018) in the Supplementary Note 1. Yield data since 2014 are AGBtree, above-ground biomass of the planted trees (kg tree−1)
shown in Extended Data Fig. 4, in which the oil palms in position 3 were AGBpalm, above-ground biomass of the oil palms (kg palm−1)
used as reference palms for the corresponding time period. DBH, tree diameter at breast height (cm)
H, height of tree or palm (m)
Per area yield. We estimated per area yield (∆Yha, kg ha−1) as the yield ρ, wood density (g cm−3).
of a given palm (kg palm−1) multiplied by a stand density-dependent Wood density for Peronema canescens (0.61 g cm−3), Parkia speciosa
expansion factor (EF) to derive estimates of per area yield (kg ha−1). (0.54 g cm−3) and Dyera polyphylla (0.36 g cm−3) was based on EFForTS
We then calculated the per area yield change between tree islands and core plot data, whereas for Archidendron sp. (0.36 g cm−3), Shorea
reference (kg ha−1; Supplementary Note 3). This approach accounts leprosula (0.44 g cm−3) and Durio zibethinus (0.516 g cm−3) it was taken
for changes in per area yield due to oil palm thinning (that is, reduced from the global wood density database92.
oil palm densities and changes in per palm yield in the tree islands) We estimated the total above-ground biomass per plot as the sum of the
but does not account for potential changes in per palm yield on the above-ground biomass of the palms and the planted trees (equation (7)).
surrounding plantation, for example, because of spillover effects19. The estimations of total AGB did not consider the necromass, litter,
An alternative analysis considering spillover effects was performed understorey vegetation and spontaneously established trees, which
at the plot level (Per island yield). were considered negligible.

Per island yield. We estimated oil palm yield changes at the tree island AGB = ( ∑ AGBtree + ∑ AGBpalm /Nin × d palm )/(A × 10, 000) (7)
scale (∆YIsland, in kg island−1); equations (1)–(4) following established
methodology33. This method considers the yield foregone owing to the AGB, total above-ground biomass per plot (t ha−1)
removal of some oil palms before the experiment, as well as changes in dpalm, density of oil palms (number of oil palms per ha) that takes
per palm yield inside the tree islands and directly adjacent to the tree into account the local neighbourhoods of the plots (also referred to
islands (at position 1, that is, spillover effects). Because the number as EF; see Supplementary Note 3)
of oil palms inside and adjacent to the tree islands and the number of A, area of the plot (m2).
removed oil palms vary depending on tree island area33, the net oil palm
yield changes are provided per plot and not per area. Even though this Tree growth. The growth of the planted trees per plot was calculated
method was initially designed to calculate the oil palm yield changes as35:
for the tree islands, here we also apply it to the four control plots to
integrate them in our synthesis analysis. BA inc,2017−2018 = ∑ (BA tree,2018 − ∑ BA tree,2017)/A. (8)

∆YIsland = YSpillover + YRemainChange − YForegone (1) BAinc, 2017–2018, total plot-level basal area increment between 2017 and
2018 (in cm2 m−2 yr−1, equivalent to m2 ha−1 yr−1)
YForegone = Nfelled × Yp_ref (2) BAtree, year, tree basal area (in cm2) derived from the basal diameter
(cm) in the specific year.
YRemainChange = Nin × (Yp_in − Yp_ref) (3)
Leaf litter input. We measured leaf litter fall (in g m−2 yr−1) using the
four seed traps installed randomly in each four quadrants of the plots
YSpillover = Nadj × (Yp_ adj − Yp_ ref) (4) from April 2017 to March 2018 (Seeds). The contents of the traps were
collected twice a month, dried at 40 °C for 4–7 days and weighted. We
∆YIsland, per island oil palm yield change (kg island−1) also sorted the leaves by species and weighted the content for the six
YSpillover, per island yield changes due to spillover effects (kg island−1) planted tree species and oil palm separately.
YRemainChange, per island yield changes inside the island (kg island−1) For each sampling date, we aggregated the values at plot level using
YForegone, per island yield foregone due to oil palm removal (kg island−1) the median per plot of the litter weight. We then excluded outliers
Ni, number of remaining oil palms inside the island defined as plot-level values outside the range of 3 standard deviations
Article
around the median of the entire data (less than 5% of the litter weight were inserted horizontally into the first 5 cm of topsoil. The soil was
data, total and per species). To get annual estimates, we summed the weighed, dried at 105 °C until constant weight and weighted again.
available plot-level values over time and divided them by the number of Calculation was done on the dry weight basis, for which the sample
sampling dates (that is, between 17 and 24, depending on the number dry weight (g) was divided by the volume of the sample (cm3) collected
of missing traps or excluded outliers). We then multiplied the obtained from the average of the five replicates. We used the mean per plot for all
values by the seed trap area (0.25 m2) to get the leaf litter fall in g m2 yr. mentioned soil variables and used the inverse of C:N ratio and soil bulk
density as measures of soil fertility and soil decompaction, respectively
Leaf litter decomposition. We installed litterbags (20 × 20 cm2, (Supplementary Table 1).
4 mm mesh size) each filled with 12 g of material: 6 g of freshly cut and
air-dried (approximately 25 °C) fronds of oil palm leaves93 and 6 g of the Vegetation structure. We measured 12 variables representing various
freshly fallen air-dried leaf litter for each tree species or their combi- aspects of the vegetation structure (Supplementary Table 3). We used
nations in experimental plot. In each plot, one litterbag was installed a terrestrial laser scanner Focus M70 (Faro Technologies) to create
in November 2017 for a duration of 6 months. Decomposition (litter three-dimensional point clouds of the vegetation at the centre of each
mass loss) was calculated as the difference between the initial litter plot in September and October 2016, as described in ref. 30. We comput-
dry mass and litter dry mass remaining after 6 months and expressed ed the (1) stand structural complexity index (SSCI) following ref. 99 and
as a percentage of decomposed material. its two components: (2) the mean fractal dimension index (MeanFRAC)
derived from cross-sections of polygons in the three-dimensional point
Water infiltration capacity. To quantify soil water infiltration capacity, cloud, which is a scale-independent and density-dependent measure
we measured saturated soil hydraulic conductivity (Kfs, cm h−1) using a of structural complexity and (3) the effective number of layers (ENL)
dual-head infiltrometer (Saturo) in March 2018 near the subplot centre that describes vertical stratification based on the Simpson Index100. ENL
in 35 (out of the 52) tree islands and in the four control plots represent- and MeanFRAC are integrated in the SSCI and all these three measures
ing conventionally managed oil palm monocultures. Owing to a broken were derived from vegetation parts above 130 cm. We also derived (4)
instrument, the 17 remaining plots were measured using a custom the understory complexity index that measures the fractal dimen-
manual double-ring infiltrometer, which tends to yield higher Kfs esti- sion of horizontal cross-sections of the point cloud between 80 and
mates than the dual-head approach because there is no correction for 180 cm height, thereby measuring the structural complexity of the
lateral flow. In three plots, Kfs was measured with both devices. We plot- understorey vegetation101. (5) Canopy gap fraction was estimated
ted these values against each other and found a close linear relationship from hemispherical photographs at plot level as described in ref. 30.
(R² = 0.98, P = 0.066); even though it was only marginally significant Drone-based photogrammetry dated from September to October 2016
because of the small sample size, we used it to correct the values from was used to further partition the canopy (in %) as (6) oil palm cover and
the 17 plots that were measured manually (Kfs_corr = 1.44 + 0.55 Kfs_double_ring) (7) tree cover as described in ref. 102. We also used the drone-based
to allow for comparability across all 56 plots. orthophotos to calculate (8) oil palm density as the number of living
oil palms per plot irrespective of the orientation of the plot relative to
Evapotranspiration. We recorded land and canopy surface tempera- the planting scheme (Supplementary Figs. 2 and 3) . For the smaller
tures using a radiometric thermal camera (FLIR Tau 2 640, FLIR Sys- plots (5 × 5 m) unaffected by thinning, the oil palm density was simply
tems) attached to a TeAx ThermalCapture module (TeAx Technology the typical planting density in conventionally managed oil palm planta-
GmbH) mounted on a multicopter drone (MK EASY Okto V3; HiSys- tions (120 planted palms per hectare). Further details on the oil palm
temsy) as described in ref. 94. Image sets were recorded four to five density calculation are given in the Supplementary Note 3. We also
times per day around noon, covering each plot once over a 9 day period calculated (9) tree density as the number of trees planted and from
encompassing varying weather conditions. Land and canopy surface natural regeneration per plot and expressed per hectare. We estimated
temperatures were the main input for modelling latent heat flux (in the portion of the ground (in percent) as (10) understorey vegetation
W m−2) and deriving evapotranspiration using the QWaterModel QGIS3 cover and (11) litter cover per subplot in February–March 2018. The
Plugin95, which is based on the DATTUTDUT energy balance model96. understorey vegetation cover included all parts of plants lower than
Measured short-wave radiation and relative humidity were used as 130 cm in height, including the trunks and other parts of the planted
further input variables to support the prediction of latent heat flux trees but excluding oil palm trunks. (12) The litter depth (cm) was meas-
and derive evapotranspiration. ured as the mean value in three randomly chosen positions inside each
subplot with a metal ruler. To extract orthogonal axes (PC1 and PC2)
Microclimate. We measured microclimate using temperature per that represent most of the variability in the vegetation structure, we
humidity loggers: hydrochron (DS1923-F5) and thermochron applied a principal component analysis on all the structural variables
(DS1922L-F5) iButtons, Maxim integrated. The loggers were installed after standardization to zero mean and unit variance.
in the middle of each plot at 1.5 m above ground and were protected
from water and direct solar radiation using handmade multiplate radia- Restoration outcomes
tion shields97. Data were collected for 1 yr (18 November 2017 until 19 Ecosystem functioning. We measured 20 variables related to seven
September 2018) every 3 h starting at midnight. As a proxy for micro- categories of ecosystem functioning including: productivity as (0) tree
climate buffering, we calculated the daily amplitude as the absolute growth (basal area increment of the planted trees in m2 ha−1 yr−1) that was
difference between values at 7:00 and 15:00 (ref. 97), aggregated using further excluded from the analysis—see Supplementary Fig. 4, (1) oil
the median value over the entire measurement period. palm yield (per island oil palm yield changes in kg of fresh fruit bunches
per island) and (2) above-ground biomass (biomass stored in the aerial
Soil properties. We determined soil total carbon content (g mg−1), parts of the planted trees and the oil palms, in t ha−1); resistance to inva-
total nitrogen content (mg g−1) and plant available phosphorus content sion as (3) native seeds (total number of arriving native seeds per m2) and
(mg g−1) using the same three soil samples as for fungi community data (4) resistance to invasive plants (100—observed cover of Clidemia hirta,
(see ‘Fungi’) and the method of determination is described in detail in %); pollination as (5) pollinators (number of sampled individuals)
in ref. 65. We then calculated the C:N ratio accounting for the molar and (6) pollination rate (fraction of flowers on phytometer plants that
mass of the elements following ref. 98, that is, 12.0107 for carbon and are pollinated, %); soil quality as (7) soil P (phosphorous content, %),
14.0067 for nitrogen. We also measured soil bulk density (g cm−3) using (8) 1/soil C:N (the molar ratio of soil C to soil N concentration) and (9)
five soil samples taken in the subplot in May 2018. Soil rings of 100 cm3 soil decompaction (inverse of soil bulk density in g cm−3); predation
and herbivory as (10) predatory invertebrates (total activity duration monocultures), island area (plot edge length in m), planted diversity
of insectivorous bats and birds, in seconds); (11) predatory arthropods and the restoration outcome (either biodiversity or ecosystem func-
(number of sampled individuals), (12) predatory soil fauna (energy tioning indicators) as single factors and tree island × indicator, island
flux, in J h−1), (13) herbivory (energy flux, in J h−1); carbon and nutrient area × indicator, planted diversity × indicator, island area × planted
cycling as (14) decomposers (energy flux, in J h−1); (15) litter decomposi- diversity and island area × planted diversity × indicator interactions.
tion (relative biomass loss of litter after 6 months in litterbags, %) and For conventionally managed oil palm plots, island area was set to 10 m
(16) litter input (biomass of leaf litter falling in traps, g m−2); water and edge length and planted diversity to zero. Each response variable
climate regulation as (17) evapotranspiration (canopy latent heat flux, (biodiversity and ecosystem functioning indicators) was standardized
in W m−2); (18) soil water infiltration capacity (saturated soil hydrau- to unit scale (between 0 and 1) as this improved the model diagnostics
lic conductivity in cm h−1) and (19) microclimate buffering (median before applying the respective linear mixed-effect models; whereas we
daily amplitude of air temperature during 1 yr, °C d−1). A more detailed used logarithmic transformations for island area and planted diversity.
summary of the 20 ecosystem functioning variables is presented in Plot was included as a random term.
Supplementary Table 2. As an alternative to the linear mixed-effect models, we applied
Kruskal–Wallis tests on each indicator of biodiversity and ecosys-
Biodiversity. We derived taxonomic diversity for soil bacteria and soil tem functioning for comparison between the 52 tree islands and the
fungi, soil fauna, herbs, trees, seeds, pollen, understorey arthropods, four conventionally managed oil palm monocultures as control plots
birds and bats. Most of the groups (arthropods, herbs, trees, birds and (Supplementary Note 5 and Supplementary Tables 8 and 9).
seeds) were sorted at the lowest possible taxonomic level (species or
morphospecies). Pollen, soil fauna and bats were sorted to higher levels, Structural equation modelling. We used piecewiseSEM109 to assess the
mainly family, order and morphotypes, respectively. Soil bacteria and influence of tree island area and tree planted diversity on biodiversity
soil fungi were analysed by DNA-based marker gene sequencing as and ecosystem functioning operating through increasing tree domi-
amplicons sequence variants or OTU, respectively. Hereafter, we refer nance, through differences in structural complexity (indirect effects)
to these different taxonomic units (species, family, order, morphotypes or through alternative mechanisms (direct effects). As a hypothetical
and OTU) as ‘species’ for simplicity. causal model, we included direct paths between island area and tree
Diversity was measured following the Hill number framework, which dominance (PC2) and between tree planted diversity and structural
allows comparison across diversity indices that weigh relative abun- complexity gradient (from open to dense and structurally complex
dances to varying extents (species richness, Shannon diversity and vegetation, PC1; Extended Data Fig. 2). Piecewise SEMs are based on
Simpson diversity) and are expressed in terms of effective numbers a set of linear equations which are evaluated individually109. For our
of species103–106. Species richness is more sensitive to locally rare spe- analyses, we included:
cies, Simpson diversity is more sensitive to locally dominant species lm(restoration outcome ≈ island area + planted diversity) (1)
and Shannon diversity favours neither rare nor dominant species. We lm(structural complexity ≈ planted diversity) (2)
show results for species richness and Simpson diversity in the main lm(tree dominance ≈ island area) (3)
text and for all indicators in the Extended Data Tables 1 and 2 and Across all restoration outcomes, the main variables were always
Supplementary Tables 4– 8. The calculations were performed using included in the linear model (1). As tree dominance and structural com-
the R packages iNext82 and vegan107. plexity are potential mechanistic pathways explaining the influence of
island area and tree planted diversity, alternative paths between them
Multidiversity and multifunctionality. Different indicators of bio- and biodiversity or ecosystem functioning were added, if they improved
diversity and ecosystem functioning were aggregated by calculating the model fit (based on modification indices, P < 0.05). Therefore,
multidiversity and multifunctionality, respectively. Following ref. 18, model selection influenced only the inclusion of structural complexity
we performed a cluster analysis to preselect indicators for achieving and tree dominance in the linear model (1). Effects of island area and
a representative measure of ‘ecosystem function multifunctionality’. planted diversity through mechanistic pathways were calculated by
As tree growth and litter input were correlated and formed a cluster, multiplying their effect on the mechanistic explanatory variable and
we excluded tree growth from the analysis (Supplementary Note 4 and the effect of the mechanistic explanatory variable on biodiversity or
Supplementary Fig. 4). Following a threshold approach108, we calculated ecosystem functioning. Mechanisms that were not captured by either
multifunctionality (and multidiversity) as the number of ecosystem of our proposed mechanistic pathways are represented by the direct
functioning (and biodiversity) indicators that cross a threshold, ex- paths between island area and tree planted diversity and biodiversity or
pressed as a certain percentage of the maximum observed values in ecosystem functioning. We tested the assumption of normality of the
our study landscape (among all 56 study plots). We calculated multi- residuals in models (1), (2) and (3) using Shapiro–Wilk normality test,
functionality and multidiversity for all thresholds from 1% to 99% and applied a suitable transformation of the response variables if needed
presented results for a 50% threshold in the main text. To reduce the (package bestNormalize v.1.6.1). The transformation concerned four
influence of extreme values, we used the mean of the three highest out of ten indicators for species richness and Shannon diversity, three
values observed in all study plots, respectively. As an alternative to the indicators out of ten indicators for Simpson diversity and 14 indicators
threshold approach, we also calculated multidiversity and multifunc- out of 19 for ecosystem functioning. Effect sizes were calculated using
tionality as the average of the indicators108. Before multidiversity and standardized coefficients. The island area and planted tree diversity
multifunctionality calculations, all the variables were standardized were log-transformed as this improved the model fit. For each SEM,
to unit scale (for biodiversity and ecosystem functioning separately). we quantified the goodness of fit using the following metrics: Fisher’s
The calculations were performed using the package multifunc in R108. C statistic and significance value based on a Chi-square test, the infor-
mation criterion (Akaike information criterion (AIC), Bayesian informa-
Statistical analysis tion criterion (BIC), corrected AIC (AICc)) and pseudo-R2 values and
Linear mixed-effect models. We used linear mixed-effect models applied the test of directed separation as implemented in the package
to test the effects of the experimental treatment on restoration out- piecewiseSEM v.2.1.0109.
comes. We fitted three separated models for biodiversity using species
richness, Shannon and Simpson diversity as response variables and Inclusion and ethics
one model for ecosystem functioning. These models included tree The research included researchers from the Indonesian institutes Jambi
island (compared to our controls of conventionally managed oil palm University and Bogor Agricultural University throughout the research
Article
process—study design, study implementation, data ownership, intel- 66. Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinformatics 10, 421
(2009).
lectual property and authorship of publications. Local and regional 67. Kõljalg, U. et al. Towards a unified paradigm for sequence-based identification of fungi.
research relevant to our study was considered in citations. Mol. Ecol. 22, 5271–5277 (2013).
68. Berkelmann, D., Schneider, D., Hennings, N., Meryandini, A. & Daniel, R. Soil bacterial
community structures in relation to different oil palm management practices. Sci. Data 7,
Reporting summary
421 (2020).
Further information on research design is available in the Nature 69. Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor.
Portfolio Reporting Summary linked to this article. Bioinformatics 34, i884–i890 (2018).
70. Zhang, J., Kobert, K., Flouri, T. & Stamatakis, A. PEAR: a fast and accurate Illumina
paired-end reAd mergeR. Bioinformatics 30, 614–620 (2014).
71. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing
Data availability reads. EMBnet.journal 17, 10–12 (2011).
72. Rognes, T., Flouri, T., Nichols, B., Quince, C. & Mahé, F. VSEARCH: a versatile open source
The raw data are available at https://data.goettingen-research-online. tool for metagenomics. PeerJ 4, e2584 (2016).
de/dataverse/crc990, with the specific link for each dataset provided in 73. Quast, C. et al. The SILVA ribosomal RNA gene database project: improved data
Supplementary Tables 1–3. The processed data are available at https:// processing and web-based tools. Nucleic Acids Res. 41, D590–D596 (2013).
74. Chen, L. et al. GMPR: a robust normalization method for zero-inflated count data with
doi.org/10.6084/m9.figshare.22320490. Seed DNA sequences are application to microbiome sequencing data. PeerJ 6, e4600 (2018).
available in NCBI Genbank under the accession numbers OM811991– 75. White, T. J., Bruns, T., Lee, S. & Taylor, J. W. in PCR Protocols: A Guide to Methods and
OM812021, OM837673–OM837724 and OM935782–OM935815. Applications (eds Innis, M. A. et al.) 315–322 (Academic Press, 1990).
76. Taberlet, P., Gielly, L., Pautou, G. & Bouvet, J. Universal primers for amplification of three
Sequencing data of the soil fungal community were deposited in the non-coding regions of chloroplast DNA. Plant Mol. Biol. 17, 1105–1109 (1991).
NCBI Sequence Read Archive (SRA) under Bioproject accession number 77. Kress, W. J. & Erickson, D. L. A two-locus global DNA barcode for land plants: the coding
PRJNA659225. The public UNITE database (https://unite.ut.ee/) v.7.2 on rbcL gene complements the non-coding trnH-psbA spacer region. PLoS ONE 2, e508
(2007).
fungal ITS sequences was used as a reference of taxonomic classifica- 78. CBOL Plant Working Group et al. A DNA barcode for land plants. Proc. Natl Acad. Sci. USA
tion. Sequence data of the bacterial communities were deposited in the 106, 12794–12797 (2009).
NCBI SRA under Bioproject accession number PRJNA841353. Sequence 79. Zhang, Z., Schwartz, S., Wagner, L. & Miller, W. A greedy algorithm for aligning DNA
sequences. J. Comput. Biol. 7, 203–214 (2000).
identification was performed by mapping all curated sequences against 80. Rembold, K., Mangopo, H., Tjitrosoedirdjo, S. S. & Kreft, H. Plant diversity, forest
the SILVA database v.132 (https://www.arb-silva.de/). dependency, and alien plant invasions in tropical agricultural landscapes. Biol. Conserv.
213, 234–242 (2017).
81. van Kleunen, M. et al. The Global Naturalized Alien Flora (GloNAF) database. Ecology 100,
e02542 (2019).
Code availability 82. Hsieh, T. C., Ma, K. H. & Chao, A. iNEXT: an R package for rarefaction and extrapolation of
species diversity (Hill numbers). Methods Ecol. Evol. 7, 1451–1456 (2016).
The R code used in the current study is available at https://doi.org/
83. Jantz, N., Homeier, J., León-Yánez, S., Moscoso, A. & Behling, H. Trapping pollen in the
10.6084/m9.figshare.22320490. tropics—comparing modern pollen rain spectra of different pollen traps and surface
samples across Andean vegetation zones. Rev. Palaeobot. Palynol. 193, 57–69 (2013).
84. Stockmarr, J. Tablets with spores used in absolute pollen analysis. Pollen Spores 13,
45. Paquette, A. et al. A million and more trees for science. Nat. Ecol. Evol. 2, 763–766 615–621 (1971).
(2018). 85. Erdtman, G. Handbook of Palynolgy: Morphology, Taxonomy, Ecology. An Introduction to
46. Drescher, J. et al. Ecological and socio-economic functions across tropical land use the Study of Pollen Grains and Spores (København Munksgaard, 1969).
systems after rainforest conversion. Phil. Trans. R. Soc. B 371, 20150275 (2016). 86. Jantz, N., Homeier, J. & Behling, H. Representativeness of tree diversity in the modern
47. Bell, T. et al. A linear model method for biodiversity—ecosystem functioning experiments. pollen rain of Andean montane forests. J. Veg. Sci. 25, 481–490 (2014).
Am. Nat. 174, 836–849 (2009). 87. Pouliot, M., Bayala, J. & Ræbild, A. Testing the shade tolerance of selected crops under
48. R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Parkia biglobosa (Jacq.) Benth. in an agroforestry parkland in Burkina Faso, West Africa.
Statistical Computing, 2018). Agrofor. Syst. 85, 477–488 (2012).
49. Darras, K. F. A. et al. Assembling cheap, high-performance microphones for recording 88. Prabowo, W. E. et al. Bird responses to lowland rainforest conversion in Sumatran
terrestrial wildlife: the Sonitor system. F1000Research 7, 1984 (2021). smallholder landscapes, Indonesia. PLoS ONE 11, e0154876 (2016).
50. Darras K. F. A., Pérez N., M., Dilong L., Hanf-Dressler T., Markolf M., Wanger T. C. ecoSound- 89. Roldán Serrano, A. & Guerra-Sanz, J. M. Quality fruit improvement in sweet pepper
web: an open-source, online platform for ecoacoustics [version 2; peer review: 2 approved]. culture by bumblebee pollination. Sci. Hort. 110, 160–166 (2006).
F1000Research, 9:1224 (2023) (https://doi.org/10.12688/f1000research.26369.2) 90. Chave, J. et al. Improved allometric models to estimate the aboveground biomass of
51. Darras, K., Furnas, B., Fitriawan, I., Mulyani, Y. & Tscharntke, T. Estimating bird detection tropical trees. Glob. Change Biol. 20, 3177–3190 (2014).
distances in sound recordings for standardizing detection ranges and distance sampling. 91. Asari, N., Suratman, M. N., Jaafar, J. & Khalid, M. Md. Estimation of above ground biomass
Methods Ecol. Evol. 9, 1928–1938 (2018). for oil palm plantations using allometric equations. 4th Int. Conf. Biol. Environ. Chem. 58,
52. del Hoyo, J., Elliott, A., Sargatal, D., Christie, D. & de Juana, E. Handbook of the Birds of the 110–114 (2013).
World Alive (Lynx Editions, 2015). 92. Zanne, A. E. et al. Data from: Towards a worldwide wood economics spectrum. Dryad
53. Westphal, C. et al. Measuring bee diversity in different European habitats and https://doi.org/10.5061/DRYAD.234 (2009).
biogeographical regions. Ecol. Monogr. 78, 653–671 (2008). 93. Krashevska, V. et al. Micro-decomposer communities and decomposition processes in
54. van der Vecht, J. The Vespinae of the Indo-Malayan and Papuan areas (Hymenoptera, tropical lowlands as affected by land use and litter type. Oecologia 187, 255–266 (2018).
Vespidae). Zool. Verh. 34, 1–82 (1957). 94. Ellsäßer, F. et al. Predicting tree sap flux and stomatal conductance from drone-recorded
55. Bohart, R. M. & Menke, A. S. Sphecid Wasps of the World: A Generic Revision (Univ. surface temperatures in a mixed agroforestry system—a machine learning approach.
California Press, 1976). Remote Sens. 12, 4070 (2020).
56. Yamane, S. A revision of the Japanese Eumenidae (Hymenoptera, Vespoidea) Insecta 95. Ellsäßer, F., Röll, A., Stiegler, C., Hendrayanto, & Hölscher, D. Introducing QWaterModel,
matsumurana. J. Res. Fac. Agric. Hokkaido Univ. 43, 1–189 (1990). a QGIS plugin for predicting evapotranspiration from land surface temperatures. Environ.
57. Goulet, H. & Huber, J. T. Hymenoptera of the World: An Identification Guide to Families Model. Softw. 130, 104739 (2020).
(Agriculture Canada, 1993). 96. Timmermans, W. J., Kustas, W. P. & Andreu, A. Utility of an automated thermal-based
58. Carpenter, J. & Nguyen, L. Keys to the genera of social wasps of South‐East Asia approach for monitoring evapotranspiration. Acta Geophys. 63, 1571–1608 (2015).
(Hymenoptera: Vespidae). Entomol. Sci. 6, 183–192 (2003). 97. Donfack, L. S. et al. Microclimate and land surface temperature in a biodiversity enriched
59. Choate, P. M. Key to the sub-Orders of Hymenoptera (Univ. Florida, 2011); https:// oil palm plantation. For. Ecol. Manag. 497, 119480 (2021).
entnemdept.ufl.edu/choate/hymenoptera.pdf. 98. Isles, P. D. F. The misuse of ratios in ecological stoichiometry. Ecology 101, e03153
60. Engel, M. S. The honey bees of Indonesia (Hymenoptera: Apidae). TREUBIA 39, 41–49 (2020).
(2012). 99. Ehbrecht, M., Schall, P., Ammer, C. & Seidel, D. Quantifying stand structural complexity
61. Kempson, D., Lloyd, M. & Ghelardi, R. A new extractor for woodland litter. Pedobiologia 3, and its relationship with forest management, tree species diversity and microclimate.
1–21 (1963). Agric. For. Meteorol. 242, 1–9 (2017).
62. Klarner, B. et al. Trophic niches, diversity and community composition of invertebrate top 100. Ehbrecht, M., Schall, P., Juchheim, J., Ammer, C. & Seidel, D. Effective number of layers:
predators (Chilopoda) as affected by conversion of tropical lowland rainforest in Sumatra a new measure for quantifying three-dimensional stand structure based on sampling with
(Indonesia). PLoS ONE 12, e0180915 (2017). terrestrial LiDAR. For. Ecol. Manag. 380, 212–223 (2016).
63. Potapov, A. M., Klarner, B., Sandmann, D., Widyastuti, R. & Scheu, S. Linking size 101. Willim, K. et al. Assessing understory complexity in beech-dominated forests (Fagus
spectrum, energy flux and trophic multifunctionality in soil food webs of tropical sylvatica L.) in Central Europe—from managed to primary forests. Sensors 19, 1684 (2019).
land-use systems. J. Anim. Ecol. 88, 1845–1859 (2019). 102. Khokthong, W. et al. Drone-based assessment of canopy cover for analyzing tree
64. Ehnes, R. B., Rall, B. C. & Brose, U. Phylogenetic grouping, curvature and metabolic mortality in an oil palm agroforest. Front. For. Glob. Change https://doi.org/10.3389/
scaling in terrestrial invertebrates. Ecol. Lett. 14, 993–1000 (2011). ffgc.2019.00012 (2019).
65. Ballauff, J. et al. Legacy effects overshadow tree diversity effects on soil fungal 103. Roswell, M., Dushoff, J. & Winfree, R. A conceptual guide to measuring species diversity.
communities in oil palm-enrichment plantations. Microorganisms 8, 1577 (2020). Oikos 130, 321–338 (2021).
104. Chao, A. et al. Rarefaction and extrapolation with Hill numbers: a framework for sampling K.D., I.G., A. Potapov, A.R., I.A., J.B., H.B., D. Berkelmann, S.B., D. Buchori, R.D., O.G., F.E., R.F.,
and estimation in species diversity studies. Ecol. Monogr. 84, 45–67 (2014). N.H., B.I., W.K., V.K., A.K., J.K., K.L., H.L., M.M., M.S.M., C.C.M.M., Y.A.M., G.B.P., H.D.P., A. Polle,
105. Jost, L. Partitioning diversity into independent alpha and beta components. Ecology 88, D.A.P., L. Sachsenmaier, S.S., F.S., A.C.S., L. Sundawati, T.T., M.W., D.H. and H.K.). Dorothea
2427–2439 (2007). Schlözer Postdoctoral Programme of the Georg-August-Universität Göttingen (N.G.-R.).
106. Hill, M. O. Diversity and evenness: a unifying notation and its consequences. Ecology 54,
427–432 (1973). Author contributions D.C.Z., N.G.-R., H.K. and D.H. conceived this work. D.C.Z. and N.G.-R.
107. Oksanen, J. et al. vegan: Community ecology package. R version 2.5-6 https:// developed the methodology. D.C.Z., N.G.-R. and M.S.M. created the software. D.C.Z., N.G.-R.
cran.r-project.org/web/packages/vegan/index.html (2018). and H.L. undertook the formal analysis. D.C.Z., F.B., K.D., A. Potapov, I.A., J.B., D. Berkelmann,
108. Byrnes, J. E. K. et al. Investigating the relationship between biodiversity and ecosystem S.B., F.E., R.F., N.H., W.K., V.K., A.K., J.K., K.L., H.L., C.C.M.M., H.D.P., D.A.P., L. Sachsenmaier,
multifunctionality: challenges and solutions. Methods Ecol. Evol. 5, 111–124 (2014). F.S. and C.A.S. conducted the investigations. D.S. obtained resources. D.C.Z. undertook data
109. Lefcheck, J. S. PiecewiseSEM: piecewise structural equation modelling in R for ecology, curation. D.C.Z., N.G.-R. and H.L. were responsible for visualization. H.K., D.H., M.W., I.G., H.B.,
evolution, and systematics. Methods Ecol. Evol. 7, 573–579 (2016). D. Buchori, R.D., O.G., M.M., Y.A.M., A. Polle, S.S., T.T., L. Sundawati and B.I. obtained funding.
D.C.Z., L. Sundawati and B.I. conducted the project administration. H.K., D.H., I.G., H.B.,
D. Buchori, R.D., O.G., M.M., Y.A.M., A. Polle, S.S., T.T., L. Sundawati, B.I. and A.R. were responsible
Acknowledgements We thank PT Humusindo for granting us access to and use of their for supervision. D.C.Z. and N.G.-R. wrote the original draft. F.B., K.D., I.G., A.P., A.R., G.B.P., D.C.,
properties, as well as A. F. C. Rodríguez, A. Angulo, L. S. Donfack, M. Ehbrecht and D. Seidel for M.W., H.K. and D.H. reviewed and edited the final manuscript.
their contribution to data; the field assistants E. J. Siahaan, K. H. Dalimunthe and Juliandi for
plot maintenance and for their technical assistance in data collection; and M. Teuscher and
Funding Open access funding provided by University of Neuchâtel.
A. Gérard for their contribution to establishing the experiment. Illustrations were provided by
F. Arndt, Formenorm Leipzig. This study was conducted using the research permit nos. 337/
SIP/FRP/E5/Dit.KI/IX/2016 (D.C.Z.), 58/SIP/IV/FR/10/2021 and 20/SIP/IV/FR/1/2022 (C.M.), 53/ Competing interests The authors declare no competing interests.
EXT/SIP/FRP/E5/Dit.KI/IX/2017 (V.K.), 2C13UA0011-S (H.L.), 54/EXT/SIP/FRP/E5/Dit.KI/IX/2017
(A. Potapov), 4233/FRP/E5/Dit.KI/XII/2017 (L. Sachsenmaier), 343/SIP/FRP/E5/Dit.KI/X/2016 Additional information
(A.K.), SK.336/KSDAE/SET/KSA.2/7/2019 (S.B. and C.A.S.), 1837/FRP/E5/Dit.KI/VIII/2017 (K.K.) Supplementary information The online version contains supplementary material available at
and 1092/FRP/SM/VIII/2015 and 1487/FRP/SM/KI VI/2016 (W.K.) from the Indonesia Ministry https://doi.org/10.1038/s41586-023-06086-5.
of Research and Technology. The experiment used in this study is part of the global Correspondence and requests for materials should be addressed to Delphine Clara Zemp.
network of tree diversity experiments TreeDivNet (https://treedivnet.ugent.be/). Deutsche Peer review information Nature thanks Bernhard Schmid, Robert Nasi and the other,
Forschungsgemeinschaft (DFG, German Research Foundation)—project ID 192626868—SFB anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer
990 and the Ministry of Research, Technology and Higher Education (Ristekdikti) in the reports are available.
framework of the collaborative German–Indonesian research project CRC990 (D.C.Z., F.B., Reprints and permissions information is available at http://www.nature.com/reprints.
Article

Extended Data Fig. 1 | Interaction between the experimental treatment indicate 95% confidence intervals with the centre line indicating marginal
(island, area, planted diversity) and biodiversity indicator. (a) Interaction means using ggeffect. Points are observed values (n = 56 study plots). Only
island × indicator. (b) Interaction area × indicator. (c) Interaction planted significant interaction terms in the analysis of variance of the linear mixed-
diversity × indicator. Lines/Solid points are linear mixed-effect model fits using effect models are shown (p < 0.05).
ggeffect: (a) centre for the error bars indicate marginal means (b-c) bands
 Oil palm dominance
Number of tree species

0 1 2 3 6

Plot size: PC 2
Edge lenght (m)

Oil
Oil palm cover
palm
5 10 20 40 density

ENL
SSCI
UCI
Tree density

Understorey Mean-
cover FRAC PC 1
Open Dense
Canopy gap Litter cover
fraction

Litter Tree
depth cover

Tree dominance
Extended Data Fig. 2 | Principal components of the vegetation structure in i.e., a measure of the geometric complexity of the vegetation structure that is
the tree islands. Each square represents one of the 52 tree islands, which vary density-dependent (1.20, −0.10), tree density (0.64, 0.23), UCI, i.e., understorey
in planted tree diversity (colour) and area (size of the square). The first two complexity index (0.93, 0.35), SSCI, i.e., stand structural complexity index
components (PC1 and PC2) explain 32% and 21% of the total variance, respectively. (0.85, 0.49) and oil palm density (0.64, 0.23). PC1 is mostly associated with dense
Vegetation structure variables included (variable loadings PC1 and PC2): Oil and structurally complex vegetation and with low understorey vegetation cover
palm cover (−0.13, 1.26), ENL, i.e., effective number of layers (−0.83, 0.65), (see Methods for details). PC2 is mostly associated with a high proportion of
understorey cover (−1.10, −0.07), canopy gap fraction (−1.12, −0.41), litter depth palms and low proportion of trees in the canopy. In the main text, we refer to
(−0.45, −0.79), tree cover (0.73, −0.99), litter cover (0.88, −0.55), MeanFRAC, PC1 as ‘structural complexity’ and minus PC2 as ‘tree dominance’.
Article

Extended Data Fig. 3 | Direct and indirect influence of tree island size and
planted tree diversity on multidiversity and multifunctionality. Results
of the structural equation models (i.e. standard coefficient estimates) for
multifunctionality and multidiversity (based on species richness, Shannon
diversity and Simpson diversity) with alternative calculation methods (average
or thresholding approach - see Methods). The effect of tree island area can be
direct (black bars) or via tree dominance (grey bars), and the effect of planted
tree diversity can be direct (dark green bars - here absent) or via structural
complexity (light green bars). Only significant estimates are presented
(p < 0.05 two-sided ANOVA, n = 52 tree islands).
Extended Data Fig. 4 | Changes in per area yield across four years for horizontal line is the mean of the reference palms (i.e. the palms in adjacent
different planted diversity levels. Relative changes in per area yield (∆Y ha/ position 3). ‘n’-values denote the number of months included in the annual
Y ha_ref in %) compared to conventionally managed oil palm monocultures boxes. Significant differences to reference as based on Mann-Whitney tests are
(‘reference’) for different planted tree species diversity: 0, 1, 2, 3, 6 species, and indicated at the bottom of each panel, with the significance levels * (p < 0.1),
all 52 tree islands combined (“total”). Plus signs represent means, horizontal ** (p < 0.05) and *** (p < 0.01).
lines medians, boxes interquartile ranges and vertical bars ranges. The dashed
Article

Extended Data Fig. 5 | Effects of the experimental treatment on per palm yield in the conventionally managed oil palm monocultures. c) There was no
and per area yield. Each dot represents an individual oil palm that was effect of the distance to the tree island edge (i.e. positions 1, 2 and 3) on the per
monitored within our study (N = 404 in total, from October 2017 to November palm yield of the adjacent oil palms (χ2 = 4.47, df = 2, p-value = 0.10). d) There
2018). In each boxplot, the horizontal line corresponds to the median; the lower was no effect of tree planting on the per palm yield of the adjacent oil palms
and upper hinges correspond to the first and third quartiles and the upper and (χ2 = 1.02, df = 1, p-value = 0.31). e) The effect of oil palm thinning on per
lower whiskers are limited by the 1.5 the interquartile ranges. The overall palm yield of the adjacent palms was marginally significant (χ2 = 3.78, df = 1,
effects of the experimental treatment on (a) per palm yield (in kg/palm) and p-value = 0.05). Multi-comparisons were conducted with post hoc Tukey tests.
(b) per area yield (in kg/ha) were significant in both cases (χ2 = 6.39, df = 2, More details of the statistical analysis are provided in section 2.3.2 of 19, where
p-value = 0.04 and χ2 = 75.35, df = 2, p-value < 0.001, respectively). Multi- models number one, two, three and four correspond here to the panels a), b),
comparison indicates that (a) per palm yield inside the islands was lower than c) and e), respectively. Significance levels are indicated by. (p < 0.1), * (p < 0.05),
per palm yield adjacent to the tree islands; (b) per area yield inside the tree ** (p < 0.01) and *** (p < 0.001).
islands were lower than per area yield adjacent to the tree islands and per area
Extended Data Table 1 | ANOVA of the linear mixed-effect models for biodiversity and ecosystem functioning

Each biodiversity and ecosystem functioning indicator was standardized to unit scale. Plot was included as a random term. Two-sided ANOVA (n = 56 study plots). numDF: degrees of freedom,
denDF: denominator degrees of freedom.
Article
Extended Data Table 2 | Summary statistics of piecewise structural equation models for each biodiversity and ecosystem
functioning indicators and multidiversity and multifunctionality

Every row is an individual model. Fisher’s C statistic, degrees of freedom (df), significance value based on a Chi-square test, Information criterion (Akaike, corrected Akaike and Bayesian), likeli-
hood degrees of freedom (K), sample size (n = 52 tree islands) and R-squared values.
Extended Data Table 3 | Result of the principal component
analysis of the vegetation structure variables.
(n = 52 tree islands)
Article
Extended Data Table 4 | List of soil fauna groups and
associated guild
 May 04, 2023

The data were processed and analyzed in R version 1.2.1335

Data and R code are available at https://doi.org/10.6084/m9.figshare.22320490

The raw data are available at https://data.goettingen-research-online.de/dataverse/crc990, with the specific link for each dataset provided in the Supplementary Tables 1 -
3. The processed data is available at https://doi.org/10.6084/m9.figshare.22320490.
Seed DNA sequences are available in NCBI Genbank under the accession numbers HYPERLINK "https://
www.ncbi.nlm.nih.gov/nuccore/OM811991.1/" OM811991- HYPERLINK "https://www.ncbi.nlm.nih.gov/nuccore/OM812021" OM812021, HYPERLINK "https://
www.ncbi.nlm.nih.gov/nuccore/OM837673" OM837673- HYPERLINK "https://www.ncbi.nlm.nih.gov/nuccore/OM837724" OM837724, and HYPERLINK "https://
www.ncbi.nlm.nih.gov/nuccore/OM935782" OM935782- HYPERLINK "https://www.ncbi.nlm.nih.gov/nuccore/OM935815" OM935815. Sequencing data of the soil
fungal community were deposited in the NCBI Sequence Read Archive (SRA) under Bioproject accession number PRJNA659225. The public UNITE database (https://
unite.ut.ee/) v7.2 on fungal ITS sequences was used as a reference of taxonomic classification. Sequence data of the bacterial communities were deposited in the NCBI
SRA under Bioproject accession number PRJNA841353. Sequence identification was performed by mapping all curated sequences against the SILVA database version 132
(https://www.arb-silva.de/).
 .....

Human research participants

Policy information about studies involving human research participants and Sex and Gender in Research.
Reporting on sex and gender N/A

Population characteristics N/A

Recruitment
N/A

Ethics oversight N/A

)

Note that full information on the approval of the study protocol must also be provided in the manuscript.

Field-specific reporting
Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before ma king your selection.

D Lite sciences D Behavioural & social sciences � Ecological, evolutionary & environmental sciences
For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf

Ecological, evolutionary & environmental sciences study design

Ali studies must disclose on these points even when the disclosure is negative.

Study description Our study was conducted in EFForTS-BEE, the Biodiversity Enrichment Experiment of the EFForTS project [Ecological and
Socioeconomic Functions of Tropical Lowland Rainforest Transformation Systems (Sumatra, lndonesia)]. EFForTS-BEE is part of the
global network of tree diversity experiments TreeDivNet (https://treedivnet.ugent.be/). The study region is characterized by a humid
tropical climate with a mean temperature of 26.7 ± 0.2 °C and an annual rainfall of 2,235 ± 381 mm and the dominant soil type is
loamy Acrisol. ln December 2013, 52 experimental plots (i.e. tree islands) were established in a conventional 140 ha oil palm
plantation. Following a random partition design, we systematically varied plot area (25, 100, 400 and 1600 m2) and tree species
diversity (O, 1, 2, 3 and 6 species). The six planted tree species (Archidendron jiringa (Jack) I.C.Nielsen. (Fabaceae), Parkia speciosa
Hassk. (Fabaceae), Durio zibethinus L. (Malvaceae), Dyera polyphylla (Miq.) Steenis (Apocynaceae), Shorea leprosula Miq.
(Dipterocarpaceae) and Peronema canescens Jack (Lamiaceae)) are native to the region and widely used for their fruits, timber or
latex. Around 40% of the oil palms located inside the tree islands were felled, with the number of felled oil palms differed depending
on the tree island area. The trees were planted between the felled and standing oil palms on a 2-m triangular grid. The tree islands
were fenced, and the management comprised a total stop of fertilizer, herbicide and pesticide application after planting. After May
2016, manual weeding was restricted to 1-m circles around the planted trees when these were shorter than the surrounding grass
layer, allowing for natural regeneration. ln addition to the 52 tree islands, we established four contrai plots in the conventional oil
palm plantation that were managed as usual, in the main text referred to as conventional monocultures. ln total, the study
comprised 56 plots. ln each study plot larger than 25 m2, one subplot of 5 m x 5 m was established in a random location at a
minimum distance of 1.5 m from the plot edge.
)

Research sample Ecosystem functioning

We measured 20 variables related to seven categories of ecosystem functioning including Productivity: (0) tree growth (basal area
increment of the planted trees in m2/ha/year) that was further excluded from the analysis - see method on multifunctionality, (1) Oil
palm yield (per island oil palm yield changes in kg of fresh fruit bunches / island), (2) aboveground biomass (biomass stored in the
aerial parts of the planted trees and the oil palms, in t/ha), Resistance to invasion: (3) native seeds (total number of arriving native
seeds / m2); (4) resistance to invasive plants (100- observed cover of Clidemia hirta, in %); Pollination: (5) pollinators (number of
sampled individuals), (6) pollination rate (fraction of flowers on phytometer plants that are pollinated, %), Soil quality: (7) soil P
(phosphorous content, %), (8) 1 / soil C:N (that is the molar ratio of soil C content to soil N content), (9) soil decompaction (inverse of
soil bulk density in g/cm3); Predation and herbivory: (10) predatory invertebrates (total activity duration of insectivorous bats and
birds, in seconds); (11) predatory arthropods (number of sampled individuals), (12) predatory soil fauna (energy flux, in J/hour), (13)
herbivory (energy flux, in J/hour); Carbon & nutrient cycling: (14) Decomposers (energy flux, in J/hour); (15) litter decomposition
(relative biomass loss of litter after 6 months in litterbags, %), (16) litter input (biomass of leaf litter falling in traps, g / m2), water and
climate regulation: (17) evapotranspiration (canopy latent heat flux, in W/m2); (18) soil water infiltration capacity (saturated soil
 57

As tree growth and leaf litter input were correlated and formed cluster, we excluded tree growth from the analysis. For leaf litter
input, we then excluded outliers defined as plot-level values outside the range of 3 standard deviations around the median of the
entire data (less than 5% of the litterweight data, total and per species).

at https://doi.org/10.6084/m9.figshare.22320490

Each variable presented in the main text had one measurement per plot, such that randomization was not applicable.

Each variable presented in the main text had one measurement per plot, such that blinding was not applicable.

◦ ◦

https://doi.org/10.6084/m9.figshare.22320490
 Materials & experimental systems Methods
n/a lnvolved in the study n/a lnvolved in the study
IZ! D Antibodies IZ! D ChlP-seq

IZ! D Eukaryotic cell lines IZ! D Flow cytometry

IZ! D Palaeontology and archaeology IZ! D MRl-based neuroimaging

D IZ! Animais and other organisms
IZ! D Clinical data
IZ! D Dual use research of concern

Animais and other research organisms

Policy information about studies involving animais: ARRIVE guidelines recommended for reporting animal research, and Sex and Gender in
Research

Laboratory animais No laboratory animais were used for this study.

-'
Only invertebrate animals (arthropods and earthworms) were collected and killed using ethanol during the study. This was
Wild animais
necessary to assess biomass/density and community composition. Soil animals were extracted from soil samples using a Kempson
extractor under a heat gradient of 40-45 grad Celsius above and 15 grad Celsius below the samples. Animals were first collected in
dimethyleneglycol - water solution (1:1) and thereafter transferred into 70-80% ethanol solution. Understorey arthropods were
collected using pan traps that stayed for 45 hours in the field. Arthropods where directly killed in the bowl that is filled with
water and one drop of scentless soap. Each pan trap was shaken off through a sieve so that the arthropods could be collected in
the sieve and then stored in a test tube containing 70% ethanol. The animals were collected by Anton Potapov with the research
permit: 349/SIP/FRP/ES/Dit.KI/X/2016 (Validity 5 October 2016 to 5 October 2017) and 54/EXT/SIP/FRP/ES/Dit.KI/IX/2017 )
(validity 04 October 2017 to 4 October 2018) and Isabelle Arimond with the research permit 370/SIP/FRP/ES/Dit. KI/X/2016
Reporting on sex
(Validity 20 October 2016 to 20 March 2017). Animal sex was not considered in the study.
After each round of fieldwork, samples were taken to the lab for identifications. Collected soil samples were transported in the
Field-collected samples lab within 2-3 days for heat extraction. All extracted invertebrates were stored in 70-80% ethanol solution and sorted to high-
rank taxa under dissecting microscope.

Ethics oversight No ethical approval was required for this study as it complies with regulations of the Collaborative Research Center 990 (project ID
192626868 - SFB 990)
Note that full information on the approval of the study protocol must also be provided in the manuscript.