Philippa Kaur

Knowlton, N. The future of coral reefs. Proc. Natl. Acad. Sci. 98, 5419–5425 (2001).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Hughes, T. P. et al. Climate change, human impacts, and the resilience of coral reefs. Science 1979(301), 929–933 (2003).Article 
ADS 

Google Scholar 
Hoegh-Guldberg, O. et al. Coral reefs under rapid climate change and ocean acidification. Science 1979(318), 1737–1742 (2007).Article 
ADS 

Google Scholar 
Eakin, C. M. et al. Monitoring coral reefs from space. Oceanography 23, 118–133 (2010).Article 

Google Scholar 
Hughes, T. P. et al. Global warming transforms coral reef assemblages. Nature 556, 492–496 (2018).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Baker, D. M., Freeman, C. J., Wong, J. C. Y., Fogel, M. L. & Knowlton, N. Climate change promotes parasitism in a coral symbiosis. ISME J. 12, 921–930 (2018).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Hughes, T. P. et al. Global warming and recurrent mass bleaching of corals. Nature 543, 373–377 (2017).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Berkelmans, R. & van Oppen, M. J. H. The role of zooxanthellae in the thermal tolerance of corals: a ‘nugget of hope’ for coral reefs in an era of climate change. Proc. R. Soc. B Biol. Sci. 273, 2305–2312 (2006).Article 

Google Scholar 
Byler, K. A., Carmi-Veal, M., Fine, M. & Goulet, T. L. Multiple symbiont acquisition strategies as an adaptive mechanism in the coral Stylophora pistillata. PLoS One 8, e59596 (2013).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Cumbo, V., van Oppen, M. & Baird, A. Temperature and Symbiodinium physiology affect the establishment and development of symbiosis in corals. Mar. Ecol. Prog. Ser. 587, 117–127 (2018).Article 
ADS 
CAS 

Google Scholar 
Mundy, C. N. & Babcock, R. C. Role of light intensity and spectral quality in coral settlement: Implications for depth-dependent settlement?. J. Exp. Mar Biol. Ecol. 223, 235–255 (1998).Article 

Google Scholar 
Gleason, D. F., Edmunds, P. J. & Gates, R. D. Ultraviolet radiation effects on the behavior and recruitment of larvae from the reef coral Porites astreoides. Mar. Biol. 148, 503–512 (2006).Article 

Google Scholar 
Yusuf, S., Zamani, N. P., Jompa, J. & Junior, M. Z. Larvae of the coral Acropora tenuis (Dana 1846) settle under controlled light intensity. IOP Conf. Ser. Earth Environ. Sci. 253, 012023 (2019).Article 

Google Scholar 
Vermeij, M. J. A., Marhaver, K. L., Huijbers, C. M., Nagelkerken, I. & Simpson, S. D. Coral larvae move toward reef sounds. PLoS One 5, e10660 (2010).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Doropoulos, C. et al. Characterizing the ecological trade-offs throughout the early ontogeny of coral recruitment. Ecol. Monogr. 86, 20–44 (2016).Article 

Google Scholar 
Morse, D. E., Hooker, N., Morse, A. N. C. & Jensen, R. A. Control of larval metamorphosis and recruitment in sympatric agariciid corals. J. Exp. Mar. Biol. Ecol. 116, 193–217 (1988).Article 

Google Scholar 
Price, N. Habitat selection, facilitation, and biotic settlement cues affect distribution and performance of coral recruits in French Polynesia. Oecologia 163, 747–758 (2010).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Ritson-Williams, R., Arnold, S. N., Paul, V. J. & Steneck, R. S. Larval settlement preferences of Acropora palmata and Montastraea faveolata in response to diverse red algae. Coral Reefs 33, 59–66 (2014).Article 
ADS 

Google Scholar 
Negri, A., Webster, N., Hill, R. & Heyward, A. Metamorphosis of broadcast spawning corals in response to bacteria isolated from crustose algae. Mar. Ecol. Prog. Ser. 223, 121–131 (2001).Article 
ADS 

Google Scholar 
Webster, N. S. et al. Metamorphosis of a scleractinian coral in response to microbial biofilms. Appl. Environ. Microbiol. 70, 1213–1221 (2004).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Erwin, P. M., Song, B. & Szmant, A. M. Settlement behavior of Acropora palmata planulae: effects of biofilm age and crustose coralline algal cover. In Proceedings of 11th International Coral Reef Symposium 24, (2008).Siboni, N. et al. Crustose coralline algae that promote coral larval settlement harbor distinct surface bacterial communities. Coral Reefs 39, 1703–1713 (2020).Article 

Google Scholar 
Petersen, L.-E. et al. Mono- and multispecies biofilms from a crustose coralline alga induce settlement in the scleractinian coral Leptastrea purpurea. Coral Reefs 40, 381–394 (2021).Article 

Google Scholar 
Jorissen, H. et al. Coral larval settlement preferences linked to crustose coralline algae with distinct chemical and microbial signatures. Sci. Rep. 11, 14610 (2021).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Tebben, J. et al. Induction of larval metamorphosis of the coral Acropora millepora by tetrabromopyrrole isolated from a pseudoalteromonas bacterium. PLoS One 6, e19082 (2011).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Tebben, J. et al. Chemical mediation of coral larval settlement by crustose coralline algae. Sci. Rep. 5, 10803 (2015).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Tran, C. & Hadfield, M. Larvae of Pocillopora damicornis (Anthozoa) settle and metamorphose in response to surface-biofilm bacteria. Mar. Ecol. Prog. Ser. 433, 85–96 (2011).Article 
ADS 

Google Scholar 
Sneed, J. M., Sharp, K. H., Ritchie, K. B. & Paul, V. J. The chemical cue tetrabromopyrrole from a biofilm bacterium induces settlement of multiple Caribbean corals. Proc. R. Soc. B Biol. Sci. 281, 20133086 (2014).Article 

Google Scholar 
Petersen, L.-E., Kellermann, M. Y., Nietzer, S. & Schupp, P. J. Photosensitivity of the bacterial pigment cycloprodigiosin enables settlement in coral larvae—light as an understudied environmental factor. Front. Mar. Sci. 8, 749070 (2021).Article 

Google Scholar 
Heyward, A. J. & Negri, A. P. Natural inducers for coral larval metamorphosis. Coral Reefs 18, 273–279 (1999).Article 

Google Scholar 
Harrington, L., Fabricius, K., Death, G. & Negri, A. Recognition and selection of settlement substrata determine post-settlement survival in corals. Ecology 85, 3428–3437 (2004).Article 

Google Scholar 
Da-Anoy, J. P., Villanueva, R. D., Cabaitan, P. C. & Conaco, C. Effects of coral extracts on survivorship, swimming behavior, and settlement of Pocillopora damicornis larvae. J. Exp. Mar. Biol. Ecol. 486, 93–97 (2017).Article 

Google Scholar 
Morse, D. E. & Morse, A. N. C. Enzymatic characterization of the morphogen recognized by Agaricia humilis (Scleractinian Coral) larvae. Biol. Bull. 181, 104–122 (1991).Article 
CAS 
PubMed 

Google Scholar 
Kitamura, M., Koyama, T., Nakano, Y. & Uemura, D. Characterization of a natural inducer of coral larval metamorphosis. J. Exp. Mar. Biol. Ecol. 340, 96–102 (2007).Article 

Google Scholar 
Kitamura, M., Schupp, P. J., Nakano, Y. & Uemura, D. Luminaolide, a novel metamorphosis-enhancing macrodiolide for scleractinian coral larvae from crustose coralline algae. Tetrahedron Lett. 50, 6606–6609 (2009).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Maru, N. et al. Relative configuration of luminaolide. Tetrahedron Lett. 54, 4385–4387 (2013).Article 
CAS 

Google Scholar 
Nietzer, S., Moeller, M., Kitamura, M. & Schupp, P. J. Coral larvae every day: Leptastrea purpurea, a brooding species that could accelerate coral research. Front. Mar. Sci. 5, 466 (2018).Article 

Google Scholar 
Moeller, M., Nietzer, S. & Schupp, P. J. Neuroactive compounds induce larval settlement in the scleractinian coral Leptastrea purpurea. Sci. Rep. 9, 2291 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Petersen, L.-E., Kellermann, M. Y. & Schupp, P. J. Secondary metabolites of marine microbes: from natural products chemistry to chemical ecology. In YOUMARES 9 – The Oceans: Our Research, Our Future: Proceedings of the 2018 Conference for Young Marine Researcher in Oldenburg, Germany (eds Jungblut, S. et al.) 159–180 (Springer International Publishing, 2020). https://doi.org/10.1007/978-3-030-20389-4_8.Chapter 

Google Scholar 
Fiegel, L. J. et al. A detailed visualization of the early development stages of Leptastrea purpurea reveals distinct bio-optical features. Front. Mar. Sci. 10, 1–10 (2023).
Google Scholar 
Strader, M. E., Aglyamova, G. V. & Matz, M. V. Molecular characterization of larval development from fertilization to metamorphosis in a reef-building coral. BMC Genom. 19, 17 (2018).Article 

Google Scholar 
Puisay, A. et al. Parental bleaching susceptibility leads to differences in larval fluorescence and dispersal potential in Pocillopora acuta corals. Mar. Environ. Res. 163, 105200 (2021).Article 
CAS 
PubMed 

Google Scholar 
Perez-Tomas, R. & Vinas, M. New insights on the antitumoral properties of prodiginines. Curr. Med. Chem. 17, 2222–2231 (2010).Article 
CAS 
PubMed 

Google Scholar 
You, Z. et al. Insights into the anti-infective properties of prodiginines. Appl. Microbiol. Biotechnol. 103, 2873–2887 (2019).Article 
CAS 
PubMed 

Google Scholar 
Kellermann, M. Y., Yoshinaga, M. Y., Valentine, R. C., Wörmer, L. & Valentine, D. L. Important roles for membrane lipids in haloarchaeal bioenergetics. Biochim. Biophys. Acta (BBA) Biomembr. 1858, 2940–2956 (2016).Article 
CAS 

Google Scholar 
Hirose, M., Yamamoto, H. & Nonaka, M. Metamorphosis and acquisition of symbiotic algae in planula larvae and primary polyps of Acropora spp.. Coral Reefs 27, 247–254 (2008).Article 
ADS 

Google Scholar 
Bollati, E. et al. Green fluorescent protein-like pigments optimize the internal light environment in symbiotic reef building corals. Elife 11, e73521 (2022).Article 
PubMed 
PubMed Central 

Google Scholar 
Palmer, C. V., Modi, C. K. & Mydlarz, L. D. Coral fluorescent proteins as antioxidants. PLoS One 4, e7298 (2009).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Alegado, R. A. et al. A bacterial sulfonolipid triggers multicellular development in the closest living relatives of animals. Elife 1, e00013 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
Woznica, A. et al. Bacterial lipids activate, synergize, and inhibit a developmental switch in choanoflagellates. Proc. Natl. Acad. Sci. 113, 7894–7899 (2016).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
He, J. et al. Bacterial nucleobases synergistically induce larval settlement and metamorphosis in the invasive mussel Mytilopsis sallei. Appl. Environ. Microbiol. 85, e01039 (2019).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Guo, H., Rischer, M., Westermann, M. & Beemelmanns, C. Two distinct bacterial biofilm components trigger metamorphosis in the colonial hydrozoan Hydractinia echinata. MBio 12, e00401 (2021).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Ross, C., Fogarty, N. D., Ritson-Williams, R. & Paul, V. J. Interspecific variation in coral settlement and fertilization success in response to hydrogen peroxide exposure. Biol. Bull. 233, 206–218 (2017).Article 
CAS 
PubMed 

Google Scholar 
Boettcher, A. A., Dyer, C., Casey, J. & Targett, N. M. Hydrogen peroxide induced metamorphosis of queen conch, Strombus gigas: tests at the commercial scale. Aquaculture 148, 247–258 (1997).Article 
CAS 

Google Scholar 
Covarrubias, L., Hernández-García, D., Schnabel, D., Salas-Vidal, E. & Castro-Obregón, S. Function of reactive oxygen species during animal development: Passive or active?. Dev. Biol. 320, 1–11 (2008).Article 
CAS 
PubMed 

Google Scholar 
Gauron, C. et al. Hydrogen peroxide (H2O2) controls axon pathfinding during zebrafish development. Dev. Biol. 414, 133–141 (2016).Article 
CAS 
PubMed 

Google Scholar  More

