Tropical rainforests around the world are rapidly being converted into cash-crop agricultural systems. The associated massive losses of plant and animal species lead to changes in arthropod food webs and the energy fluxes therein. These changes are poorly understood, in particular in the extremely biodiverse canopies of tropical ecosystems. Using canopy fogging followed by stable isotope and energy flux analyses, we show that land-use conversion from rainforest to rubber and oil palm plantations not only causes a drastic reduction in energy fluxes of up to 75% but also shifts fluxes among trophic groups. While rainforests featured high levels of both herbivory and algae-microbiology, and a balanced ratio of herbivory to predation, relative fluxes were shifted towards predation in rubber and towards herbivory in oil palm plantations, indicating profound shifts in ecosystem functioning. Our results highlight that the ongoing loss of animal biodiversity and biomass in tropical canopies degrades animal-driven functions and restructures canopy food webs.

Study region and experimental design

Canopy arthropods were collected from research plots (n =32) of the EFForTS (Ecological and Socioeconomic Functions of Tropical Lowland Rainforest Transformation Systems) project in the lowlands of Jambi Province, Sumatra, Indonesia (Supplementary Fig. S1; Drescher et al. 2016), a region dominated by smallholder cash crop agriculture (Clough et al. 2016). Eight 50 m ⨯ 50 m plots each were established in three of the most dominant agricultural land-use systems in Jambi Province, i.e. monocultures of rubber (Hevea brasiliensis), oil palm (Elaeis guineensis), and jungle rubber (rubber agroforest system in gradually degrading rainforest remnants; Rembold et al. 2017), as well as in lowland rainforest as control habitat (Supplementary Fig. S2). The 32 research plots were located within and around two lowland rainforest reserves, i.e. the Bukit Duabelas National Park and the Harapan Rainforest, which are covered by degraded primary rainforest (Margono et al. 2014). The plots inside and near Bukit Duabelas National Park are referred to as belonging to the “Bukit Duabelas landscape”, while the research plots inside and near Harapan rainforest are affiliated to the “Harapan landscape” (Drescher et al. 2016).

Sampling and analyses of canopy arthropods and leaves

Arthropods were collected by canopy fogging, using the SwingFog® SN50 fogger (Swingtec GmbH, Germany) to apply a mixture of 50 mL DECIS 25 EC® (Bayer Crop Science, deltamethrine 25g/L) and four litres petroleum white oil to tree and palm canopies (Fig. 1). Three target canopies were fogged per 50 m ⨯ 50 m plot, located randomly but >20 m apart from each other; fallen trees and canopy gaps were avoided. Canopy fogging was conducted twice in the dry and rainy seasons of 2013 and 2013/14, respectively. Fogging was conducted in dry conditions and took about 20 min for each ~4L mix of while oil with insecticide. It began no later than 7:00 am for the first within-plot replicate and stopped before 9:00 am for the third replicate to avoid turbulences caused by rising daytime temperature above the canopy. Sixteen 1 m x 1 m fogging trays were arranged around a central tree or palm in a circle which was within 1 m from the trunk with its inner edge (Fig. S2 A, B). Each tray was fitted with a 250 mL plastic bottle containing 100 mL 96% EtOH (Fig. S2 C, D). Insecticide coverage was optimized by the person operating the fogger with the assistance from two observers on opposing sides of the target canopy. This ensured complete coverage of the canopy in rubber and oil palm; coverage in jungle rubber and rainforest was lower due to the insecticide plume dissolving with increased height. Typically, the fogging tray circles were underneath one oil palm, 2–3 in rubber trees and >5 trees in jungle rubber and rainforest. Two hours after fogging, stunned or dead arthropods were collected from the fogging trays, cleaned from debris, the EtOH exchanged and the samples stored at -20°C. Collected arthropods were sorted to orders, of which we used the 12 major ones in this study: Acarina, Araneae, Blattodea (without termites), Coleoptera, Collembola, Diptera, Hemiptera, Hymenoptera (Formicidae, Braconidae), Lepidoptera, Orthoptera, Psocoptera and Thysanoptera. Other, large-bodied arthropod taxa (e.g., Myriapoda, Mantodea, Phasmatodea, Scorpiones) were excluded from the analysis due to inefficiency of the insecticide (J. Drescher, pers. obs.). Termites were not included as their highly patchy distribution (0 – 66,084 individuals per plot across the two seasons) reflects that they were inadequately sampled (Jones & Eggleton 2000). Hymenoptera in this study are represented by Formicidae (ants) and Braconidae (a dominating family of parasitoid wasps; Azhar et al. 2022). Within Coleoptera, the four most abundant families Chrysomelidae, Curculionidae, Elateridae and Staphylinidae (Kasmiatun et al. 2022), which contract in feeding strategies, were separately analysed from ‘other Coleoptera’. Average body length and width were measured for 15 (or if <15 in the sample as many as present) randomly selected individuals from each order/family at plot-level in the dry and rainy season. Biomasses (fresh weight) per square metre were calculated using taxon-specific regressions (Supplementary Table S1). To calibrate stable isotope measurements of canopy arthropods to their potential basal resource, five randomly selected leaves per plot at head height were collected. In rainforest, leaves were from five tree species; in jungle rubber, two were from rubber trees and three from other tree species. In rubber, one leaf each was collected from the five rubber trees closest to the plot centre. In oil palm, we collected 40 cm² leaf material from five different oil palms closest to the plot centre.

