The biodiversity enrichment experiment
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, Indonesia))19. EFForTS-BEE is part of the global network of tree diversity experiments TreeDivNet45 (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 Acrisol46. In December 2013, 52 experimental plots (tree islands) were established in a conventionally managed 140 ha oil palm plantation. Following a random partition design47, we systematically varied island area (25, 100, 400 and 1,600 m2) and planted diversity (zero, one, two, three and six tree 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 latex35. Around 40% of the oil palms located inside the tree islands were felled, with the number of felled oil palms differing depending on the tree island area33. 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. In addition to the 52 tree islands, we established four control plots in the oil palm plantation that have a fixed area (100 m2) and that were managed conventionally (that is, no oil palm was felled, no tree was planted and application of fertilizer, herbicide and pesticide was as usual), in the main text referred to as conventionally oil palm monocultures. In total, the experiment comprises 56 study plots19. In each study plot larger than 25 m2, one subplot of 5 × 5 m2 was established in a random location at a minimum distance of 1.5 m from the plot edge.
Field measurements
We conducted an interdisciplinary field campaign from October 2016 to October 2018, that is, 33–57 months after establishment of the experiment. At this early stage of the experiment, the tree islands already differed in their structural complexity30 and the planted trees reached up to 16 m height35. In all the 56 study plots, several indicators related to biodiversity, ecosystem functioning and structure were measured using standardized procedures and constant sampling areas at the level of the plot or subplot (Supplementary Tables 1–3). Only trees were sampled at unequal areas (that is, all trees present in the plots were sampled) and were therefore standardized using rarefaction curves (see ‘Trees’). Oil palm yield was continuously measured since the beginning of the experiment at the level of individual palm but the data were then aggregated over space and time (see ‘Per area yield’ and ‘Per island yield’). Each variable presented in the main text had one measurement per plot, such that blinding and randomization were not applicable. No statistical methods were used to predetermine sample size. The data were processed and analysed in R v.1.2.1335 (ref. 48).
Birds and bats
We recorded audible and ultrasound in March 2017 using automated sound recorders (SM2Bat+ recorders, Wildlife acoustics; with an acoustic SMX-II microphone on the left channel and one full-spectrum Sonitor Parus49 microphone on the right channel), strapped to wooden poles at a height of 1.5 m in the centre of the plot. On consecutive days, we extracted sound recordings for sampling birds and insectivorous bats. We used two stereo 15 min recordings starting 15 min before sunset and two 15 min stereo recordings starting at sunrise, sampled at 22.05 kHz for birds. We used two 40 min mono sound recordings from the right channel, extracted from consecutive nights, starting 20 min after sunset, sampled at 384 kHz for bats. Twelve sound recorders were installed simultaneously in 12 randomly chosen plots. The recordings were annotated in ecoSound-web50 to extract the duration of each bird vocalization and bat pass and bird detection distances were estimated using reference sound transmission sequences51. We assigned birds to species according to Birdlife International taxonomy. Owing to the lack of standard protocols and reference collections for Southeast Asia, we could not identify bats to species and used sonotypes instead. We appended feeding guild information to each bird species52 all detectable bats were echolocating and thus considered insectivorous bats. Only bird vocalizations detected within a 28 m radius were included, which corresponds to the diameter of the largest study plot (40 × 40 m2). We used the maximum number of individuals detected simultaneously in all recordings per plot as a conservative proxy of abundance per bird species or bat sonotype.
Understorey arthropods
Arthropods were sampled in the understorey vegetation during October 2016 to January 2017. Each plot was sampled three times with six pan traps per plot exposed for 45 h. Traps were made of white plastic soup bowls covered with yellow ultraviolet spray-paint53 and were filled with water and one drop of regular soap. They were fixed in a holding system in groups of three at the height of the surrounding plants and these systems then equally distributed in distance from edge and to each other. All arthropods were preserved in 70% ethanol. Subsequently, all individuals were identified to higher taxonomic groups and morphospecies. The taxonomic groups Hymenoptera, Lepidoptera and Araneae were categorized into functional groups (pollinators, predators and parasitoids) using different identification keys54,55,56,57,58,59,60 (Supplementary Table 10). Predators and parasitoids were merged into the single functional group ‘predators’.
Soil fauna
During October–November 2016, in each plot, four soil samples of 16 × 16 cm2 were taken randomly within the subplot with a spade. Samples included litter (if present) and soil down to a depth of 5 cm. Animals were extracted using a gradient heat extractor61 and collected in dimethyleneglycol–water solution (1:1) and thereafter transferred into 70% ethanol62. All extracted animals were counted and sorted into 28 taxonomic groups (in most cases orders) allowing for functional group classifications63; Extended Data Table 4. We calculated community metabolism of all animals that were classified as detritivores, herbivores and predators in a sample by using mean group- and ecosystem-specific estimates derived from ref. 63. The estimates are based on measurements of more than 5,000 individuals of soil animals across eight different oil palm plantations in the same region; to estimate community metabolism, individual body masses were recalculated to metabolic rates using group-specific regressions from ref. 64. Community metabolism was calculated by summing up metabolic rates of all individuals; we used the mean per plot across four samples for each functional group (detritivores, herbivores and predators) for the analysis. We also computed taxonomic diversity as the number of taxonomic groups (in most cases orders) present in each plot for the analysis.
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 before soil collection. The soil was sieved through a 50 × 50 mm2 sieve and roots were separated from soil. The fungal community was assessed using Illumina next-generation sequencing (Illumina) of the ITS2 marker region. The detailed protocol for amplification, amplicon sequencing and generation of fungal operational taxonomic units (OTU) is described in ref. 65. OTUs were classified taxonomically using the BLAST (blastn, v.2.7.1) algorithm66 and the UNITE v.7.2 (UNITE_public_01.12.2017.fasta) reference database67.
Prokaryotes
In May 2017 three cores of topsoil (10 cm depth) were taken in each subplot. Soil cores were then mixed, homogenized and freed from roots. A total of 5 ml of RNAprotect Bacteria reagent (Qiagen) was added to 5 g of soils to prevent nucleic acid degradation. DNA and RNA were extracted from 1 g of soil by using the Qiagen RNeasy PowerSoil Total RNA Kit and the RNeasy PowerSoil DNA Elution Kit (Qiagen). The V3–V4 region of the 16S ribosomal RNA gene was amplified and sequenced as described in ref. 68. Paired-end sequences were quality filtered with fastp (v.0.20.0)69 and merged with PEAR v.0.9.11 (ref. 70). Remaining primer sequences were clipped with cutadapt v.2.5 (ref. 71). Size filtering, dereplication, denoising and chimaera removal was performed with vsearch v.2.12.0 (ref. 72). Curated sequences were then classified by mapping each sequence against the SILVA database with the BLAST73. Counts were normalized by using the GMPR normalization74.
Seeds
We installed four seed traps in each of the 56 study plots for 1 yr, that is, between 1 April 2017 and 29 March 2018. The traps were built using fine-mesh cloth attached to a squared structure made in PVC pipes of size 50 cm × 50 cm fixed at 1 m from the ground. The traps were installed at random locations in each of the four quadrants of each plot, at a minimum distance of 1 m from the plot edge. The contents of the traps were collected twice a month, dried at approximately 40 °C during 3–7 days. All the seeds were carefully extracted from the samples, counted and separated by morphospecies using hand lenses (×10 magnification) and a microscope (Leica photomicroscope with ×400 magnification) for very small seeds. Molecular identification of the morphospecies was implemented using three universal plant DNA barcodes (matK, rbcL and ITS2)75,76,77,78 and taxonomic assignments were made using BLASTn search against the NCBI Genbank reference sequence database79. Sequences obtained from the barcode loci were deposited in NCBI Genbank under the accession numbers OM811991–OM812021, OM837673–OM837724 and OM935782–OM935815. We classified each morphotype as native or non-native species using available literature80,81 (http://www.plantsoftheworldonline.org/). We derived the native seed density (number of native seeds per m2) as the total number of native seeds over the entire sampling duration per plot, which was used as an indicator of ecosystem functioning (see ‘Ecosystem functioning’). Seed diversity, calculated on the basis of the Hill number frameworks and used as indicators of biodiversity (see ‘Biodiversity’), was derived from the pooled samples per plot over the entire sampling duration for all seeds (native and non-native).
Herbs
All non-woody terrestrial vascular plants (for example, angiosperm herbs and vines, ferns, but not epiphytes) in the subplot were recorded from February until March 2018. They were classified to species or morphospecies and herbaceous cover (in absolute per cent ratios from 1% to 100%) was estimated by two people. Epiphytes growing on the stems of trees or palms were excluded, whereas vine species that rooted in the ground and climbed up stems of trees or palms were included. Herbarium specimens were collected and stored in the laboratory of Jambi University. All names were checked following The Plant List 2013, v.1.1 (http://www.theplantlist.org).
Trees
All planted trees were surveyed in January to February 2018 as part of a yearly inventory35. Furthermore, we surveyed all free-standing woody plants (trees, shrubs and bamboos) that colonized the plots with a length of ≥130 cm from April until August 2018. For each species or morphospecies, one voucher specimen was collected, dried and pressed according to standard procedure. In the main text, we refer to the colonized woody plants as ‘trees’, unless stated otherwise. Because the number of sampled trees largely varied according to the tree island area, we standardized the diversity estimates using rarefaction curves (R package iNEXT)82 to 24 individuals, which represent the median number of individuals per plot.
Pollen
To collect pollen/spore rain, Behling pollen traps83 were installed from June until October 2018. Each trap consists of a plastic tube which is placed about 30 cm above the ground and is held by a fixing pole. The tube is filled with 5 ml of liquid glycerol, synthetic cotton and, on the top, it is covered by a mosquito net to reduce disturbance from animals or litter and prevent the cotton from being removed. In tropical regions heavy rainfalls occur, thus it is necessary to prevent the pollen from pouring out of the pollen trap. In the Behling trap, glycerol is used, which has a higher density compared to water. Consequently, the incoming rainfall can flow out of the trap without taking the pollen, which is trapped in the synthetic cotton and in the glycerol83. The Behling traps were modified to mimic the surrounding environment and maximize recovery. In total, 168 pollen traps were installed in the 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 processed and analysed. Firstly, each pollen trap was washed with distilled water through a 2 mm mesh sieve to remove large size materials. Afterwards, the pollen traps were sieved through a 150 µm mesh sieve to exclude medium-sized materials from the samples. Two Lycopodium tablets were added as markers to each sample to estimate palynomorph concentrations84 and the Erdtman acetolysis85 was applied, to remove cellulose material. Residues were mounted in glycerol jelly for pollen visualization, identification and counting. Pollen and spore analyses were carried out using light microscopy. All identified pollen and spore types were photographed using a Leica photomicroscope with ×400 magnification. For each trap, a total sum of at least 100 pollen grains were counted. Pollen and spore grains can rarely be identified to species level and the level of taxonomic identification varies for different groups of plants. Consequently, a reduction to the family level has been proposed for studies involving analysis of palynological diversity in the tropics86.
Pollination
We assessed pollination rate on chilli pepper plants (Capsicum annuum) as phytometer plants, selected for potential shade tolerance87, widespread home garden cultivation in this region88 and the potential role pollination can play in fruit quality and yield89. We raised 1,500 individuals of a locally available variety of C. annuum from seed outside of the study plots. During the growth period outside the study plots, we applied NPK fertilizer and pesticide (imidacloprid, deltamethrin, mancozeb and abamectin) following local practices to standardize growing conditions and control pest damage before transfer to field sites. We halted fertilizer and pesticide application 1 week before placement in the plots and only watered as conditions required thereafter. In February 2018, we selected 224 healthy individuals of comparable size to transfer to the 56 study plots (four plants per plot). The four chilli plants were placed, still in their pots, at the centre of each plot for 5 weeks for a period of open pollination and monitoring, followed by 3 weeks for fruit harvesting. We removed any flowers before placement in the field, so pollinated flowers and developing fruits were assumed to result from pollination within the study plots. During the period before final harvest, each plot was revisited once per week and the number of flowers were counted per plant. The rate of successful pollination was estimated from the fruit to flower ratio, which was the total number of harvested fruits divided by the total number of flowers observed per plot.
Per palm yield
We followed the conventional harvesting procedure established by the plantation manager of PT Humusindo and measured the weight of the fresh fruit bunches directly after harvest using a portable scale. We measured the per palm yield (kg per palm) of all palms inside the 52 tree islands (N = 214) and one palm per control plot with conventional management (N = 4). To obtain a more solid estimate for the conventional plantation, we measured the per palm yield of 30 more reference palms that were evenly distributed across the conventional plantation at approximately an equal distance to each tree island and whose neighbourhood is characteristic of conventionally managed oil palm monocultures (Supplementary Note 2 and Supplementary Figs. 1–3). To examine potential changes in yield in the conventionally managed oil palm plantation surrounding the tree islands (‘spillover effects’), we measured the per palm yield of three oil palms adjacent to each tree island, at increasing distance to the island’s edge (at position number 1, 2 and 3)33. For direct comparison with earlier findings, we analysed the data following established methodology33. The tests were based on linear mixed-effect models with the annual yield of individual oil palms as the response variable and the plot identity as the random effect. Pairwise comparison was conducted with a post hoc Tukey test. Because our results indicate that only the palm directly adjacent to the tree island (in position number 1) was affected by the experimental treatment (Extended Data Fig. 5a, in agreement with ref. 33), we do not consider the palms in position 2 and 3 in the yield calculation per island (Per island yield). The per palm yields in and adjacent to the 56 study plots have been continuously monitored since the establishment of the experiment in December 2013. The extra 30 palms were established in December 2016. For consistency with other indicators (Extended Data Tables 1–3), we reported yield data for 1 yr (November 2017 until 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 shown in Extended Data Fig. 4, in which the oil palms in position 3 were used as reference palms for the corresponding time period.
Per area yield
We estimated per area yield (∆Yha, kg ha−1) as the yield of a given palm (kg palm−1) multiplied by a stand density-dependent expansion factor (EF) to derive estimates of per area yield (kg ha−1). We then calculated the per area yield change between tree islands and reference (kg ha−1; Supplementary Note 3). This approach accounts for changes in per area yield due to oil palm thinning (that is, reduced oil palm densities and changes in per palm yield in the tree islands) but does not account for potential changes in per palm yield on the surrounding plantation, for example, because of spillover effects19. An alternative analysis considering spillover effects was performed at the plot level (Per island yield).
Per island yield
We estimated oil palm yield changes at the tree island scale (∆YIsland, in kg island−1); equations (1)–(4) following established methodology33. This method considers the yield foregone owing to the removal of some oil palms before the experiment, as well as changes in per palm yield inside the tree islands and directly adjacent to the tree islands (at position 1, that is, spillover effects). Because the number of oil palms inside and adjacent to the tree islands and the number of 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 method was initially designed to calculate the oil palm yield changes for the tree islands, here we also apply it to the four control plots to integrate them in our synthesis analysis.
$${varDelta Y}_{{rm{Island}}}={Y}_{{rm{Spillover}}}+{Y}_{{rm{RemainChange}}}-{Y}_{{rm{Foregone}}}$$
(1)
$${Y}_{{rm{F}}{rm{o}}{rm{r}}{rm{e}}{rm{g}}{rm{o}}{rm{n}}{rm{e}}}={N}_{{rm{f}}{rm{e}}{rm{l}}{rm{l}}{rm{e}}{rm{d}}}times {Y}_{{rm{p}}{rm{_}}{rm{r}}{rm{e}}{rm{f}}}$$
(2)
$${Y}_{{rm{R}}{rm{e}}{rm{m}}{rm{a}}{rm{i}}{rm{n}}{rm{C}}{rm{h}}{rm{a}}{rm{n}}{rm{g}}{rm{e}}}={N}_{{rm{i}}{rm{n}}}times ({Y}_{{rm{p}}{rm{_}}{rm{i}}{rm{n}}}-{Y}_{{rm{p}}{rm{_}}{rm{r}}{rm{e}}{rm{f}}})$$
(3)
$${Y}_{{rm{Spillover}}}={N}_{{rm{adj}}}times left({Y}_{{rm{p}}_{rm{adj}}}-{Y}_{{rm{p}}_{rm{ref}}}right)$$
(4)
∆YIsland, per island oil palm yield change (kg island−1)
YSpillover, per island yield changes due to spillover effects (kg island−1)
YRemainChange, per island yield changes inside the island (kg island−1)
YForegone, per island yield foregone due to oil palm removal (kg island−1)
Ni, number of remaining oil palms inside the island
Nad, number of oil palms directly adjacent to the island (that is, adjacent position 1)
Nfelled, number of removed oil palms in the island
Yp_ref, median per palm yield of the reference palms in the conventionally managed oil palm plantation (kg palm−1)
Yp_in median per palm yield inside the tree island (kg palm−1)
Yp_adj, per palm yield directly adjacent to the tree island (that is, adjacent position 1 (kg palm−1)).
Above-ground biomass
For all the planted trees, we measured the basal diameter (at 10 cm above ground), the diameter at breast height (130 cm above ground) and the tree height in 2017 as part of a yearly inventory35. In January and February 2017, we also measured the height of the oil palms at meristem, that is, the point of attachment of the young leaves to the oil palm trunk30. We estimated above-ground biomass of the trees (equation (5)) and the oil palms (equation (6)) using the respective allometric equations of refs. 90,91:
$${{rm{AGB}}}_{{rm{tree}}}=0.0673times {(rho times {{rm{DBH}}}^{2}times H)}^{0.976}$$
(5)
$${{rm{AGB}}}_{{rm{palm}}}=71.797times H-7.0872$$
(6)
AGBtree, above-ground biomass of the planted trees (kg tree−1)
AGBpalm, above-ground biomass of the oil palms (kg palm−1)
DBH, tree diameter at breast height (cm)
H, height of tree or palm (m)
ρ, wood density (g cm−3).
Wood density for Peronema canescens (0.61 g cm−3), Parkia speciosa (0.54 g cm−3) and Dyera polyphylla (0.36 g cm−3) was based on EFForTS core plot data, whereas for Archidendron sp. (0.36 g cm−3), Shorea leprosula (0.44 g cm−3) and Durio zibethinus (0.516 g cm−3) it was taken from the global wood density database92.
We estimated the total above-ground biomass per plot as the sum of the above-ground biomass of the palms and the planted trees (equation (7)). The estimations of total AGB did not consider the necromass, litter, understorey vegetation and spontaneously established trees, which were considered negligible.
$${rm{A}}{rm{G}}{rm{B}}=(,sum {{rm{A}}{rm{G}}{rm{B}}}_{{rm{t}}{rm{r}}{rm{e}}{rm{e}}}+sum {{rm{A}}{rm{G}}{rm{B}}}_{{rm{p}}{rm{a}}{rm{l}}{rm{m}}},/{N}_{{rm{i}}{rm{n}}}times {d}_{{rm{p}}{rm{a}}{rm{l}}{rm{m}}},)/(Atimes 10,000)$$
(7)
AGB, total above-ground biomass per plot (t ha−1)
dpalm, density of oil palms (number of oil palms per ha) that takes into account the local neighbourhoods of the plots (also referred to as EF; see Supplementary Note 3)
A, area of the plot (m2).
Tree growth
The growth of the planted trees per plot was calculated as35:
$${{rm{B}}{rm{A}}}_{{rm{i}}{rm{n}}{rm{c}},2017-2018}=sum ({{rm{B}}{rm{A}}}_{{rm{t}}{rm{r}}{rm{e}}{rm{e}},2018}-sum {{rm{B}}{rm{A}}}_{{rm{t}}{rm{r}}{rm{e}}{rm{e}},2017})/A.$$
(8)
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)
BAtree, year, tree basal area (in cm2) derived from the basal diameter (cm) in the specific year.
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 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 also sorted the leaves by species and weighted the content for the six planted tree species and oil palm separately.
For each sampling date, we aggregated the values at plot level using the median per plot of the litter weight. 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 litter weight data, total and per species). To get annual estimates, we summed the available plot-level values over time and divided them by the number of sampling dates (that is, between 17 and 24, depending on the number of missing traps or excluded outliers). We then multiplied the obtained values by the seed trap area (0.25 m2) to get the leaf litter fall in g m2 yr.
Leaf litter decomposition
We installed litterbags (20 × 20 cm2, 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 freshly fallen air-dried leaf litter for each tree species or their combinations in experimental plot. In each plot, one litterbag was installed in November 2017 for a duration of 6 months. Decomposition (litter mass loss) was calculated as the difference between the initial litter dry mass and litter dry mass remaining after 6 months and expressed as a percentage of decomposed material.
Water infiltration capacity
To quantify soil water infiltration capacity, we measured saturated soil hydraulic conductivity (Kfs, cm h−1) using a dual-head infiltrometer (Saturo) in March 2018 near the subplot centre in 35 (out of the 52) tree islands and in the four control plots representing conventionally managed oil palm monocultures. Owing to a broken instrument, the 17 remaining plots were measured using a custom manual double-ring infiltrometer, which tends to yield higher Kfs estimates than the dual-head approach because there is no correction for lateral flow. In three plots, Kfs was measured with both devices. We plotted these values against each other and found a close linear relationship (R² = 0.98, P = 0.066); even though it was only marginally significant because of the small sample size, we used it to correct the values from the 17 plots that were measured manually (Kfs_corr = 1.44 + 0.55 Kfs_double_ring) to allow for comparability across all 56 plots.
Evapotranspiration
We recorded land and canopy surface temperatures using a radiometric thermal camera (FLIR Tau 2 640, FLIR Systems) attached to a TeAx ThermalCapture module (TeAx Technology GmbH) mounted on a multicopter drone (MK EASY Okto V3; HiSystemsy) as described in ref. 94. Image sets were recorded four to five times per day around noon, covering each plot once over a 9 day period encompassing varying weather conditions. Land and canopy surface temperatures were the main input for modelling latent heat flux (in W m−2) and deriving evapotranspiration using the QWaterModel QGIS3 Plugin95, which is based on the DATTUTDUT energy balance model96. Measured short-wave radiation and relative humidity were used as further input variables to support the prediction of latent heat flux and derive evapotranspiration.
Microclimate
We measured microclimate using temperature per humidity loggers: hydrochron (DS1923-F5) and thermochron (DS1922L-F5) iButtons, Maxim integrated. The loggers were installed in the middle of each plot at 1.5 m above ground and were protected from water and direct solar radiation using handmade multiplate radiation shields97. Data were collected for 1 yr (18 November 2017 until 19 September 2018) every 3 h starting at midnight. As a proxy for microclimate buffering, we calculated the daily amplitude as the absolute difference between values at 7:00 and 15:00 (ref. 97), aggregated using the median value over the entire measurement period.
Soil properties
We determined soil total carbon content (g mg−1), total nitrogen content (mg g−1) and plant available phosphorus content (mg g−1) using the same three soil samples as for fungi community data (see ‘Fungi’) and the method of determination is described in detail in ref. 65. We then calculated the C:N ratio accounting for the molar mass of the elements following ref. 98, that is, 12.0107 for carbon and 14.0067 for nitrogen. We also measured soil bulk density (g cm−3) using five soil samples taken in the subplot in May 2018. Soil rings of 100 cm3 were inserted horizontally into the first 5 cm of topsoil. The soil was weighed, dried at 105 °C until constant weight and weighted again. Calculation was done on the dry weight basis, for which the sample dry weight (g) was divided by the volume of the sample (cm3) collected from the average of the five replicates. We used the mean per plot for all mentioned soil variables and used the inverse of C:N ratio and soil bulk density as measures of soil fertility and soil decompaction, respectively (Supplementary Table 1).
Vegetation structure
We measured 12 variables representing various aspects of the vegetation structure (Supplementary Table 3). We used a terrestrial laser scanner Focus M70 (Faro Technologies) to create three-dimensional point clouds of the vegetation at the centre of each plot in September and October 2016, as described in ref. 30. We computed the (1) stand structural complexity index (SSCI) following ref. 99 and its two components: (2) the mean fractal dimension index (MeanFRAC) derived from cross-sections of polygons in the three-dimensional point cloud, which is a scale-independent and density-dependent measure of structural complexity and (3) the effective number of layers (ENL) that describes vertical stratification based on the Simpson Index100. ENL and MeanFRAC are integrated in the SSCI and all these three measures were derived from vegetation parts above 130 cm. We also derived (4) the understory complexity index that measures the fractal dimension of horizontal cross-sections of the point cloud between 80 and 180 cm height, thereby measuring the structural complexity of the understorey vegetation101. (5) Canopy gap fraction was estimated from hemispherical photographs at plot level as described in ref. 30. Drone-based photogrammetry dated from September to October 2016 was used to further partition the canopy (in %) as (6) oil palm cover and (7) tree cover as described in ref. 102. We also used the drone-based 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 the planting scheme (Supplementary Figs. 2 and 3) . For the smaller plots (5 × 5 m) unaffected by thinning, the oil palm density was simply the typical planting density in conventionally managed oil palm plantations (120 planted palms per hectare). Further details on the oil palm density calculation are given in the Supplementary Note 3. We also calculated (9) tree density as the number of trees planted and from natural regeneration per plot and expressed per hectare. We estimated the portion of the ground (in percent) as (10) understorey vegetation cover and (11) litter cover per subplot in February–March 2018. The understorey vegetation cover included all parts of plants lower than 130 cm in height, including the trunks and other parts of the planted trees but excluding oil palm trunks. (12) The litter depth (cm) was measured as the mean value in three randomly chosen positions inside each subplot with a metal ruler. To extract orthogonal axes (PC1 and PC2) that represent most of the variability in the vegetation structure, we applied a principal component analysis on all the structural variables after standardization to zero mean and unit variance.
Restoration outcomes
Ecosystem functioning
We measured 20 variables related to seven categories of ecosystem functioning including: productivity as (0) tree growth (basal area increment of the planted trees in m2 ha−1 yr−1) that was further excluded from the analysis—see Supplementary Fig. 4, (1) oil 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 parts of the planted trees and the oil palms, in t ha−1); resistance to invasion as (3) native seeds (total number of arriving native seeds per m2) and (4) resistance to invasive plants (100—observed cover of Clidemia hirta, in %); pollination as (5) pollinators (number of sampled individuals) and (6) pollination rate (fraction of flowers on phytometer plants that are pollinated, %); soil quality as (7) soil P (phosphorous content, %), (8) 1/soil C:N (the molar ratio of soil C to soil N concentration) and (9) soil decompaction (inverse of soil bulk density in g cm−3); predation and herbivory as (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 h−1), (13) herbivory (energy flux, in J h−1); carbon and nutrient cycling as (14) decomposers (energy flux, in J h−1); (15) litter decomposition (relative biomass loss of litter after 6 months in litterbags, %) and (16) litter input (biomass of leaf litter falling in traps, g m−2); water and climate regulation as (17) evapotranspiration (canopy latent heat flux, in W m−2); (18) soil water infiltration capacity (saturated soil hydraulic conductivity in cm h−1) and (19) microclimate buffering (median daily amplitude of air temperature during 1 yr, °C d−1). A more detailed summary of the 20 ecosystem functioning variables is presented in Supplementary Table 2.
Biodiversity
We derived taxonomic diversity for soil bacteria and soil fungi, soil fauna, herbs, trees, seeds, pollen, understorey arthropods, birds and bats. Most of the groups (arthropods, herbs, trees, birds and seeds) were sorted at the lowest possible taxonomic level (species or morphospecies). Pollen, soil fauna and bats were sorted to higher levels, mainly family, order and morphotypes, respectively. Soil bacteria and soil fungi were analysed by DNA-based marker gene sequencing as amplicons sequence variants or OTU, respectively. Hereafter, we refer to these different taxonomic units (species, family, order, morphotypes and OTU) as ‘species’ for simplicity.
Diversity was measured following the Hill number framework, which allows comparison across diversity indices that weigh relative abundances to varying extents (species richness, Shannon diversity and Simpson diversity) and are expressed in terms of effective numbers of species103,104,105,106. Species richness is more sensitive to locally rare species, Simpson diversity is more sensitive to locally dominant species and Shannon diversity favours neither rare nor dominant species. We show results for species richness and Simpson diversity in the main text and for all indicators in the Extended Data Tables 1 and 2 and Supplementary Tables 4– 8. The calculations were performed using the R packages iNext82 and vegan107.
Multidiversity and multifunctionality
Different indicators of biodiversity and ecosystem functioning were aggregated by calculating multidiversity and multifunctionality, respectively. Following ref. 18, we performed a cluster analysis to preselect indicators for achieving a representative measure of ‘ecosystem function multifunctionality’. As tree growth and litter input were correlated and formed a cluster, we excluded tree growth from the analysis (Supplementary Note 4 and Supplementary Fig. 4). Following a threshold approach108, we calculated multifunctionality (and multidiversity) as the number of ecosystem functioning (and biodiversity) indicators that cross a threshold, expressed as a certain percentage of the maximum observed values in our study landscape (among all 56 study plots). We calculated multifunctionality and multidiversity for all thresholds from 1% to 99% and presented results for a 50% threshold in the main text. To reduce the influence of extreme values, we used the mean of the three highest values observed in all study plots, respectively. As an alternative to the threshold approach, we also calculated multidiversity and multifunctionality as the average of the indicators108. Before multidiversity and multifunctionality calculations, all the variables were standardized to unit scale (for biodiversity and ecosystem functioning separately). The calculations were performed using the package multifunc in R108.
Statistical analysis
Linear mixed-effect models
We used linear mixed-effect models to test the effects of the experimental treatment on restoration outcomes. We fitted three separated models for biodiversity using species richness, Shannon and Simpson diversity as response variables and one model for ecosystem functioning. These models included tree island (compared to our controls of conventionally managed oil palm monocultures), island area (plot edge length in m), planted diversity and the restoration outcome (either biodiversity or ecosystem functioning indicators) as single factors and tree island × indicator, island area × indicator, planted diversity × indicator, island area × planted diversity and island area × planted diversity × indicator interactions. For conventionally managed oil palm plots, island area was set to 10 m edge length and planted diversity to zero. Each response variable (biodiversity and ecosystem functioning indicators) was standardized to unit scale (between 0 and 1) as this improved the model diagnostics before applying the respective linear mixed-effect models; whereas we used logarithmic transformations for island area and planted diversity. Plot was included as a random term.
As an alternative to the linear mixed-effect models, we applied Kruskal–Wallis tests on each indicator of biodiversity and ecosystem functioning for comparison between the 52 tree islands and the four conventionally managed oil palm monocultures as control plots (Supplementary Note 5 and Supplementary Tables 8 and 9).
Structural equation modelling
We used piecewiseSEM109 to assess the influence of tree island area and tree planted diversity on biodiversity and ecosystem functioning operating through increasing tree dominance, through differences in structural complexity (indirect effects) or through alternative mechanisms (direct effects). As a hypothetical causal model, we included direct paths between island area and tree dominance (PC2) and between tree planted diversity and structural complexity gradient (from open to dense and structurally complex vegetation, PC1; Extended Data Fig. 2). Piecewise SEMs are based on a set of linear equations which are evaluated individually109. For our analyses, we included:
lm(restoration outcome ≈ island area + planted diversity) (1)
lm(structural complexity ≈ planted diversity) (2)
lm(tree dominance ≈ island area) (3)
Across all restoration outcomes, the main variables were always included in the linear model (1). As tree dominance and structural complexity are potential mechanistic pathways explaining the influence of island area and tree planted diversity, alternative paths between them and biodiversity or ecosystem functioning were added, if they improved the model fit (based on modification indices, P < 0.05). Therefore, model selection influenced only the inclusion of structural complexity and tree dominance in the linear model (1). Effects of island area and planted diversity through mechanistic pathways were calculated by multiplying their effect on the mechanistic explanatory variable and the effect of the mechanistic explanatory variable on biodiversity or ecosystem functioning. Mechanisms that were not captured by either of our proposed mechanistic pathways are represented by the direct paths between island area and tree planted diversity and biodiversity or ecosystem functioning. We tested the assumption of normality of the residuals in models (1), (2) and (3) using Shapiro–Wilk normality test, applied a suitable transformation of the response variables if needed (package bestNormalize v.1.6.1). The transformation concerned four out of ten indicators for species richness and Shannon diversity, three indicators out of ten indicators for Simpson diversity and 14 indicators out of 19 for ecosystem functioning. Effect sizes were calculated using standardized coefficients. The island area and planted tree diversity were log-transformed as this improved the model fit. For each SEM, we quantified the goodness of fit using the following metrics: Fisher’s C statistic and significance value based on a Chi-square test, the information criterion (Akaike information criterion (AIC), Bayesian information criterion (BIC), corrected AIC (AICc)) and pseudo-R2 values and applied the test of directed separation as implemented in the package piecewiseSEM v.2.1.0109.
Inclusion and ethics
The research included researchers from the Indonesian institutes Jambi University and Bogor Agricultural University throughout the research process—study design, study implementation, data ownership, intellectual property and authorship of publications. Local and regional research relevant to our study was considered in citations.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
www.nature.com
https://www.nature.com/articles/s41586-023-06086-5