EthicsNo ethics approval was required for this study. Our study was conducted in a state-owned industrial oil palm plantation where we established a cooperation with the estate owner to access the site and collect data. No endangered or protected species were sampled. Research permits were obtained from the Ministry of Research, Technology and Higher Education, and sample collection and sample export permits were obtained from the Ministry of Environment and Forestry of the Republic of Indonesia.Study area and experimental designOur study was conducted in a state-owned industrial oil palm plantation (PTPN VI) located in Jambi, Indonesia (1.719° S, 103.398° E, 73 m above sea level). Initial planting of oil palms within the 2,025 ha plantation area started in 1998 and ended in 2002; planting density was 142 palms ha−1, spaced 8 m apart in each row and between rows, and palms were ≥16 years old during our study period of 2016–2020. The study sites have a mean annual temperature of 27.0 ± 0.2 °C and a mean annual precipitation of 2,103 ± 445 mm (2008–2017, Sultan Thaha Airport, Jambi). The management practices in large-scale oil palm plantations typically result in three contrasting management zones: (1) a 2 m radius around the base of the palm that was weeded (four times a year) and raked before fertilizer application, hereafter called the ‘palm circle’; (2) an area occurring every second inter-row, where pruned senesced palm fronds were piled up, hereafter called ‘frond piles’; and (3) the remaining area of the plantation where less weeding (two times a year) and no fertilizer were applied, hereafter called ‘inter-rows’.Within this oil palm plantation, we established a management experiment in November 2016 with full factorial treatments of two fertilization rates × two weeding practices: conventional fertilization rates at PTPN VI and other large-scale plantations (260 kg N–50 kg P–220 kg K ha−1 yr−1), reduced fertilization rates based on quantified nutrient export by harvest (136 kg N–17 kg P ha−1 yr−1–187 kg K ha−1 yr−1), herbicide and mechanical weeding15. The reduced fertilization treatment was based on quantified nutrient export from fruit harvest, calculated by multiplying the nutrient content of fruit bunches with the long-term yield data of the plantation. Fertilizers were applied yearly in April and October following weeding and raking of the palm circle. The common practice at PTPN VI and other large-scale plantations on acidic Acrisol soils is to apply lime and micronutrients, and these were unchanged in our management experiment. Before each N–P–K fertilizer application, dolomite and micronutrients were applied to the palm circle in all treatment plots using the common rates)52: 426 kg ha−1 yr−1 dolomite and 142 kg Micro-Mag ha−1 yr−1 (containing 0.5% B2O3, 0.5% CuO, 0.25% Fe2O3, 0.15% ZnO, 0.1% MnO and 18% MgO). Herbicide treatment was carried out using glyphosate in the palm circle (1.50 l ha−1 yr−1, split into four applications per year) and in the inter-rows (0.75 l ha−1 yr−1, split into two applications per year). Mechanical weeding was done using a brush cutter in the same management zones and at the same frequency as the herbicide treatment.The 22 factorial design resulted in four treatment combinations: conventional fertilization with herbicide treatment, reduced fertilization with herbicide treatment, conventional fertilization with mechanical weeding and reduced fertilization with mechanical weeding. The four treatments were randomly assigned on 50 m × 50 m plots replicated in four blocks, totalling 16 plots. The effective measurement area was the inner 30 m × 30 m area within each replicate plot to avoid any possible edge effects. For indicators (below) that were measured within subplots, these subplots were distributed randomly within the inner 30 m × 30 m of a plot. All replicate plots were located on flat terrain and on an Acrisol soil with a sandy clay loam texture.Ecosystem functions and multifunctionalityOur study included multiple indicators for each of the eight ecosystem functions23, described in details below (Supplementary Tables 1 and 2). All the parameters were expressed at the plot level by taking the means of the subplots (that is, biological parameters) or the area-weighted average of the three management zones per plot (that is, soil parameters). (1) Greenhouse gas (GHG) regulation was indicated by NEP, soil organic C (SOC) and soil GHG fluxes. (2) Erosion prevention was signified by the understory vegetation cover during the four-year measurements. (3) Organic matter decomposition was indicated by leaf litter decomposition and soil animal decomposer activity. (4) Soil fertility was signified by gross N mineralization rate, effective cation exchange capacity (ECEC), base saturation and microbial biomass N. (5) Pollination potential was designated by pan-trapped arthropod abundance and nectar-feeding bird activity. As such, it does not quantify the pollination potential for oil palm, which is mainly pollinated by a single weevil species, but rather as a proxy for a general pollination potential for other co-occurring plants. (6) Water filtration (the capacity to provide clean water) was indicated by leaching losses of the major elements. (7) Plant refugium (the capacity to provide a suitable habitat for plants) as signified by the percentage ground cover of invasive plants to the total ground cover of understory vegetation during the four-year measurements. (8) Biological control (the regulation of herbivores via predation) was indicated by insectivorous bird and bat activities and the soil arthropod predator activity.All the ecosystem functions were merged into a multifunctionality index using the established average and threshold approaches12. For average multifunctionality, we first averaged the z-standardized values (Statistics) of indicators for each ecosystem function and calculated the mean of the eight ecosystem functions for each plot. For threshold multifunctionality, this was calculated from the number of functions that exceeds a set threshold, which is a percentage of the maximum performance level of each function12; we investigated the range of thresholds from 10% to 90% to have a complete overview. The maximum performance was taken as the average of the three highest values for each indicator per ecosystem function across all plots to reduce effect of potential outliers. For each plot, we counted the number of indicators that exceeded a given threshold for each function and divided by the number of indicators for each function12.Indicators of GHG regulationWe calculated annual NEP for each plot as: net ecosystem C exchange – harvested fruit biomass C (ref. 16), whereby net ecosystem C exchange = Cout (or heterotrophic respiration) – Cin (or net primary productivity)53. The net primary productivity of oil palms in each plot was the sum of aboveground biomass production (aboveground biomass C + frond litter biomass C input + fruit biomass C) and belowground biomass production. Aboveground biomass production was estimated using allometric equations developed for oil palm plantations in Indonesia54, using the height of palms measured yearly from 2019 to 2020. Annual frond litter biomass input was calculated from the number and dry mass of fronds pruned during harvesting events of an entire year in each plot and was averaged for 2019 and 2020. Aboveground biomass production was converted to C based on C concentrations in wood and leaf litter55. Annual fruit biomass C production (which is also the harvest export) was calculated from the average annual yield in 2019 and 2020 and the measured C concentrations of fruit bunches. Belowground root biomass and litter C production were taken from previous work in our study area55, and it was assumed constant for each plot. Heterotrophic respiration was estimated for each plot as: annual soil CO2 C emission (below) × 0.7 (based on 30% root respiration contribution to soil respiration from a tropical forest in Sulawesi, Indonesia56) + annual frond litter biomass C input × 0.8 (~80% of frond litter is decomposed within a year in this oil palm plantation8). SOC was measured in March 2018 from composite samples collected from two subplots in each of the three management zones per plot down to 50 cm depth. Soil samples were air dried, finely ground and analysed for SOC using a CN analyser (Vario EL Cube, Elementar Analysis Systems). SOC stocks were calculated using the measured bulk density in each management zone, and values for each plot were the area-weighted average of the three management zones (18% for palm circle, 15% for frond piles and 67% for inter-rows)15,22.From July 2019 to June 2020, we conducted monthly measurements of soil CO2, CH4 and N2O fluxes using vented, static chambers permanently installed in the three management zones within two subplots per plot11,57. Annual soil CO2, CH4 and N2O fluxes were trapezoidal interpolations between measurement periods for the whole year, and values for each plot were the area-weighted average of the three management zones (above).Indicators of erosion preventionDiversity and abundance of vascular plants were assessed once a year from 2016 to 2020 before weeding in September–November. In five subplots per plot, we recorded the occurrence of all vascular plant species and estimated the percent cover of the understory vegetation. The percentage cover and plant species richness of each measurement year were expressed in ratio to that of 2016 to account for initial differences among the plots before the start of the experiment. For example, percentage cover in 2017 was:$$mathrm{Cover}_{2017} = frac{{left( {mathrm{Cover}_{2017} – mathrm{Cover}_{2016}} right)}}{{mathrm{Cover}_{2016}}}$$The values from five subplots were averaged to represent each plot.Indicators of organic matter decompositionLeaf litter decomposition was determined using litter bags (20 cm × 20 cm with 4 mm mesh size) containing 10 g of dry oil palm leaf litter8. Three litter bags per plot were placed on the edge of the frond piles in December 2016. After eight months of incubation in the field, we calculated leaf litter decomposition as the difference between initial litter dry mass and litter dry mass following incubation. Soil animal decomposer activity is described below (Soil arthropods).Indicators of soil fertilityAll these indicators were measured in February–March 2018 in the three management zones within two subplots per plot22. Gross N mineralization rate in the soil was measured in the top 5 cm depth on intact soil cores incubated in situ using the 15N pool dilution technique58. ECEC and base saturation were measured in the top 5 cm depth as this is the depth that reacts fast to changes in management22. The exchangeable cation concentrations (Ca, Mg, K, Na, Al, Fe, Mn) were determined by percolating the soil with 1 mol l−1 of unbuffered NH4Cl, followed by analysis of the percolates using an inductively coupled plasma-atomic emission spectrometer (ICP-AES; iCAP 6300 Duo view ICP Spectrometer, Thermo Fisher Scientific). Base saturation was calculated as the percentage exchangeable bases (Mg, Ca, K and Na) on ECEC. Microbial biomass N was measured from fresh soil samples using the fumigation-extraction method59. The values for each plot were the mean of the two subplots that were the area-weighted average of the three management zones (above)15,22.Indicators of general pollination potentialFluorescent yellow pan traps were used to sample aboveground arthropods (to determine pollinator communities60) in November 2016, September 2017 and June 2018. The traps were attached to a platform at the height of the surrounding vegetation within a 2 × 3 grid centred in the inter-rows of each plot in six clusters of three traps, totalling 18 traps per plot. Traps were exposed in the field for 48 h. We stored all trapped arthropods in 70% ethanol and later counted and identified to order and species level. The abundance of trapped arthropods in 2017 and 2018 were calculated as the ratio to the abundance in 2016 to account for initial differences among the plots before the start of the experiment. The activity of nectar-feeding birds is described below (Birds and bats).Indicators of water filtrationElement leaching losses were determined from analyses of soil-pore water sampled monthly at 1.5 m depth using suction cup lysimeters (P80 ceramic, maximum pore size 1 μm; CeramTec) over the course of one year (2017–2018)15. Lysimeters were installed in the three management zones within two subplots per plot. Dissolved N was analysed using continuous flow injection colorimetry (SEAL Analytical AA3, SEAL Analytical), whereas these other elements were determined using ICP-AES. The values for each plot were the mean of the two subplots that were the area-weighted average of the three management zones15,22.Indicators of plant refugiumIn five subplots per plot, the percentage cover and species richness of invasive understory plant species were assessed once a year from 2016 to 2020 before weeding in September–November. We defined invasive species as those plants non-native to Sumatra61 and among the ten dominant species (excluding oil palm) in the plantation for each year. The percentage cover of invasive understory plant species of each measurement year was expressed in a ratio to that of 2016 to account for initial differences among the plots before the start of the experiment. The values for each plot were represented by the average of five subplots.Indicators of biological controlThe activities of insectivorous birds and bats are described below (Birds and bats). In five subplots per plot, soil invertebrates were collected (Soil arthropods), counted, identified to taxonomic order level and subsequently classified according to their trophic groups that include predators60. The values from five subplots were average to represent each plot.BiodiversityBiodiversity was measured by the taxonomic richness of seven multitrophic groups, described in details below (Supplementary Tables 1 and 2).Understory plant species richnessThe method is described above (Indicators of erosion prevention), using the number of species as an indicator (Supplementary Table 2).Soil microorganism richnessThis was determined in May 2017 by co-extracting RNA and DNA from three soil cores (5 cm diameter, 7 cm depth) in five subplots per plot62. While DNA extraction describes the entire microbial community, RNA represents the active community. The v3–v4 region of the 16S rRNA gene was amplified and sequenced with a MiSeq sequencer (Illumina). Taxonomic classification was done by mapping curated sequences against the SILVA small subunit (SSU) 138 non-redundant (NR) database63 with the Basic Local Alignment Search Tool (BLASTN)64.Soil arthropod order richnessFor determination of soil arthropods, we collected soil samples (16 cm × 16 cm, 5 cm depth) in five subplots per plot in October–November 2017. We extracted the animals from the soil using a heat-gradient extractor65, collected them in dimethyleneglycol-water solution (1:1) and stored in 80% ethanol. The extracted animals were counted and identified to taxonomic order level61. They were also assigned to the trophic groups decomposers, herbivores and predators based on the predominant food resources recorded in previous reviews and a local study66,67. Orders with diverse feeding habits were divided into several feeding groups, for example, Coleoptera were divided into mostly predatory families (Staphylinidae, Carabidae), herbivorous families (for example, Curculionidae) and decomposer families (for example, Tenebrionidae). The total number of individuals per taxonomic group in each subplot was multiplied by the group-specific metabolic rate, which were summed to calculate soil animal decomposer activity. The values from five subplots were average to represent each plot.Aboveground arthropod order and insect family richnessIn addition to the fluorescent yellow pan traps described above (Indicators of general pollination potential), sweep net and Malaise trap samplings were conducted in June 2018, which targeted the general flying and understory dwelling arthropod communities. Sweep net sampling was conducted within the understory vegetation along two 10 m long transects per plot, with ten sweeping strokes performed per transect. In each plot, we installed a single Malaise trap between two randomly chosen palms and exposed it for 24 h. Arthropods were counted, identified to taxonomic order level and the insects to taxonomic family level and values from the three methods were summed to represent each plot.Birds and batsBirds and bats passing at each replicate plot were sampled in September 2017 using SM2Bat + sound recorders (Wildlife Acoustics) with two microphones (SMX-II and SMX-US) placed at a height of 1.5 m in the middle of each plot68. We assigned the bird vocalization to species with Xeno-Canto69 and the Macaulay library70. Insectivorous bat species richness was computed by dividing them into morphospecies based on the characteristics of their call (call frequency, duration, shape). In addition, we gathered information on proportional diet preferences of the bird species using the EltonTrait database71. We defined birds feeding on invertebrates (potential biocontrol agents) as the species with a diet of at least 80% invertebrates and feeding on nectar (potential pollinators), if the diet included at least 20% of nectar.Economic indicatorsWe used six indicators linked to the level and stability of yield and profit: yield, lower fifth quantile of the yield per palm per plot, shortfall probability, management costs, profit and relative gross margin. We assessed fruit yield by weighing the harvested fruit bunches from each palm within the inner 30 m × 30 m area of each plot. The harvest followed the schedule and standard practices of the plantation company: each palm was harvested approximately every ten days and the lower fronds were pruned. For each plot, we calculated the average fruit yield per palm and scaled up to a hectare, considering the planting density of 142 palms per ha. Because the palms in each plot have different fruiting cycles and were harvested continuously, the calculation of an annual yield may lead to misleading differences between treatments. Therefore, we calculated the cumulative yield from the beginning of the experiment to four years (2017–2020), which should account for the inter- and intra-annual variations in fruit production of the palms in the plots and thus allowing for comparison among treatments. As effects of management practices on yield may be delayed46, we also calculated the cumulative yield during two consecutive years (2017–2018 and 2019–2020) and checked for treatment effects on yield and profit indicators separately for these two periods.We computed risk indicators on the cumulative yield and on the yield between the two periods. We used the lowest fifth quantile of the yield per palm per plot (left side of the distribution) to indicate the production of the palms with lowest performance. Also, we determined the yield shortfall probability (lower partial moment 0th order), defined as the share of palms that fell below a predefined threshold of yield; the thresholds chosen were 630 kg−1 per palm for cumulative yield and 300 kg−1 per palm per year for the two-year yield, which corresponded to 75% of the average yield.Revenues and costs were calculated as cumulative values during four years of the experiment (2017–2020) using the same prices and costs for all the years. This was because we were interested in assessing the economic consequences of different management treatments, and they might be difficult to interpret when changes in prices and costs between calendar years are included, which are driven by external market powers rather than the field-management practices. For the same reason, we abstained from discounting profits. Given the usually high discount rates applied to the study area, slight differences in harvesting activities between calendar years or months might lead to high systematic differences between the management treatments, which are associated with the variation in work schedule within the plantation rather than the actual difference among management treatments. Revenues were calculated from the yield and the average price of the fruit bunches in 2016 and 201761. Material costs were the sum of the costs of fertilizers, herbicide and gasoline for the brush cutter. Labour costs were calculated from the minimum wage in Jambi and the time (in labour hours) needed for the harvesting, fertilizing and weeding operations, which were recorded in 2017 for each plot. The weeding labour included the labour for raking the palm circle before fertilization, which was equal in all treatments, and the weeding in the palm circle and inter-rows either with herbicide or brush cutter. In addition, we included the time to remove C. hirta, which must be removed mechanically from all plots once a year, calculated from the average weed-removal time in the palm circle and the percentage cover of C. hirta in each plot for each year. We then calculated the profit as the difference between revenues and the total management costs and the relative gross margin as the gross profit proportion of the revenues.StatisticsTo test for differences among management treatments for each ecosystem function and across indicators of biodiversity, the plot-level value of each indicator was first z standardized (z = (actual value − mean value across plots) / standard deviation)4. This prevents the dominance of one or few indicators over the others, and z standardization allows several distinct indicators to best characterize an ecosystem function or biodiversity4. Standardized values were inverted (multiplied by −1) for indicators of which high values signify undesirable effect (that is, NEP, soil N2O and CH4 fluxes, element leaching losses, invasive plant cover, yield shortfall, management costs) for intuitive interpretations. For a specific ecosystem function (Supplementary Figs. 1 and 2) and across indicators of biodiversity (Fig. 2), linear mixed-effects (LME) models were used to assess differences among management treatments (fertilization, weeding and their interaction) as fixed effects with replicate plots and indicators (Supplementary Tables 2 and 3) as random effects. The significance of the fixed effects was evaluated using ANOVA72. The LME model performance was assessed using diagnostic residual plots73. As indicator variables may systematically differ in their responses to management treatments, we also tested the interaction between indicator and treatment (Table 1). For testing the differences among management treatments across ecosystem functions (that is, multifunctionality; Fig. 1), we used for each replicate plot the average of z-standardized indicators of each ecosystem function and ranges of thresholds (that is, number of functions that exceeds a set percentage of the maximum performance of each function12; Supplementary Fig. 3). The LME models had management treatments (fertilization, weeding and their interaction) as fixed effects and replicate plots and ecosystem functions as random effects; the interaction between ecosystem function and treatment were also tested to assess if there were systematic differences in their responses to management treatments (Table 1). As we expected that the type of weeding will influence ground vegetation, we tested for differences in ground cover of understory vegetation, measured from 2016 to 2020, using LME with management treatments as fixed effect and replicate plots and year as random effects. Differences among management treatments (fertilization, weeding and their interaction) in yield and profit indicators, which were cumulative values over four years (Fig. 3) or for two separate periods (2017–2018 and 2019–2020; Supplementary Fig. 4), were assessed using linear model ANOVA (Table 1). For clear visual comparison among management treatments across ecosystem functions, multitrophic groups for biodiversity, and yield and profit indicators, the fifth and 95th percentiles of their z-standardized values were presented in a petal diagram (Fig. 4 and Supplementary Fig. 5). Data were analysed using R (version 4.0.4), using the R packages ‘nlme’ and ‘influence.ME’73.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More