Stable isotope ratios of carbon and nitrogen (¹³C/¹²C and ¹⁵N/¹⁴N), and C-to-N ratios were measured for pooled leaves from each plot (Supplementary Table S2) and for 15 (or if <15 in the sample as many as present) randomly selected and pooled individuals per plot for each order/family both in the dry and rainy season, resulting in 1073 measured samples of canopy arthropods and 32 measured samples of canopy leaves. Samples were measured with a coupled system of an elemental analyser and a mass spectrometer (Reineking et al. 1993). Isotopic signatures were expressed using the δ notation, with δX (‰) = (R_sample - R_standard)/R_standard x 1000; X represents the target isotope, and R_sample and R_standard represent ¹³C/¹²C and ¹⁵N/¹⁴N ratios of samples and standard, respectively. PD belemnite (PDB) and atmospheric nitrogen served as the primary standard for δ¹³C and δ¹⁵N, respectively. Acetanilide (C₈H₉NO, Merck, Darmstadt) was used for internal calibration.

Statistical analyses

All statistical analyses were conducted with R version 4.1.0 (2021-05-18) "Camp Pontanezen". To account for different isotopic baselines among different land-use types, bulk isotope values (δ¹³C and δ¹⁵N) of canopy arthropod orders were normalized to trees or oil palms representing primary producers by subtracting the isotopic signatures of leaves from the respective plots. These leaf-calibrated isotope data are denoted as Δ¹³C and Δ¹⁵N. Using these values, we calculated biomass-weighted isotopic metrics of canopy arthropod orders for each plot level community including mean, minimum, maximum, divergence, dispersion, evenness and uniqueness (Cucherousset and Villéger 2015). To estimate biomass distribution among trophic levels as a measure of trophic redundancy, we assigned ‘Δ¹⁵N classes’ from the highest to the lowest Δ¹⁵N values, each with a span of 3 ‰ (equivalent to approximately one trophic level) and summed the biomass of taxa in each class for each plot. To test for differences in abundance, biomass and isotopic composition between canopy arthropod orders, land-use systems, landscapes and seasons, and for differences in biomass distribution among trophic levels and land-use systems, landscapes and seasons, we constructed linear mixed effects models in R using the packages ‘lme4’ (Bates et al. 2015) and ‘lmerTest’ (Kuznetsova et al. 2017). Data distribution was inspected using the ‘fitdistrplus’ package (Delignette-Muller and Dutang 2015) and data were log-transformed if necessary. ‘Plot’ was included in the models as random effect. To test for differences among mean, minimum and maximum Δ¹³C and Δ¹⁵N values, we used the R packages ‘nlme’ (Pinheiro et al. 2021) and ‘mass’ (Ripley et al. 2019) to construct linear models with landscape, land-use and season as fixed factors.

For reconstructing the trophic structure of the studied taxa and calculating energy fluxes among them, we generated predator-prey adjacency matrices for each plot in both the dry and rainy season based on (1) bulk stable isotope composition, (2) optimum predator-prey mass-ratios (PPMR) and (3) biomass-based preferences (Potapov 2022). Bulk isotope composition was used to calculate ‘optimum’ prey or food resource for each animal group by taking into account a trophic enrichment of 2.3 ‰ for δ¹⁵N and of 1 ‰ for δ¹³C between prey and predator (Tiunov 2007). Taxa with δ¹⁵N values below those of plant leaves were assumed to mainly feed on algae/microbes that have lower δ¹⁵N values than plants (Potapov et al. 2019). PPMR was used as a characteristic that reflects size-based predation (i.e., small predator feeds on small prey and large predator can also feed on larger prey) and optimum foraging strategy (i.e., balancing the energetic profit and handling efforts; Brose et al. 2008), commonly used in food-web ecology (Brose et al. 2019). The optimum PPMR was set to 100, implying that typical prey has 100 times less mass than the predator (Brose et al. 2008). Since this value is derived from laboratory experiments and modelling, we allowed for a very broad range for “optimum” prey (PPMR width), i.e. body mass range of the optimum prey was set to be triple the body mass range of the predator, representing a large niche. Parasitoids (parasitoid wasps/Braconidae and Diptera) and ants were ‘allowed’ to feed on larger prey and the range for potential prey was set two times wider than for other groups due to parasitic lifestyle/pack hunting (Potapov 2022). Biomass-based preferences were set up assuming that prey preference scaled with available prey biomass (Gauzens et al. 2019). The three optimum prey (adjacency) matrices above were multiplied to obtain a final food web matrix for each plot in each season representing feeding preferences among food web nodes (taxa). The food web matrices were subsequently used to calculate energy fluxes per square metre at plot-level using the R package ‘Fluxweb’ (Gauzens et al. 2019). When applying the ‘fluxing’ function, biomass preferences were set to ‘false’ as they were already accounted for in the food web matrices. Biomass losses were set to ‘true’, as metabolic losses of taxa were defined per unit of biomass. Per capita metabolic rates in W based on metabolic theory scaling (Brown 2004) were calculated assuming a constant temperature of 25.2 °C and using general coefficients for invertebrates (Jochum et al. 2021). The efficiency level was set to ‘predator’, i.e. the efficiencies with which the predator/consumer assimilates consumed prey were used. Temperature-corrected assimilation efficiencies of food for predators (0.915) and herbivores (0.573) were calculated using parameters from Lang et al. (2017) and the mean annual temperature measured by meteo-stations across all studied plots, i.e. 25.2 °C (Drescher et al. 2016). We assumed assimilation efficiencies of algae-microbivores to be similar to herbivores. To infer ecological functions, the fluxes to herbivores were summed up as herbivory, the fluxes to algae-microbivores were summed up as algae-microbivory and the fluxes to predators were summed up as predation.

We additionally ran a sensitivity analysis to test if assumptions used for energy flux calculations will influence the conclusions of our study (Supplementary Text S1 and Fig. S4). To inspect differences between ecological functions depending on land-use, region and season, we again constructed linear mixed effects models with plot as random term (see above). 

All statistical analyses were conducted with R version 4.1.0 (2021-05-18) "Camp Pontanezen".