Ritz, K. & Young, I. M. Interactions between soil structure and fungi. Mycologist 18, 52–59 (2004).Article 

Google Scholar 
Schimel, J. P. & Schaeffer, S. M. Microbial control over carbon cycling in soil. Front. Microbiol. 3, 348 (2012).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Six, J., Bossuyt, H., Degryze, S. & Denef, K. A history of research on the link between (micro)aggregates, soil biota, and soil organic matter dynamics. Soil Tillage Res. 79, 7–31 (2004).Article 

Google Scholar 
van der Heijden, M. G. A. & Wagg, C. Soil microbial diversity and agro-ecosystem functioning. Plant Soil 363, 1–5 (2013).Article 
CAS 

Google Scholar 
Winter, S. et al. Effects of vegetation management intensity on biodiversity and ecosystem services in vineyards: a meta-analysis. J. Appl. Ecol. 55, 2484–2495 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Belmonte, S. A. et al. Effect of long-term soil management on the mutual interaction among soil organic matter, microbial activity and aggregate stability in a vineyard. Pedosphere 28, 288–298 (2018).Article 
CAS 

Google Scholar 
Bronick, C. J. & Lal, R. Soil structure and management: a review. Geoderma 124, 3–22 (2005).Article 
ADS 
CAS 

Google Scholar 
Kratschmer, S. et al. Enhancing flowering plant functional richness improves wild bee diversity in vineyard inter-rows in different floral kingdoms. Ecol. Evol. 11, 7927–7945 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Constancias, F. et al. Microscale evidence for a high decrease of soil bacterial density and diversity by cropping. Agron. Sustain. Dev. 34, 831–840 (2014).Article 
CAS 

Google Scholar 
Schmidt, R., Gravuer, K., Bossange, A. V., Mitchell, J. & Scow, K. Long-term use of cover crops and no-till shift soil microbial community life strategies in agricultural soil. PLoS One 13, e0192953 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Vink, S. N., Chrysargyris, A., Tzortzakis, N. & Salles, J. F. Bacterial community dynamics varies with soil management and irrigation practices in grapevines (Vitis vinifera L.). Appl. Soil Ecol. 158, 103807 (2021).Article 

Google Scholar 
Pingel, M., Reineke, A. & Leyer, I. A 30-years vineyard trial: plant communities, soil microbial communities and litter decomposition respond more to soil treatment than to N fertilization. Agr. Ecosyst. Environ. 272, 114–125 (2019).Article 
CAS 

Google Scholar 
Sharma-Poudyal, D., Schlatter, D., Yin, C., Hulbert, S. & Paulitz, T. Long-term no-till: a major driver of fungal communities in dryland wheat cropping systems. PLoS One 12, e0184611 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Hungria, M., Franchini, J. C., Brandão-Junior, O., Kaschuk, G. & Souza, R. A. Soil microbial activity and crop sustainability in a long-term experiment with three soil-tillage and two crop-rotation systems. Appl. Soil. Ecol. 42, 288–296 (2009).Article 

Google Scholar 
Pascault, N. et al. In situ dynamics of microbial communities during decomposition of wheat, rape, and alfalfa residues. Microb. Ecol. 60, 816–828 (2010).Article 
PubMed 

Google Scholar 
Tresch, S. et al. Litter decomposition driven by soil fauna, plant diversity and soil management in urban gardens. Sci. Total Environ. 658, 1614–1629 (2019).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Faust, S., Koch, H.-J., Dyckmans, J. & Joergensen, R. G. Response of maize leaf decomposition in litterbags and soil bags to different tillage intensities in a long-term field trial. Appl. Soil. Ecol. 141, 38–44 (2019).Article 

Google Scholar 
Liu, Y.-R. et al. New insights into the role of microbial community composition in driving soil respiration rates. Soil Biol. Biochem. 118, 35–41 (2018).Article 
CAS 

Google Scholar 
Yang, C., Liu, N. & Zhang, Y. Soil aggregates regulate the impact of soil bacterial and fungal communities on soil respiration. Geoderma 337, 444–452 (2019).Article 
ADS 
CAS 

Google Scholar 
Delgado-Baquerizo, M. et al. Microbial diversity drives multifunctionality in terrestrial ecosystems. Nat. Commun. 7, 10541 (2016).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Bruggisser, O. T., Schmidt-Entling, M. H. & Bacher, S. Effects of vineyard management on biodiversity at three trophic levels. Biol. Cons. 143, 1521–1528 (2010).Article 

Google Scholar 
Lienhard, P. et al. Pyrosequencing evidences the impact of cropping on soil bacterial and fungal diversity in Laos tropical grassland. Agron. Sustain. Dev. 34, 525–533 (2014).Article 

Google Scholar 
Schnoor, T. K., Lekberg, Y., Rosendahl, S. & Olsson, P. A. Mechanical soil disturbance as a determinant of arbuscular mycorrhizal fungal communities in semi-natural grassland. Mycorrhiza 21, 211–220 (2011).Article 
PubMed 

Google Scholar 
Kazakou, E. et al. A plant trait-based response-and-effect framework to assess vineyard inter-row soil management. Bot. Lett. 163, 373–388 (2016).Article 

Google Scholar 
Svensson, J. R., Lindegarth, M., Jonsson, P. R. & Pavia, H. Disturbance-diversity models: What do they really predict and how are they tested?. Proc. Biol. Sci. 279, 2163–2170 (2012).PubMed 
PubMed Central 

Google Scholar 
Bao, T. et al. Moderate disturbance increases the PLFA diversity and biomass of the microbial community in biocrusts in the Loess Plateau region of China. Plant Soil 451, 499–513 (2020).Article 
CAS 

Google Scholar 
Liu, J. et al. Soil carbon content drives the biogeographical distribution of fungal communities in the black soil zone of northeast China. Soil Biol. Biochem. 83, 29–39 (2015).Article 
CAS 

Google Scholar 
Cotton, J. & Acosta-Martínez, V. Intensive tillage converting grassland to cropland immediately reduces soil microbial community size and organic carbon. Agric. Environ. Lett. 3, 180047 (2018).Article 

Google Scholar 
Poeplau, C. et al. Temporal dynamics of soil organic carbon after land-use change in the temperate zone – carbon response functions as a model approach. Glob. Change Biol. 17, 2415–2427 (2011).Article 
ADS 

Google Scholar 
Burns, K. N. et al. Vineyard soil bacterial diversity and composition revealed by 16S rRNA genes: differentiation by vineyard management. Soil Biol. Biochem. 103, 337–348 (2016).Article 
CAS 

Google Scholar 
Steiner, M. et al. Local conditions matter: minimal and variable effects of soil disturbance on microbial communities and functions in European vineyards. PLoS One 18, e0280516 (2023).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Zeng, J. et al. Nitrogen fertilization directly affects soil bacterial diversity and indirectly affects bacterial community composition. Soil Biol. Biochem. 92, 41–49 (2016).Article 
CAS 

Google Scholar 
Fierer, N. & Jackson, R. B. The diversity and biogeography of soil bacterial communities. Proc. Natl. Acad. Sci. U.S.A. 103, 626–631 (2006).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Eisenhauer, N. Plant diversity effects on soil microorganisms: spatial and temporal heterogeneity of plant inputs increase soil biodiversity. Pedobiologia 59, 175–177 (2016).Article 

Google Scholar 
Porazinska, D. L. et al. Plant diversity and density predict belowground diversity and function in an early successional alpine ecosystem. Ecology 99, 1942–1952 (2018).Article 
PubMed 

Google Scholar 
Prober, S. M. et al. Plant diversity predicts beta but not alpha diversity of soil microbes across grasslands worldwide. Ecol. Lett. 18, 85–95 (2015).Article 
PubMed 

Google Scholar 
Sun, Y.-Q., Wang, J., Shen, C., He, J.-Z. & Ge, Y. Plant evenness modulates the effect of plant richness on soil bacterial diversity. Sci. Total Environ. 662, 8–14 (2019).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Kuzyakov, Y. Priming effects: interactions between living and dead organic matter. Soil Biol. Biochem. 42, 1363–1371 (2010).Article 
CAS 

Google Scholar 
Huo, C., Luo, Y. & Cheng, W. Rhizosphere priming effect: a meta-analysis. Soil Biol. Biochem. 111, 78–84 (2017).Article 
CAS 

Google Scholar 
Dimassi, B. et al. Effect of nutrients availability and long-term tillage on priming effect and soil C mineralization. Soil Biol. Biochem. 78, 332–339 (2014).Article 
CAS 

Google Scholar 
Prescott, C. E. Litter decomposition: What controls it and how can we alter it to sequester more carbon in forest soils?. Biogeochemistry 101, 133–149 (2010).Article 
CAS 

Google Scholar 
Petraglia, A. et al. Litter decomposition: effects of temperature driven by soil moisture and vegetation type. Plant Soil 435, 187–200 (2019).Article 
CAS 

Google Scholar 
Vukicevich, E., Lowery, T., Bowen, P., Úrbez-Torres, J. R. & Hart, M. Cover crops to increase soil microbial diversity and mitigate decline in perennial agriculture. A review. Agron. Sustain. Dev. (2016).Bani, A. et al. The role of microbial community in the decomposition of leaf litter and deadwood. Appl. Soil. Ecol. 126, 75–84 (2018).Article 

Google Scholar 
Bonanomi, G., Capodilupo, M., Incerti, G., Mazzoleni, S. & Scala, F. Litter quality and temperature modulate microbial diversity effects on decomposition in model experiments. Community Ecol. 16, 167–177 (2015).Article 

Google Scholar 
Daebeler, A. et al. Pairing litter decomposition with microbial community structures using the Tea Bag Index (TBI). SOIL Discuss. [preprint]; 10.5194/soil-2021-110 (2021).Keuskamp, J. A., Dingemans, B. J. J., Lehtinen, T., Sarneel, J. M. & Hefting, M. M. Tea Bag Index: a novel approach to collect uniform decomposition data across ecosystems. Methods Ecol. Evol. 4, 1070–1075 (2013).Article 

Google Scholar 
Schaller, K. Praktikum zur Bodenkunde und Pflanzenernährung. Hochschule Geisenheim, (2000).Klindworth, A. et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 41, e1 (2013).Article 
CAS 
PubMed 

Google Scholar 
Ihrmark, K. et al. New primers to amplify the fungal ITS2 region–evaluation by 454-sequencing of artificial and natural communities. FEMS Microbiol. Ecol. 82, 666–677 (2012).Article 
CAS 
PubMed 

Google Scholar 
Schoch, C. L. et al. SI: Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi. Proc. Natl. Acad. Sci. U.S.A. 109, 6241–6246 (2012).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data. Available at https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (2010).Joshi, N. A. & Fass, J. N. sickle – A Windowed Adaptive Trimming Tool for FASTQ Files Using Quality. Available at https://github.com/najoshi/sickle (2011).Schloss, P. D. et al. Introducing mothur: Open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541 (2009).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Quast, C. et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596 (2013).Article 
CAS 
PubMed 

Google Scholar 
Westcott, S. L. & Schloss, P. D. OptiClust, an improved method for assigning amplicon-based sequence data to operational taxonomic units. mSphere 2, e00073 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Cole, J. R. et al. Ribosomal database project: data and tools for high throughput rRNA analysis. Nucleic Acids Res. 42, D633–D642 (2014).Article 
CAS 
PubMed 

Google Scholar 
Gweon, H. S. et al. PIPITS: an automated pipeline for analyses of fungal internal transcribed spacer sequences from the Illumina sequencing platform. Methods Ecol. Evol. 6, 973–980 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Kõljalg, U. et al. Towards a unified paradigm for sequence-based identification of fungi. Mol. Ecol. 22, 5271–5277 (2013).Article 
PubMed 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing. Available at https://www.R-project.org/ (2019).McMurdie, P. J. & Holmes, S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8, e61217 (2013).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Haegeman, B. et al. Robust estimation of microbial diversity in theory and in practice. ISME J. 7, 1092–1101 (2013).Article 
PubMed 
PubMed Central 

Google Scholar 
Scheu, S. Automated measurement of the respiratory response of soil microcompartments: Active microbial biomass in earthworm faeces. Soil Biol. Biochem. 24, 1113–1118 (1992).Article 

Google Scholar 
Mori, T. Validation of the Tea Bag Index as a standard approach for assessing organic matter decomposition: a laboratory incubation experiment. Ecol. Ind. 141, 109077 (2022).Article 
CAS 

Google Scholar 
Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D. & R Core Team. nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1–142. Available at https://CRAN.R-project.org/package=nlme (2019).Lenth, R. Emmeans: Estimated Marginal Means, Aka Least-Squares Means. R package version 1.4.4. Available at https://CRAN.R-project.org/package=emmeans (2020).Lefcheck, J. S. piecewiseSEM: Piecewise structural equation modelling in r for ecology, evolution, and systematics. Methods Ecol. Evol. 7, 573–579 (2016).Article 

Google Scholar 
Grace, J. B., Anderson, T. M., Olff, H. & Scheiner, S. M. On the specification of structural equation models for ecological systems. Ecol. Monogr. 80, 67–87 (2010).Article 

Google Scholar 
Shipley, B. A new inferential test for path models based on directed acyclic graphs. Struct. Equ. Model. 7, 206–218 (2000).Article 
MathSciNet 

Google Scholar  More

Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. Science 293, 2248–2251 (2001).Article 
CAS 

Google Scholar 
Seebacher, F., White, C. R. & Franklin, C. E. Physiological plasticity increases resilience of ectothermic animals to climate change. Nat. Clim. Change 5, 61–66 (2015).Article 

Google Scholar 
Havird, J. C. et al. Distinguishing between active plasticity due to thermal acclimation and passive plasticity due to Q10 effects: why methodology matters. Funct. Ecol. 34, 1015–1028 (2020).Article 

Google Scholar 
Dillon, M. E., Wang, G. & Huey, R. B. Global metabolic impacts of recent climate warming. Nature 467, 704–706 (2010).Article 
CAS 

Google Scholar 
White, C. R., Alton, L. A., Bywater, C. L., Lombardi, E. J. & Marshall, D. J. Metabolic scaling is the product of life history optimization. Science 377, 834–839 (2022).Article 
CAS 

Google Scholar 
Savage, V. M., Gilloly, J. F., Brown, J. H. & Charnov, E. L. Effects of body size and temperature on population growth. Am. Nat. 163, 429–441 (2004).Article 

Google Scholar 
Bernhardt, J. R., Sunday, J. M. & O’Connor, M. I. Metabolic theory and the temperature–size rule explain the temperature dependence of population carrying capacity. Am. Nat. 192, 687–697 (2018).Article 

Google Scholar 
Damuth, J. Population density and body size in mammals. Nature 290, 699–700 (1981).Article 

Google Scholar 
Schuster, L., Cameron, H., White, C. R. & Marshall, D. J. Metabolism drives demography in an experimental field test. Proc. Natl Acad. Sci. USA 118, e2104942118 (2021).Article 
CAS 

Google Scholar 
Amarasekare, P. & Coutinho, R. M. The intrinsic growth rate as a predictor of population viability under climate warming. J. Anim. Ecol. 82, 1240–1253 (2013).Article 

Google Scholar 
Amarasekare, P. & Savage, V. A framework for elucidating the temperature dependence of fitness. Am. Nat. 179, 178–191 (2012).Article 

Google Scholar 
Lande, R. Risks of population extinction from demographic and environmental stochasticity and random catastrophes. Am. Nat. 142, 911–927 (1993).Article 

Google Scholar 
Comeault, A. A. & Matute, D. R. Temperature-dependent competitive outcomes between the fruit flies Drosophila santomea and Drosophila yakuba. Am. Nat. 197, 312–323 (2021).Article 

Google Scholar 
Davis, A. J., Jenkinson, L. S., Lawton, J. H., Shorrocks, B. & Wood, S. Making mistakes when predicting shifts in species range in response to global warming. Nature 391, 783–786 (1998).Article 
CAS 

Google Scholar 
Davis, A. J., Lawton, J. H., Shorrocks, B. & Jenkinson, L. S. Individualistic species responses invalidate simple physiological models of community dynamics under global environmental change. J. Anim. Ecol. 67, 600–612 (1998).Article 

Google Scholar 
Gilman, S. E., Urban, M. C., Tewksbury, J., Gilchrist, G. W. & Holt, R. D. A framework for community interactions under climate change. Trends Ecol. Evol. 25, 325–331 (2010).Article 

Google Scholar 
Janča, M. & Gvoždík, L. Costly neighbours: heterospecific competitive interactions increase metabolic rates in dominant species. Sci. Rep. 7, 5177 (2017).Article 

Google Scholar 
Pettersen, A. K., Hall, M. D., White, C. R. & Marshall, D. J. Metabolic rate, context-dependent selection, and the competition–colonization trade-off. Evol. Lett. 4, 333–344 (2020).Article 

Google Scholar 
DeLong, J. P., Hanley, T. C. & Vasseur, D. A. Competition and the density dependence of metabolic rates. J. Anim. Ecol. 83, 51–58 (2014).Article 

Google Scholar 
Reid, D., Armstrong, J. D. & Metcalfe, N. B. Estimated standard metabolic rate interacts with territory quality and density to determine the growth rates of juvenile Atlantic salmon. Funct. Ecol. 25, 1360–1367 (2011).Article 

Google Scholar 
Ayala, F. J. in Essays in Evolution and Genetics in Honor of Theodosius Dobzhansky (eds Hecht, M. K. & Steere, W. C.) 121–158 (Springer, 1970).Atkinson, W. D. & Shorrocks, B. Aggregation of larval Diptera over discrete and ephemeral breeding sites: the implications for coexistence. Am. Nat. 124, 336–351 (1984).Article 

Google Scholar 
McKenzie, J. A. & McKechnie, S. W. A comparative study of resource utilization in natural populations of Drosophila melanogaster and D. simulans. Oecologia 40, 299–309 (1979).Article 
CAS 

Google Scholar 
Alton, L. A. et al. Developmental nutrition modulates metabolic responses to projected climate change. Funct. Ecol. 34, 2488–2502 (2020).Article 

Google Scholar 
Mitchell, K. A. & Hoffmann, A. A. Thermal ramping rate influences evolutionary potential and species differences for upper thermal limits in Drosophila. Funct. Ecol. 24, 694–700 (2010).Article 

Google Scholar 
Overgaard, J., Kristensen, T. N., Mitchell, K. A. & Hoffmann, A. A. Thermal tolerance in widespread and tropical Drosophila species: does phenotypic plasticity increase with latitude? Am. Nat. 178, S80–S96 (2011).Article 

Google Scholar 
Kellermann, V. et al. Comparing thermal performance curves across traits: how consistent are they? J. Exp. Biol. 222, jeb193433 (2019).Article 

Google Scholar 
Terblanche, J. S., Clusella-Trullas, S. & Chown, S. L. Phenotypic plasticity of gas exchange pattern and water loss in Scarabaeus spretus (Coleoptera: Scarabaeidae): deconstructing the basis for metabolic rate variation. J. Exp. Biol. 213, 2940–2949 (2010).Article 

Google Scholar 
Tewksbury, J. J., Huey, R. B. & Deutsch, C. A. Putting the heat on tropical animals. Science 320, 1296–1297 (2008).Article 
CAS 

Google Scholar 
Bos, M., Burnet, B., Farrow, R. & Woods, R. A. Mutual facilitation between larvae of the sibling species Drosophila melanogaster and D. simulans. Evolution 31, 824–828 (1977).Article 
CAS 

Google Scholar 
Arthur, W. On the complexity of a simple environment: competition, resource partitioning and facilitation in a two-species Drosophila system. Phil. Trans. R. Soc. B 313, 471–508 (1986).
Google Scholar 
Hodge, S., Mitchell, P. & Arthur, W. Factors affecting the occurrence of facilitative effects in interspecific interactions: an experiment using two species of Drosophila and Aspergillus niger. Oikos 87, 166–174 (1999).Article 

Google Scholar 
Bath, E., Morimoto, J. & Wigby, S. The developmental environment modulates mating-induced aggression and fighting success in adult female Drosophila. Funct. Ecol. 32, 2542–2552 (2018).Article 

Google Scholar 
Thibert, J., Farine, J. P., Cortot, J. & Ferveur, J. F. Drosophila food-associated pheromones: effect of experience, genotype and antibiotics on larval behavior. PLoS ONE 11, e0151451 (2016).Article 

Google Scholar 
Chown, S. L. et al. Scaling of insect metabolic rate is inconsistent with the nutrient supply network model. Funct. Ecol. 21, 282–290 (2007).Article 

Google Scholar 
Becker, R. A., Wilks, A. R. & Brownrigg, R. mapdata: extra map databases. R version 2.3.0 https://CRAN.R-project.org/package=mapdata (2018).R Core Team R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2022).Brooks, M. E. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9, 378–400 (2017).Article 

Google Scholar 
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
Bolker, B. & R Development Core Team bbmle: tools for general maximum likelihood estimation. R version 1.0.25 https://CRAN.R-project.org/package=bbmle (2022).Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest package: tests in linear mixed effects models. J. Stat. Softw. 82, 1–26 (2017).Article 

Google Scholar 
Fox, J. & Weisberg, S. An R Companion to Applied Regression 3rd edn (Sage, 2019).Hartig, F. DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. R version 0.4.6 https://CRAN.R-project.org/package=DHARMa (2022).Messamah, B., Kellermann, V., Malte, H., Loeschcke, V. & Overgaard, J. Metabolic cold adaptation contributes little to the interspecific variation in metabolic rates of 65 species of Drosophilidae. J. Insect Physiol. 98, 309–316 (2017).Article 
CAS 

Google Scholar 
Chamberlain, S. et al. rgbif: interface to the global biodiversity information facility API. R version 3.7.3 https://CRAN.R-project.org/package=rgbif (2022).Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
Hijmans, R. J. raster: geographic data analysis and modeling. R version 3.6-3 https://CRAN.R-project.org/package=raster (2022).Alton, L. A. & Kellermann, V. Data for “Interspecific interactions alter the metabolic costs of climate warming”. Zenodo https://doi.org/10.5281/zenodo.7475922 (2023).White, C. R. et al. Geographical bias in physiological data limits predictions of global change impacts. Funct. Ecol. 35, 1572–1578 (2021).Article 

Google Scholar  More

Terraces are a land type that is defined by its shape. They have a distinct morphological structure and edge features that distinguish them from other land types. In this study, we define terraces as agricultural land with strip or wavy sections built on slopes greater than 2° along the contour direction. Figure 1 depicts Google Maps satellite images of terraces in the Loess Plateau region. Terraces can be distinguished from other features in remote sensing images based on their colour, morphology, texture, and structure. Terraces can be distinguished from construction land, water, glaciers, and deserts by their colours. Figure 1b–d shows terraces that are primarily green and yellow. Furthermore, terraces are generally distributed along the contour direction, and can therefore be identified based on their morphology. Terraced field ridges curve downward and resemble strips in Fig. 1b,d or circles or ovals in Fig. 1c rather than a neat grid-like distribution. These features differ in morphology from the flat land shown in Fig. 1h. Based on texture and structure, the field area of terraces can be identified based on their strong edge features, as shown in Fig. 1b–d. The edges of terraces have dark stripes caused by oblique illumination received from the sun, and the field ridge of terraces often intercepts part of the sunlight due to their height. Sloping cultivated land, as shown in Fig. 1g, has no evident terraced wall. The outline of sloping cultivated land in the high-resolution image is curved, with no prominent edge features. These findings are critical differences distinguishing terraces and sloping land in high-resolution images.Fig. 1The spatial location of the Loess Plateau and images of various types of cultivated land. (a) The spatial location of the Loess Plateau and Spatial distribution of various cultivated land types images, (b) wide strip-mounted terraces in Longxi, (c) circular wide terraces in central Yulin, (d) high resolution image of Zhuanglang County in July 2019, (e) Zhuanglang County in February 2020, (f) narrow terraces in Shangbao, Chongyi, Jiangxi Province, (g) sloping cropland in Zhenjiang Town, Laibin, Guangxi, and (h) horizontal cropland in the North China Plain.Full size imageDeep learning-based terrace extraction modelThe DLTEM is a terrace extraction model that uses deep learning algorithms and other supplementary information. Initially, a preliminary terrace distribution map was obtained using a deep learning algorithm. It was then combined with the spectral and digital elevation model (DEM) elevation information to fine-tune the results. The final spatial distribution of the terraces was produced by manual correction (Fig. 2). Traditional land classification models or methods typically superimpose spectral, elevation, and morphological texture information from remote sensing images together for training, such as random forest, which is easily ignored in training since morphological texture information accounts for a relatively small amount of the total information. This leads to significant errors while identifying land classes with textural characteristics. In contrast, the DLTEM focuses on morphological texture information from remote sensing images and classifies it into land classes, followed by auxiliary correction through additional information. Thus, this method is more suitable to extract terraces enriched with texture structure information.Fig. 2Flow chart of the deep learning-based terrace extraction model.Full size imageThe UNet++ network is a classic deep learning algorithm that is uniquely unrivaled in extracting colour, morphology, texture, and structure features from images and applying them for classification. In comparison with other Convolutional Neural Network (CNN) classification models (e.g., Fully Convolutional Networks (FCN)), it has high classification accuracy, fast computation speed, strong robustness, and provides variable importance metrics. Therefore, in this study, the UNet++ network was adopted as the network framework for deep learning; the primary data source used was high-resolution satellite imagery from 2019. DEM (SRTM v4.1) data were used to obtain the elevation information and GlobeLand30 data were used to obtain the spectral information. The results were corrected to construct the final map of the distribution of terraces in the Loess Plateau.Study areaThe Loess Plateau, one of China’s four major plateaus, is located in northern central China (34°–40° N and 103°–114° E) (Fig. 1). It is covered by a thick loess layer that ranges in thickness from 50 to 80 m, and is the world’s largest loess deposition area, covering 648,700 km2. The altitude of the Loess Plateau ranges from 800 to 3,000 m, its average annual temperature is 6–14 °C, and its average annual precipitation is 200–700 mm. Since ancient times, the Loess Plateau has been used for agriculture because of its fine grains, fluffy soil texture, and rich soluble mineral nutrients, all of which are conducive to crop cultivation. However, long-term unsustainable land use caused the degradation of the vegetation cover in the Loess Plateau. Moreover, the land is degrading due to considerable nutrient loss caused by long-term water erosion in conjunction with natural conditions, such as arid climate, loose soil, concentrated and heavy rainfall. The fragmented ground in the region has made it susceptible to soil erosion. It has also become the primary source of Yellow River sediment as a result of the massive flow of eroded sediment into the Yellow River, posing a serious threat to the economic and social development of the lower Yellow River basin.Terracing is one of the main measures used to enhance crop yield and conserve soil and water in the region. Since the 1980s, the Chinese government has implemented many large-scale slope-to-terrace projects in the Loess Plateau. Especially in recent years, the outline of the comprehensive management plan for the Loess Plateau area (2010–2030) has been promulgated with a planned area of 2.608 million hectares for slope to terrace conversion, making it the core area of slope to terrace conversion projects in the country.Data preparationAlthough high-resolution satellite images can be an important data source for the spatial distribution of terraces on the Loess Plateau, they are not ideal for terraces classification. On the one hand, a higher resolution image requires more storage space. On the other hand, it reduces the efficiency, prolongs the interpretation time, and increases the noise in the image, affecting the interpretation accuracy. Most of the terraces on the Loess Plateau are wider than 7 m (Fig. 1b–d). These are wide terraces in comparison with the narrow terraces of southern China (Fig. 1f), which are less than 2 m wide. Furthermore, it is also easy to mistake the fish-scale pits constructed for soil and water conservation for terraces because of their similarity in form. However, as the width of their field surface is less than 1.5 m, remote sensing images with a 2 m resolution can effectively prevent the false extraction of such features. Based on the actual situation of this study area, we chose a high-resolution image with a spatial resolution of 1.89 m from Google Maps 16 level as the data source. The colour, texture, and morphological features of terraces in the images show seasonal variations. In autumn and winter, the weather is dry, and the vegetation is less shaded in the Loess Plateau. During this time, even the edge features become more visible and easier to identify. As a result, we selected images from October 2018 to February 2019 whenever possible (Fig. 1c,d).Deep learning network selectionLand classification is the extraction of land types from remote sensing images using image segmentation techniques. As the key technology of image segmentation, the Fully Convolutional Network (FCN) classifies images at the pixel level. FCN follows the network structure pattern of encoding and decoding, which adopts AlexNet as the encoder of the network and then employs transposed convolution to up-sample the feature map output from the final convolutional layer of the encoder to the resolution of the input image to achieve pixel-level image segmentation. However, due to the large error in image pixel boundary localization, Ronneberger et al.29 improved the FCN structure in 2015 by expanding the capacity of the network decoder by adding a contracting path to the encoding and decoding modules to achieve more accurate pixel boundary localisation29. The U-Net network is commonly used in medical image processing because it requires a small number of training samples and is effective in classifying objects with a fixed structure and limited semantic information. This network is comparable to natural image semantic segmentation such as Deeplab v3+, which has a smaller number of model parameters and the same effect.Since the texture and morphological features of terraces and human organs have certain similarities, they are primarily manifested by simple semantic information contained within the terrace images themselves. Thus, high-level semantic information and low-level features of such images become more important. However, high-resolution images are more complicated and variable than medical image patterns, and errors in terrace extraction edge identification using the U-Net network, such as boundary segmentation of terraces and flatlands, still occur. To fully utilize the semantic information of the network, we adopted a nested U-Net architecture, namely the UNet++ network proposed by Zhou et al.28. The network integrates long-connected and short-connected architectures to capture features at different levels by adding a shallower U-Net structure and integrates them via feature superposition to make the scale difference of feature maps smaller when fused to enhance the correct rate of image segmentation edges. However, because the U-Net++ network increases the number of model parameters, this study adopted the sparse matrix approach to accelerate model training and decrease the number of parameters.Data pre-processingData pre-processing is a prerequisite for UNet++ network training, that is, valid input according to the standard format annotation before training can be performed. Since the UNet++ network proposed by Zhou et al.28. is primarily used for medical images, which have characteristics such as fixed image structure, no spatial information, and less pattern variation, labelling medical images is comparatively easier using this method. In contrast, high-resolution remote sensing images have a large number of rasters, many pattern changes, irregular image structure, and spatial information. Therefore, determining how to better annotate high-resolution remote sensing images and reduce the annotation workload becomes critical. First, we vectorized the training sample area and generated the terrace vector dataset using ArcGIS with a high-resolution remote sensing image as the primitive map. Second, we converted the terrace vector dataset into raster data. The information of the raster had to be identical to that of the primitive map, including the size of the raster, its processing range, and its coordinate system. The output was converted to TIFF format to complete the image annotation. Since the raster size input to UNet++ network training is a fixed size, it is much smaller than the original image. To simplify the process of inputting the original image and its annotation information, we added an image import module to DLTEM, which was a sliding window of 400*400, and read the image automatically by setting the corresponding judgement conditions. Finally, the entire high-resolution image was processed automatically into the model in accordance with the established rules for training.The goal of the data enhancement was to improve the universality and robustness of the UNet++ network training results. As mentioned above, the high-resolution images taken simultaneously often included clouds or other anomalies in some areas, as the images were stitched together using multiple sources of data fusion. This can easily form evident stitching traces (Fig. 1c,d) due to the different shooting times and image quality of various data sources, i.e., brightness, saturation, and colour contrast of the images. Thus, the model trained on the original image data has strong limitations, and in many scenes, there are notable matrix-type misclassification regions due to image differences, making extraction work challenging. Therefore, in this study, we first adjusted the brightness, grayscale, and contrast of the training data after input to enhance its colour feature recognition ability. We then altered the scaling of the image, and rotated and transformed the training image from 0° to 360° to enhance morphological feature recognition and the accuracy of the training network in terrace extraction.Parameter settingThe network parameter setting is the most critical hyperparameter for UNet++ network training. They are mainly divided into input image size, batch size, learning rate, number of iterations, objective function, gradient descent strategy, momentum, decay rate, and activation function. Among them, we set the image size to 400*400 pixels based on the actual situation of the terraced area, where the UNet++ network has four scaling times, and the image size must be a multiple of 16. The batch size primarily affects the convergence of the model. If the batch limit is set to one, the model is easily affected by the random perturbation phenomenon and cannot converge to find the optimal solution. Since the batch size is determined by the size of the video memory, the value of the batch is limited by equipment constraints. The model in this study used a 2080Ti video card with 11 GB of video memory, and the batch was set to 8. The learning rate, gradient descent strategy, and objective function play a role in whether the network can find the best classification model better and faster. The learning rate was set to 0.001 for the first 500 generations, with the goal of achieving fast convergence to the target region. The learning rate was then set to 0.0001 for 500–1,000 generations, and the model was fine-tuned by choosing a smaller learning rate to find the model with the highest classification accuracy. Adam was chosen for the gradient descent strategy. The momentum and adaptive learning rate were used to increase the convergence rate. The cross-entropy classification loss function was chosen as the objective function to improve the differentiation between terraced and non-terraced areas. Momentum, decay rate, and activation function were all adopted from the previous default settings of the UNet++ network.Data correctionIn this study, we primarily used high-resolution images from Google Earth as the data source to extract the distribution of terraces on the Loess Plateau. Because this image source only contains a large amount of texture structure information and no vegetation information, it is easy to misjudge and misclassify features with the same morphological structure and edge features, such as permanent snow and ice, water bodies, bare land, and artificial surfaces. Vegetation information was generally processed based on waveband data from multispectral/hyperspectral images. It requires topographic correction, atmospheric correction, radiometric calibration, de-clouding, and other operational processes, which are extremely sophisticated30.GlobeLand30 is a 30 m spatial resolution global surface coverage dataset developed by the National Geomatics Center of China. The most recent GlobeLand30 dataset (v2020) has been updated with data sources from 2017 to the present. Its extensive data sources enable effective reduction of the impacts of cloud cover, with an overall accuracy of 85.72%. The classification accuracy of permanent snow and ice, water bodies, bare land, and artificial surfaces of this dataset is as high as 75.79%, 84.70%, 81.76%, and 86.70%, respectively. Since the update time of v2020 data is similar to that of high-resolution images, it can be used as correction data for vegetation information31.Since the training image data are two-dimensional planar data with no elevation or slope information (Fig. 1g), certain flat fields with visible field bumps are easily misclassified as terraces. The Space Shuttle Radar Topography Mission (SRTM v4.1) DEM has a spatial resolution of 30 m and ranges from 60° N to 56° S, completely covering the Loess Plateau32,33. In this study, these data were treated as terrain correction data. The amendment standard corrects the areas that have been extracted as terraces below 2° to non-terraced areas according to the requirements of the Ministry of Natural Resources of China.The spatial resolution of our extracted terraces is 1.89 m, whereas the spatial resolution of GlobeLand30 and DEM as correction data sources is 30 m, which is difficult to meet the requirements of data processing. Hence, we up-sample the two correction data sources, and then used multi-source data fusion. First, we extracted and up-sampled the terraced areas of glaciers, rivers, and deserts from GlobeLand30 to a spatial resolution of 1.89 m. Secondly, we up-sampled the DEM to 1.89 m using spatial interpolation for its raster centre as the true value of the region and performed a slope calculation for the up-sampled DEM. Further, the spatial distribution maps of glaciers, rivers, deserts, and slope maps of the Loess Plateau with the same resolution as the spatial distribution maps of terraces were available. Finally, we superimposed these images, used the terrace range in the TDMLP as a mask, and assessed the pixels in the mask area one by one. If a pixel belonged to permanent snow and ice, a water body, bare land, or an artificial surface, or had a slope less than 2°, it was modified to the background value. Otherwise, the original value was retained.We made artificial corrections to the data based on the extracted results for the arid areas of the Loess Plateau as well as for the flatter basins, given that these areas do not feature terraces.Training and validation dataFor supervised classification, the selection of sample areas and sample features is crucial. The focus and core of any land classification work is representative and effective training sample selection. To obtain a better sample area selection, we considered the selection of sample areas from three perspectives, i.e., colour texture features, topographic features, and spatial distance of the training samples. First, the terraces in this study are in agricultural land, including cultivated land, woodland, grassland, and other types of land; thus, different types of land will present different texture details. At the same time, high-resolution images from Google Earth are mosaicked. Because of the different acquisition times, the same region and land type will have visible colour differences and stitching traces, which is more common in the Loess Plateau region. Therefore, these factors should be considered in the selection of training samples as much as possible to improve the generability of the model and the correct rate of its extraction. Second, the state of the terraces varies according to topographic features. Among them, gradient, direction, altitude, and climate are the most significant factors. Terraces can be categorised as shallow-slope or steep-slope terraces. Based on slope aspect, altitude, and climate characteristics, they can also be categorised as either easy to identify or hard to identify. Thus, the sample should be inclusive of these types of terraces. According to the first law of geography, terraces in different spatial locations have different morphologies. Therefore, the spatial location of the samples should also be at a certain distance.In summary, we selected one county in each region based on the geomorphic zoning characteristics of the Loess Plateau. In addition, we added one more in the area where the density of terraces may be higher. Finally, we selected the whole area of seven counties (Fig. 3) as the training sample area distribution, covering 2.18% of the overall Loess Plateau area. The colour morphological features, topographic features, spatial location, and imaging quality of terrace images in these regions are highly representative. This method was unique from other classification methods. Most of the traditional methods are based on the single-pixel information of feature layers such as random forests, which tend to ignore the neighbouring information around the point, and thus are subject to misclassification and under classification for land types with outstanding texture information. In our study, we adopted the visual interpretation of the whole domain, which can cover the neighbourhood information of each pixel point more comprehensively. To ensure the uniformity and correctness of visual interpretation, the terraces in the training area were visually interpreted by seven interpreters after uniform professional training. For the disputed and uncertain areas, the seven interpreters carried out interactive interpretation and scoring according to the interpretation results. Finally, two other interpretation experts made the final review and corrections. The interpretation results of the training area were re-examined and revised based on the results of the later interpretations.Fig. 3Distribution of training sample areas and validation sites in terraces on the Loess Plateau.Full size imageTo better assess and compare the validity and correctness of the terraced agricultural area datasets on the Loess Plateau in quantitatively, the validation dataset was divided into two parts: a per-pixel point-based validation set and a field validation dataset of terraces with location information. The extracted datasets were comprehensively evaluated in terms of both pixel scale and field validation.We constructed a single-pixel validation point that evaluates the TDMLP. We applied the Icosahedral Snyder Equal Area Discrete Global Grid created by ArcGIS. Based on this strategy, the study area was partitioned into 972 regions (Fig. 3). To better validate the terrace classification results (excluding non-terrace classes), we placed more validation points within the grid where the terrace distribution is more concentrated. First, we calculated the proportion of terraces in each hexagonal grid to the total area of the hexagonal grid. Second, we separated the terraces into four levels according to the proportion of terraces to the whole grid area as 0–20%, 20–50%, 50–80%, and 80–100% and the number of validation points was 10, 20, 40, and 50, respectively.Since the proportion of the extracted terraced area to the total area was only 14%, direct random point deployment would have led to fewer terraced validation sets and thus would have affected the final data evaluation. Therefore, in the deployment strategy, we ensured that the validation points distributed in the extracted terraces in each grid account for at least one-fifth of the total number of validation points, but for the grid with a smaller proportion of terraces or even 0, this practice was meaningless. Hence, we stipulated that in the grid with a proportion of terraces ≤1%, direct random scattering was to be performed. The final scattered verification points in the terraced and non-terraced areas were 5,194 and 6,226, respectively, with a ratio close to 1:1 for easy verification. The spatial distribution is shown in Fig. 3.We validated the spatial distribution map of terraces on the Loess Plateau from 14 April 2021 to 1 May 2021 and constructed a field validation dataset of terraces with location information. Considering the longitudinal, latitudinal, and vertical heterogeneities of the Loess Plateau, the verification route was divided into two sections, north to south and east to west, to more comprehensively cover all regions of the Loess Plateau. The verification route started at Hohhot in the northeast of the Loess Plateau. It passed through the Datong Basin, followed the Yellow River to the south and the Weihe Plain, and then travelled westward through Mount Liupan to the westernmost part of the Loess Plateau. The route was through 54 counties/districts in 16 cities and six provinces on the Loess Plateau, with a total distance of 3,680 km, covering 15.8% of the counties on the Loess Plateau (total of 341 counties). We also surveyed and sampled the verification points approximately every 5 km along the route and collected data from a total of 815 sample points, covering various types of terraces on the Loess Plateau. The results are shown in Fig. 3. More

All procedures involving animals were conducted in accordance with the guidelines and regulations from Institutional Animal Care and Use Committee (IACUC) of the University of Florida (protocol #201509019). Tis manuscript is reported in accordance with ARRIVE guidelines.Site descriptionThis study was carried out at the North Florida Research and Education Center, in Marianna, FL (30°46′35″N 85°14′17″W, 51 m.a.s.l). The trial was performed in two experimental years (2019 and 2020) in a greenhouse.The soil used was collected from a pasture of rhizoma peanut (Arachis glabrata Benth.) and Argentine bahiagrass (Paspalum notatum Flügge) as the main forages. Without plant and root material, only soil was placed into buckets, as described below in the bucket assemblage section. Soil was classified as Orangeburg loamy sand (fine-loamy-kaolinitic, thermic Typic Kandiudults), with a pHwater of 6.7, Mehlich-1-extratable P, K, Mg and Ca concentrations of 41, 59, 63, 368 mg kg−1, respectively. Average of minimum and maximum daily temperature and relative humidity in the greenhouse for September and November (September for beetle trial due seasonal appearance of beetles, and October and November to the Pear Millet trial) in 2019 and 2020 were 11 and 33 °C, 81%; 10 and 35 °C, 77%, respectively.Biological material determinationTo select the species of beetles, a previous dung beetle sampling was performed in the grazing experiment in the same area (grass and legume forage mixture) to determine the number of dung beetle species according to the functional groups as described by Conover et al.44. Beetles were pre-sampled from March 2017 to June 2018, where Tunnelers group were dominant and represented by Onthophagus taurus (Schreber), Digitonthophagus gazella (Fabricius), Phanaeus vindex (MacLeay), Onthophagus oklahomensis (Brown), and Euniticellus intermedius (Reiche). Other species were present but not abundant, including Aphodius psudolividus (Linnaeus), Aphodius carolinus (Linnaeus), and Canthon pilularius (Linnaeus) identified as Dweller and Roller groups, respectively. The pre-sampling indicated three species from the Tunneler group were more abundant, and thereby, were chosen to compose the experimental treatments (Fig. 4).Figure 4Most abundant dung beetle species in Marianna, FL used in the current study. Credits: Carlos C.V. García.Full size imageBeetles collection and experimental treatmentsThree species of common communal dung beetles were used: O. taurus (1), D. gazella (2), and P. vindex (3). Treatments included two treatments containing only soil and soil + dung without beetles were considered as Control 1 (T1) and Control 2 (T2), respectively. Isolated species T3 = 1, T4 = 2, T5 = 3 and their combinations T6 = 1 + 2 and T7 = 1 + 2 + 3. Dung beetles were trapped in the pasture with grazing animals using the standard cattle-dung-baited pitfall traps, as described by Bertone et al.41. To avoid losing samples due to cattle trampling, 18 traps were randomized in nine paddocks (two traps per paddock) and installed protected by metal cages, and after a 24-h period, beetles were collected, and the traps removed. Table 1 shows the number of dung beetles, their total mass (used to standardize treatments) per treatment, and the average mass per species. To keep uniformity across treatments we kept beetle biomass constant across species at roughly 1.7 to 1.8 g per assemblage (Table 1). Twenty-four hours after retrieving the beetles from the field traps, they were separated using an insect rearing cage, classified, and thereafter stored in small glass bottles provided with a stopper and linked to a mesh to keep the ventilation and maintaining the beetles alive.Table 1 Total number and biomass of dung beetles per treatment.Full size tableBuckets assemblageThe soil used in the buckets was collected from the grazing trial in two experimental years (August 2019 and August 2020) across nine paddocks (0.9 ha each). The 21 plastic buckets had a 23-cm diameter and 30-cm (0.034 m2) and each received 10 kg of soil (Fig. 5). At the bottom of the recipient, seven holes were made for water drainage using a metallic mesh with 1-mm diameter above the surface of the holes to prevent dung beetles from escaping. Water was added every four days to maintain the natural soil conditions at 60% of the soil (i.e., bucket) field capacity (measured with the soil weight and water holding capacity of the soil). Because soil from the three paddocks had a slightly different texture (sandy clay and sandy clay loam), we used them as the blocking factor.Figure 5Bucket plastic bucket details for dung beetle trial.Full size imageThe fresh dung amount used in the trial was determined based on the average area covered by dung and dung weight (0.05 to 0.09 m2 and 1.5 to 2.7 kg) from cattle in grazing systems, as suggested by Carpinelli et al.45. Fresh dung was collected from Angus steers grazing warm-season grass (bahiagrass) pastures and stored in fridge for 24 h, prior to start the experiment. A total of 16.2 kg of fresh dung was collected, in which 0.9 kg were used in each bucket. After the dung application, dung beetles were added to the bucket. To prevent dung beetles from escaping, a mobile plastic mesh with 0.5 mm diameter was placed covering the buckets before and after each evaluation. The experiment lasted for 24 days in each experimental year (2019 and 2020), with average temperature 28 °C and relative humidity of 79%, acquired information from the Florida Automated Weather Network (FAWN).Chamber measurementsThe gas fluxes from treatments were evaluated using the static chamber technique46. The chambers were circular, with a radius of 10.5 cm (0.034 m2). Chamber bases and lids were made of polyvinyl chloride (PVC), and the lid were lined with an acrylic sheet to avoid any reactions of gases of interest with chamber material (Fig. 6). The chamber lids were covered with reflective tape to provide insulation, and equipped with a rubber septum for sampling47. The lid was fitted with a 6-mm diameter, 10-cm length copper venting tube to ensure adequate air pressure inside the chamber during measurements, considering an average wind speed of 1.7 m s−148,49. During measurements, chamber lids and bases were kept sealed by fitting bicycle tire inner tubes tightly over the area separating the lid and the base. Bases of chambers were installed on top of the buckets to an 8-cm depth, with 5 cm extending above ground level. Bases were removed in the last evaluation day (24th) of each experimental year.Figure 6Static chamber details and instruments for GHG collection in the dung beetle trial.Full size imageGas fluxes measurementsThe gas fluxes were measured at 1000 h following sampling recommendations by Parkin & Venterea50, on seven occasions from August 28th to September 22nd in both years (2019 and 2020), being days 0, 1, 2, 3, 6, 12, and 24 after dung application. For each chamber, gas samples were taken using a 60-mL syringe at 15-min intervals (t0, t15, and t30). The gas was immediately flushed into pre-evacuated 30-mL glass vials equipped with a butyl rubber stopper sealed with an aluminium septum (this procedure was made twice per vial and per collection time). Time zero (t0) represented the gas collected out of the buckets (before closing the chamber). Immediately thereafter, the bucket lid was tightly closed by fitting the lid to the base with the bicycle inner tube, followed by the next sample deployment times.Gas sample analyses were conducted using a gas chromatograph (Trace 1310 Gas Chromatograph, Thermo Scientific, Waltham, MA). For N2O, an electron capture detector (350 °C) and a capillary column (J&W GC packed column in stainless steel tubing, length 6.56 ft (2 M), 1/8 in. OD, 2 mm ID, Hayesep D packing, mesh size 80/100, pre-conditioned, Agilent Technologies) were used. Temperature of the injector and columns were 80 and 200 °C, respectively. Daily flux of N2O-N (g ha−1 day−1) was calculated as described in Eq. (1):$${text{F}}, = ,{text{A}}*{text{dC}}/{text{dt}}$$
(1)
where F is flux of N2O (g ha−1 day−1), A is the area of the chamber, and dC/dt is the change of concentration in time calculated using a linear method of integration by Venterea et al.49.Ammonia volatilization measurementAmmonia volatilization was measured using the open chamber technique, as described by Araújo et al.51. The ammonia chamber was made of a 2-L volume polyethylene terephthalate (PET) bottle. The bottom of the bottle was removed and used as a cap above the top opening to keep the environment controlled, free of insects and other sources of contamination. An iron wire was used to support the plastic jar. A strip of polyfoam (250 mm in length, 25 mm wide, and 3 mm thick) was soaked in 20 ml of acid solution (H2SO4 1 mol dm−3 + glycerine 2% v/v) and fastened to the top, with the bottom end of the foam remaining inside the plastic jar. Inside each chamber there was a 250-mm long wire designed with a hook to support it from the top of the bottle, and wire basket at the bottom end to support a plastic jar (25 mL) that contained the acid solution to keep the foam strip moist during sampling periods (Fig. 7). The ammonia chambers were placed installed in the bucket located in the middle of each experimental block after the last gas sampling of the day and removed before the start of the next gas sampling.Figure 7Mobile ammonia chamber details for ammonia measurement in dung beetle trial. Adapted from Araújo et al.51.Full size imageNutrient cyclingPhotographs of the soil and dung portion of each bucket were taken twenty-four hours after the last day of gas flux measurement sampling to determine the dung removal from single beetle species and their combination. In the section on statistical analysis, the programming and statistical procedures are described. After this procedure, seeds of pearl millet were planted in each bucket. After 5 days of seed germination plants were thinned, maintaining four plants per bucket. Additionally, plants were clipped twice in a five-week interval, with the first cut occurring on October 23rd and the second cut occurring on November 24th, in both experimental years. Before each harvest, plant height was measured twice in the last week. In the harvest day all plants were clipped 10 cm above the ground level. Samples were dried at 55 °C in a forced-air oven until constant weight and ball-milled using a Mixer Mill MM 400 (Retsch, Newton, PA, USA) for 9 min at 25 Hz, and analyzed for total N concentration using a C, H, N, and S analyzer by the Dumas dry combustion method (Vario Micro Cube; Elementar, Hanau, Germany).Statistical analysisTreatments were distributed in a randomized complete block design (RCBD), with three replications. Data were analyzed using the Mixed Procedure from SAS (ver. 9.4., SAS Inst., Cary, NC) and LSMEANS compared using PDIFF adjusted by the t-test (P  More

Test organisms and exposuresIn this study, we used test organisms and reagents according to the Acute Toxicity Test Method of Daphnia magna Straus(Cladocera, Crustacea); ES 04704.1b29. Daphnia magna were fostered at the National Institute of Environmental Research and were adopted. During the test, adult female Daphnia magna over two weeks of age, cultured over several generations, were transferred to a freshly prepared container the day before the test. Daphnia magna are neonates for less than 24 h after birth29. To maintain the sensitivity of the organism, young individuals less than 24 h old that reproduced the following day were used. Individuals of a similar size were selected for the test. Daphnia magna was fed YCT, which is a mixture of green algae in Chlorella sp., yeast, Cerophy II(R), and trout chow. Sufficient amounts of prey were supplied 2 h before the test to minimize the effects of prey during the test. The test medium was prepared by dissolving KCl (8 mg/L), (text {MgSO}_4) (120 mg/L), (text {CaSO}_4 cdot 2 text {H}_2 text {O} ) (120 mg/L), and (text {NaHCO}_3) (192 mg/L) in deionized water.Automatic high-throughput Daphnia magna tracking systemTo build an automatic high-throughput Daphnia magna tracking system, we equipped the system with a video analysis algorithm as well as flow cells (Fig. 1). In the tracking system, six flow cells filled with culture medium were installed in the device. Each flow cell contained 10 Daphnia magna. Subsequently, to automatically measure the state of Daphnia magna, the six flow cells were photographed at 15 frames per second using a camera (Industrial Development Systems imaging) equipped with a CMOSIS sensor capable of infrared imaging. A red light close to the infrared spectrum was placed at the back of the flow cells for uniform illumination and to minimize stress on Daphnia magna. To capture the size and movement of the Daphnia magna as accurately as possible, the camera was set to a frame rate of 15 fps and a resolution of 2048 (times ) 1088 (2.23 MB), using a 12 mm lens. The distance between the flow cell and the camera was set to 16 cm. To measure the number of mobile Daphnia magna, their lethality, and swimming inhibition automatically and simultaneously, one camera for every two cells was used to collect the status data of Daphnia magna. For assessing ecotoxicity, the video analysis system used images obtained from the six flow cells to track each Daphnia magna and estimate key statistics such as the number of mobile individuals, average distance, and radius of activity.Figure 1New automatic high-throughput video tracking system for behavioral analysis using Daphnia magna as a model organismFull size imageThe automatic high-throughput video tracking system in the ecotoxicity measuring device was designed to continuously measure the ecotoxicity of Daphnia magna (Fig. 2). Daphnia magna moves faster at high temperatures and is less active at low temperatures. Thus, a constant temperature module that can be set to an appropriate Daphnia magna habitat temperature (20 ± 2 (^{circ })C) was added to create a suitable culture environment for Daphnia magna29. Natural pseudo-light ((lambda >590) nm, 3000 k) was installed on the upper part of the detector for proper habitat light intensity (500 Lux–1000 Lux). The size of the flow cell was set as small as possible while observing the movement of the Daphnia magna. An automatic feeding system was installed so that food could be injected during the replacement cycle. The six independent multi-flow cells were designed with an automatic dilution injection module; therefore, these flow cells were diluted to six different concentrations (100%, 50%, 25%, 12.5%, 6.25%, and 0%).Figure 2Schematic representation of the automatic high-throughput video tracking systemFull size imageAutomatic tracking algorithmThe CPU used for Daphnia magna tracking was Intel i5-9300H @ 2.40 GHz, with 8 GB of memory and Windows 10 Pro 64-bit operating system. In this experiment, the algorithms were trained using 12 Daphnia magna videos and tested using an additional four Daphnia magna videos. Subsequently, the detection and tracking methods were compared. The videos, each of which had a duration of 30 s, were captured at a rate of 15 frames per second. Generally, for long-time or real-time videos, the following factors must be considered in tracking Daphnia magna: automatic binarization between the object and background, effective classification of Daphnia magna or noise, and the speed of the algorithm. Therefore, to develop an efficient tracking algorithm, we propose the following tracking process (Fig. 3A). In this process, each frame is initially converted into an image and the background is identified from the obtained video (Fig. 3B). The background is the average of the frames over the previous 20 s, and the tracking system takes 20 s to capture the first background image. The background is subtracted from the image for object detection (Fig. 3C). The objects include Daphnia magna and noise such as droplets and sediment. The difference between the background and frame images is binarized, and each area of the binarized values is regarded as an object. Conventionally, the binarized values are manually generated using specific thresholds. In this study, the images are automatically binarized using k-means clustering to select the threshold value. After binarization, several machine learning methods are used to classify the objects as Daphnia magna or noise (Fig. 3D). For a faster tracking algorithm, we use simple machine learning methods such as random forest (RF) and support vector machine (SVM). The predicted Daphnia magna are tracked using SORT24, which is a fast and highly accurate tracking algorithm (Fig. 3E). Finally, based on the tracked results, statistics for assessing ecotoxicity, such as the number of mobile individuals, average distance, and radius of activity, are estimated to evaluate the toxicity of the aquatic environment.Figure 3Automatic Daphnia magna tracking algorithm process. (A) Overview of automatic tracking algorithm process. (B) Image extraction step. (C) Background subtraction step. (D) Daphnia magna detection step. (E) Daphnia magna tracking step.Full size imagek-means clustering for automatic background subtractionMany tracking algorithms assume that the background is fixed. With fixed backgrounds, the difference between the frame and background can be used to identify objects. However, automatically selecting the precise threshold value for image pixel binarization becomes one of the key problems in identifying objects. The proposed method applies k-means clustering to the pixel values of the subtracted image30, and the center value of each calculated cluster mean is selected as the threshold value (Fig. 4). In the k-means clustering method, grouping is repeatedly performed using the distance between data points31. For binarization, two groups are formed. Let (mu _1 (t)) be the mean of pixels less than the threshold and (mu _2(t)) be the mean of pixels greater than the threshold. At first, (mu _1(t), mu _2(t)) are randomly initialized. Subsequently, each pixel is grouped into a closer mean of each group. The above steps are repeated several times until the group experiences a few changes. Finally, the threshold is calculated as an average of the two means.Figure 4Example of automatic threshold value setting for binarization between objects and background using k-means clusteringFull size imageClassification methodsObject detection based solely on the subtraction between the background and frame images may have low accuracy. As the background in the proposed process is the average value of the frame images, noise may occur. Although this noise is removed by threshold selection in binarization, using only the threshold selection is not efficient for long or real-time videos. Therefore, additional noise must be classified and removed using machine learning models, requiring the construction of a database. In the database, the obtained objects are manually labeled as noise or Daphnia magna and are called ground truth. For classification, the resized 8 (times ) 8 image of each object is stored in the database. The resized image is transformed into a feature using the Sobel edge detection algorithm32 and entered as inputs to the classification models. In this study, classification models such as RF33 SVM34 were used.RF is a model that integrates several decision tree models35. All training data are sampled with a replacement for training each decision tree model. The decision tree model is trained to split intervals of each independent variable by minimizing the gini index (Eq. 1) or entropy index (Eq. 2). The gini index and entropy index denote the impurity within the intervals.$$begin{aligned} G= & {} 1- sum _{i=1}^{c} p_i ^2 end{aligned}$$
(1)
$$begin{aligned} E= & {} – sum _{i=1}^{c} p_i log_2 p_i end{aligned}$$
(2)
where (p_i) is a probability within i-th interval, and c is the number of intervals. For better performance, the RF selects independent variables of training data randomly. This step serves to reduce the correlation of each model. If predictions of each decision tree are uncorrelated, then the variance of an integrated prediction of models is smaller than the variance of each model. RF integrates several model predictions using the voting method. An advantage of the RF method is that it avoids overfitting because the model uses the average of many predictions.SVM is a model designed to search for a hyperplane to maximize the distance, or margin, between support vectors. The hyperplane refers to the plane that divides two different groups, and the support vector represents the closest vector to the hyperplane. Let (D=({textbf{x}}_i, y_i), i=1, ldots , n, {textbf{x}}_i in {mathbb {R}}^p, y_n in { -1,1 }) be training data. Suppose that the training data are completely separated linearly by a hyperplane; then, the hyperplane is expressed as Eq. 3.$$begin{aligned} {textbf{w}}^T {textbf{x}} + b = 0, end{aligned}$$
(3)
where ({textbf{w}}) is a weight vector of the hyperplane, and b is a bias. The weight vector is updated by minimizing Eq. 4.$$begin{aligned} L = {1 over 2} {textbf{w}}^T {textbf{w}} text { subject to } y_i ({textbf{w}}^T {textbf{x}} + b) ge 1 end{aligned}$$
(4)
We can transform Eqs. 4 to  5 by using the Lagrange multiplier method.$$begin{aligned} L^* = {1 over 2} {textbf{w}}^T {textbf{w}} – sum _{i=1}^n a_i { y_i ({textbf{w}}^T x_i + {-}) – 1 }, end{aligned}$$
(5)
where (a_i) is the Lagrange multiplier. We can efficiently solve Eq. 5 using a dual form. Furthermore, Eq. 5 can be solved in a case where it is not completely separated using a slack variable and a kernel trick can be used to estimate the nonlinear hyperplane.SORT trackerSORT, one of the frameworks for solving the multiple object tracking (MOT) problem, aims to achieve efficient real-time tracking24. The SORT method framework is created by combining the estimation step and the association step. The estimation step forecasts the next position of each predicted Daphnia magna. The association step matches the forecasting position and next true position of each predicted Daphnia magna. In the estimation step, the SORT framework uses the Kalman filter to forecast the position of the predicted Daphnia magna in the next frame. The position of each predicted Daphnia magna is expressed as Eq. 6.$$begin{aligned} {textbf{x}} = [u,v,s,r,{dot{u}}, {dot{v}}, {dot{s}}]^T end{aligned}$$
(6)
where u and v are the center positions of each predicted Daphnia magna, s is the scale size of the bounding box, and r is the aspect ratio of the bounding box. ({dot{u}}), ({dot{v}}), and ({dot{s}}) are the amounts of change in each variable. In the association step, to associate the forecasting position and true position, the framework adopts the intersection-over-union (IOU)36 as the association metric. The Hungarian algorithm is loaded into the SORT framework to perform fast and efficient Daphnia magna association prediction. In this study, a mixed metric of IOU36 and Euclidean distance37 was used instead of only the IOU that is used in SORT (Eq. 7) for more efficient association.$$begin{aligned} C_{ij} = (1-lambda ) {max_d – d_{ij} over max_d} + lambda cdot IOU_{ij} end{aligned}$$
(7)
where (d_{ij}) is the Euclidean distance between the i-th predicted Daphnia magna in the before frame and the j-th predicted Daphnia magna in the next frame, and (lambda ) is the weight of (IOU_{ij}). (IOU_{ij}) is the IOU between the i-th predicted Daphnia magna in the before-frame and the j-th predicted Daphnia magna in the next frame.MetricsThe binary confusion matrix consists of true positive (TP), true negative (TN), false positive (FP), and false negative (FN)38. TP is the number of cases where the predicted Daphnia magna matches the actual Daphnia magna, TN is the number of cases where the objects predicted as noise are actual noise, FP is the number of cases where the predicted Daphnia magna differs from the actual Daphnia magna, and FN is the number of cases where the objects predicted as noise are not actual noise. In this study, accuracy, recall, precision, and F1 scores (Eq. 8) were used as the metrics for comparing the machine learning methods.$$begin{aligned} begin{aligned} Accuracy&= {TP + FP over TP + TN + FP + FN} \ Recall&= {TP over TP + TN} \ Precision&= {TP over TP + FP} \ F1 score&= 2 times {Precision times Recall over Precision + Recall} end{aligned} end{aligned}$$
(8)
Standard MOT metrics to evaluate tracking performance include multi-object tracking accuracy (MOTA) and multi-object tracking precision (MOTP). An important task of MOT is to identify and track the same object across two frames. Identification (ID) precision (IDP), ID recall (IDR), ID F1 measure (IDF1), and ID switches (IDs) may be used as measures for evaluating the identification and tracking of the same objects39,40.Data analysisThe toxicity test using Daphnia magna was performed following the Korean official Acute Toxicity Test Method29. The test medium was prepared by dissolving KCl (8 mg/L), (text {MgSO}_4) (120 mg/L), (text {CaSO}_4 cdot 2 text {H}_2 text {O} ) (120 mg/L), and (text {NaHCO}_3) (192 mg/L) in deionized water. Considering that Daphnia magna are neonates for less than 24 h after birth29, five neonates were exposed to 50 mL of different concentrations of heavy metals such as Potassium dichromate, Copper(II) sulfate pentahydrate, and Lead(II) sulfate (6.25, 12.5, 25, 50, and 100%) and 50 mL of culture media. Potassium dichromate is a common inorganic reagent used as an oxidizing agent in chemical industries. Copper(II) sulfate pentahydrate is a trace material widely used in industrial processes and agriculture. A significant amount of copper is emitted in semiconductor manufacturing processes, which adversely impacts the aquatic ecosystem. When present as an ion in water, copper can be acutely toxic to aquatic organisms such as Daphnia magna. Lead(II) sulfate is another nonessential and nonbiodegradable heavy metal. It is highly toxic to numerous organisms even at low concentrations and can accumulate in aquatic ecosystems41. Twenty Daphnia magna (four replicates of five each) were exposed to each test solution for 24 h. The term “immobility” means that the Daphnia magna remains stationary after exposure to chemicals such as Potassium dichromate, Copper(II) sulfate pentahydrate, and Lead(II) sulfate. In this study, immobility was used as an endpoint identifier, and the number of mobile Daphnia magna were counted to evaluate the EC50 values for the samples using the ToxCalc 5.0 program (Tidepoll Software, USA).The locomotory responses of Daphnia magna were tested after 0, 12, 18, and 24 h of exposure at different concentrations. Potassium dichromate ((text {K}_2text {Cr}_2text {O}_7)) at 2 mg/L was connected to the Daphnia magna tracking system, and standard toxic substances were automatically diluted to 100%, 50%, 25%, 12.5%, and 6.25%. The automatic high-throughput Daphnia magna tracking system automatically measured the tracking results of a 1-minute-long video at hourly intervals. The average moving distance for 20 s of each Daphnia magna in each chamber was analyzed using a repeated measures ANOVA (RMANOVA). RMANOVA was used for the analysis of data obtained by repeatedly measuring the same Daphnia magna42. It analyzes the concentration effect excluding the time effect at each hour. The time effect means the change in average distance per 20 s. RMANOVA was implemented using the agricolae package of the R 4.0.4 program43. To remove the noise affecting RMANOVA, the Daphnia magna that remained stationary for 20 s or more were removed from the observations. In this study, we used the significance level at 5%. More

A hippo swims in Colombia’s Magdalena River, near where Pablo Escobar’s compound was located.Credit: Fernando Vergara/AP/Shutterstock

Colombian environment minister Susana Muhamad has triggered fear among researchers that she will protect, rather than reduce, a growing population of invasive hippos that threaten the country’s natural ecosystems and biodiversity. Although she did not directly mention the hippos — a contentious issue in Colombia — Muhamad said during a speech in late January that her ministry would create policies that prioritize animal well-being, including the creation of a new division of animal protection.
Landmark Colombian bird study repeated to right colonial-era wrongs
The hippos escaped from drug-cartel leader Pablo Escobar’s estate after he died in 1993. Left alone, the male and three females that Escobar had illegally imported from a US zoo established themselves in Colombia’s Magdalena River and some small lakes nearby — part of the country’s main watershed. After years of breeding, the ‘cocaine hippos’ have multiplied to about 150 individuals, scientists estimate.Given that the hippos (Hippopotamus amphibius) — considered the largest invasive animal in the world — have no natural predators in Colombia and have been mating at a steady rate, their population could reach 1,500 in 16 years, according to a modelling study published in 20211. “I do not understand what the government is waiting for to act,” says Nataly Castelblanco Martínez, a Colombian conservation biologist at the Autonomous University of Quintana Roo in Chetumal, Mexico, and co-author of the study. “If we don’t do anything, 20 years from now the problem will have no solution.”Researchers have called for a strict management plan that would eventually reduce the wild population to zero, through a combination of culling some animals and capturing others, then relocating them to facilities such as zoos. But the subject of what to do with the hippos has polarized the country, with some enamoured by the animals’ charisma and value as a tourist attraction and others concerned about the threat they pose to the environment and local fishing communities.‘A bit surreal’Several studies and observations suggest how destructive it could be to allow the Colombian hippo population to explode. A 2019 paper2, for example, showed that, compared with lakes without hippos, those where the animals have taken up residence contain more nutrients and organic matter that favour the growth of cyanobacteria — aquatic microbes associated with toxic algal blooms. These blooms can reduce water quality and cause mass fish deaths, affecting local fishing communities.

A sign near Doradal, Colombia, warns passersby of the danger of invasive hippos.Credit: Juancho Torres/Anadolu Agency via Getty

Other scientists have predicted that the hippos could displace endangered species that are native to the Magdalena River, such as the Antillean manatee (Trichechus manatus manatus), by outcompeting them for food and space. They caution that traffic accidents and attacks on people caused by the hippos will become more common. And they warn that wildlife traffickers are already taking advantage of the situation by illegally selling baby hippos — a trend that could intensify.“It’s a bit surreal,” says Jorge Moreno Bernal, a vertebrate palaeontologist at the University of the North in Barranquilla, Colombia. “This is just a taste of what may come.”When Colombian authorities first recognized the speed at which the hippo population was growing, during the 2000s, they acted to reduce their numbers. But in 2009, when photos appeared online after soldiers gunned down Pepe, Escobar’s fugitive male hippo, the outcry from animal-rights activists and others plunged the environment ministry into an “institutional paralysis”, says Sebastián Restrepo Calle, an ecologist at Javeriana University in Bogotá.Researchers say that the hippos don’t belong in Colombia — they are native to sub-Saharan Africa. Simulations run by Castelblanco Martínez and her colleagues suggest that to reduce the population to zero by 2033, about 30 hippos would need to be removed from the wild population per year1. No other course of action, including sterilization or castration, would eradicate them, according to the modelling of various management scenarios, says Castelblanco Martínez.The cost of inactionThe worry now is that, instead of basing decisions on evidence and expertise in conservation, the government is listening to popular opinion, says Restrepo Calle. Neither Muhamad nor representatives of the environment ministry replied to Nature’s requests for comment.
Ancient stone tools suggest early humans dined on hippo
“Why prioritize one species over our own ecosystems?” — especially a species that isn’t native, asks Alejandra Echeverri, a Colombian conservation scientist at Stanford University in California. Along with her colleagues, Echeverri published a study last month showing that Colombia has few policies governing invasive species compared with its overall number of biodiversity policies3.Animals-rights advocates, meanwhile, argue that they aren’t ignoring environmental concerns. Luis Domingo Gómez Maldonado, an animal-rights activist and specialist in animal law at Saint Thomas University in Bogotá, says “It’s not about saving the hippos on a whim,” but rather about solving the issue while also giving the hippos justice. “My indisputable position is: let’s save as many individuals as possible, let’s do it ethically.”Researchers, too, say they have the animals’ best interests at heart. “Even if [advocates] don’t see it, we care about the hippos,” Castelblanco Martínez says. “The more time that passes, the more hippos will either have to be culled, castrated or captured.”The question is whether environmental authorities will act swiftly to draft and enforce a management plan that is both ethical and effective. Should they sit on the issue for too long, Castelblanco Martínez warns, rural communities that are most affected by the hippos might take matters into their own hands.If the government doesn’t cull them, she says, people will use shotguns to do it. More