Philippa Kaur

1.
Kroeker, K. J. et al. Impacts of ocean acidification on marine organisms: quantifying sensitivities and interaction with warming. Glob. Change Biol. 19, 1884–1896 (2013).
Article  Google Scholar 
2.
Pfister, C. A. et al. Historical baselines and the future of shell calcification for a foundation species in a changing ocean. Proc. Royal Soc. B-Biol. Sci. 283, https://doi.org/10.1098/rspb.2016.0392 (2016).

3.
Przeslawski, R., Byrne, M. & Mellin, C. A review and meta-analysis of the effects of multiple abiotic stressors on marine embryos and larvae. Glob. Change Biol. 21, 2122–2140 (2015).
Article  Google Scholar 

4.
Waldbusser, G. G. & Salisbury, J. E. Ocean acidification in the coastal zone from an organism’s perspective: multiple system parameters, frequency domains, and habitats. Annu. Rev. Mar. Sci. 6, 221–247 (2014).
Article  Google Scholar 

5.
Helmuth, B., Mieszkowska, N., Moore, P. & Hawkins, S. J. Living on the edge of two changing worlds: forecasting the responses of rocky intertidal ecosystems to climate change. Annu. Rev. Ecol. Evolution Syst. 37, 373–404 (2006).
Article  Google Scholar 

6.
Kordas, R. L., Donohue, I. & Harley, C. D. G. Herbivory enables marine communities to resist warming. Sci. Adv. 3, https://doi.org/10.1126/sciadv.1701349 (2017).

7.
Kroeker, K. J., Kordas, R. L. & Harley, C. D. G. Embracing interactions in ocean acidification research: confronting multiple stressor scenarios and context dependence. Biol. Lett. 13, https://doi.org/10.1098/rsbl.2016.0802 (2017).

8.
Jellison, B. M. & Gaylord, B. Shifts in seawater chemistry disrupt trophic links within a simple shoreline food web. Oecologia 190, 955–967 (2019).
Article  Google Scholar 

9.
Wootton, J. T., Pfister, C. A. & Forester, J. D. Dynamic patterns and ecological impacts of declining ocean pH in a high-resolution multi-year dataset. Proc. Natl Acad. Sci. USA 105, 18848–18853 (2008).
CAS  Article  Google Scholar 

10.
Burrows, M. T. et al. Global-scale species distributions predict temperature-related changes in species composition of rocky shore communities in Britain. Glob. Change Biol. 26, 2093–2105 (2020).
Article  Google Scholar 

11.
Pershing, A. J. et al. Slow adaptation in the face of rapid warming leads to collapse of the Gulf of Maine cod fishery. Science 350, 809–812 (2015).
CAS  Article  Google Scholar 

12.
Wang, Z. H. A. et al. The marine inorganic carbon system along the Gulf of Mexico and Atlantic coasts of the United States: insights from a transregional coastal carbon study. Limnol. Oceanogr. 58, 325–342 (2013).
CAS  Article  Google Scholar 

13.
Colton, H. S. On some varieties of Thais lapillus in the Mount Desert region, a study of individual ecology. Proc. Acad. Natl Sci. Phila. 68, 440–454 (1916).
Google Scholar 

14.
Fisher, J. A. D. Exploring ecology’s attic: overlooked ideas on intertidal food webs. Bull. Ecol. Soc. Am. 86, 145–151 (2005).

15.
Connell, J. H. Effects of competition, predation by Thais lapillus, and other factors on natural populations of barnacle Balanus balanoides. Ecol. Monogr. 31, 61–104 (1961).
Article  Google Scholar 

16.
Lubchenco, J. Plant species diversity in a marine intertidal community—importance of herbivore food preference and algal competitive abilities. Am. Naturalist 112, 23–39 (1978).
Article  Google Scholar 

17.
Lubchenco, J. & Menge, B. A. Community development and persistence in a low rocky intertidal zone. Ecol. Monogr. 48, 67–94 (1978).
Article  Google Scholar 

18.
Petraitis, P. S. & Dudgeon, S. R. Variation in recruitment and the establishment of alternative community states. Ecology 96, 3186–3196 (2015).
CAS  Article  Google Scholar 

19.
Lord, J. P., Harper, E. M. & Barry, J. P. Ocean acidification may alter predator-prey relationships and weaken nonlethal interactions between gastropods and crabs. Mar. Ecol. Prog. Ser. 616, 83–94 (2019).
CAS  Article  Google Scholar 

20.
Sorte, C. J. B. et al. Long-term declines in an intertidal foundation species parallel shifts in community composition. Glob. Change Biol. 23, 341–352 (2017).
Article  Google Scholar 

21.
Salisbury, J. E. & Jönsson, B. F. Rapid warming and salinity changes in the Gulf of Maine alter surface ocean carbonate parameters and hide ocean acidification. Biogeochemistry 141, 401–418 (2018).
CAS  Article  Google Scholar 

22.
Southward, A. J. 40 years of changes in species composition and population density of barnacles on a rocky shore near Plymouth. J. Mar. Biol. Assoc. UK 71, 495–513 (1991).
Article  Google Scholar 

23.
Hawkins, S. J., Southward, A. J. & Genner, M. J. Detection of environmental change in a marine ecosystem—evidence from the western English Channel. Sci. Total Environ. 310, 245–256 (2003).
CAS  Article  Google Scholar 

24.
Hawkins, S. J. et al. Consequences of climate-driven biodiversity changes for ecosystem functioning of North European rocky shores. Mar. Ecol. Prog. Ser. 396, 245–259 (2009).
Article  Google Scholar 

25.
Mieszkowska, N., Burrows, M. T., Pannacciulli, F. G. & Hawkins, S. J. Multidecadal signals within co-occurring intertidal barnacles Semibalanus balanoides and Chthamalus spp. linked to the Atlantic Multidecadal Oscillation. J. Mar. Syst. 133, 70–76 (2014).
Article  Google Scholar 

26.
Findlay, H. S., Kendall, M. A., Spicer, J. I. & Widdicombe, S. Future high CO2 in the intertidal may compromise adult barnacle Semibalanus balanoides survival and embryonic development rate. Mar. Ecol. Prog. Ser. 389, 193–202 (2009).
Article  Google Scholar 

27.
Thomas, A. C., Townsend, D. W. & Weatherbee, R. Satellite-measured phytoplankton variability in the Gulf of Maine. Continental Shelf Res. 23, 971–989 (2003).
Article  Google Scholar 

28.
Commito, J. A., Jones, B. R., Jones, M. A., Winders, S. E. & Como, S. After the fall: legacy effects of biogenic structure on wind-generated ecosystem processes following mussel bed collapse. Diversity-Basel 11, https://doi.org/10.3390/d11010011 (2019).

29.
Beal, B. F. et al. Spatial variability in recruitment of an infaunal bivalve: experimental effects of predator exclusion on the softshell clam (Mya arenaria L.) along three tidal estuaries in southern Maine, USA. J. Shellfish Res. 37, 27 (2018).
Article  Google Scholar 

30.
Roman, J. Diluting the founder effect: cryptic invasions expand a marine invader’s range. Proc. R. Soc. B-Biol. Sci. 273, 2453–2459 (2006).
Article  Google Scholar 

31.
Vadas, R. L.,Sr. & Elner, R. W. Plant-animal interactions in the north-west Atlantic. In Plant-Animal Interactions in the Marine Benthos (eds. John, D. M., Hawkins, S. J. & Price, J. H.). Systematics Association Special Vol. 46, 33–60 (Oxford: Clarendon Press, 1992).

32.
Wethey, D. S. et al. Response of intertidal populations to climate: effects of extreme events versus long term change. J. Exp. Mar. Biol. Ecol. 400, 132–144 (2011).
Article  Google Scholar 

33.
Maggs, C. A. et al. Evaluating signatures of glacial refugia for North Atlantic benthic marine taxa. Ecology 89, S108–S122 (2008).
Article  Google Scholar 

34.
Wares, J. P. & Cunningham, C. W. Phylogeography and historical ecology of the North Atlantic intertidal. Evolution 55, 2455–2469 (2001).
CAS  Article  Google Scholar 

35.
Chapman, J. W., Carlton, J. T., Bellinger, M. R. & Blakeslee, A. M. H. Premature refutation of a human-mediated marine species introduction: the case history of the marine snail Littorina littorea in the Northwestern Atlantic. Biol. Invasions 9, 995–1008 (2007).
Article  Google Scholar 

36.
Cloern, J. E. & Jassby, A. D. Drivers of changes in estuarine-coast ecosystems: discoveries from four decades of study in San Francisco Bay. Rev. Geophys. 50, 33 (2012).
Article  Google Scholar 

37.
Kemp, W. M. et al. Eutrophication of Chesapeake Bay: historical trends and ecological interactions. Mar. Ecol. Prog. Ser. 303, 1–29 (2005).
Article  Google Scholar 

38.
Kübler, J. E. & Dudgeon, S. R. Prediciting effects of ocean acidification and warming on algae lacking carbon concentrating mechanisms. PLoS ONE 10, https://doi.org/10.1371/journal.pone.0132806 (2015).

39.
Petraitis, P. S., Liu, H. & Rhile, E. C. Densities and cover data for intertidal organisms in the Gulf of Maine, USA, from 2003 to 2007. Ecology 89, 588 http://esapubs.org/archive/ecol/E089/032 (2008).
Article  Google Scholar 

40.
Petraitis, P. S., Liu, H. & Rhile, E. C. Barnacle, fucoid, and mussel recruitment in the Gulf of Maine, USA, from 1997 to 2007. Ecology 90, 571 http://esapubs.org/archive/ecol/E090/039/ (2009).
Article  Google Scholar 

41.
Petraitis, P. S. & Vidargas, N. A. Marine intertidal organisms found in experimental clearings on sheltered shores in the Gulf of Maine, USA. Ecology 87, 796 http://esapubs.org/archive/ecol/E087/047/ (2006).
Article  Google Scholar 

42.
Lewis, J. R. The Ecology of Rocky Shores (English University Press, 1964).

43.
Horton, T. et al. World register of marine species (WoRMS). http://www.marinespecies.org (2010).

44.
Petraitis, P. S. & Dudgeon, S. R. Data and R scripts for analyses of declines in invertebrate species from the Gulf of Maine, USA, 1997–2018 ver 1. Environmental Data Initiative https://doi.org/10.6073/pasta/7ad4e0c85ff585199cb2a007eb291c9d (2020).

45.
Hadfield, J. D. MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package. J. Stat. Softw. 33, 1–22 (2010).
Article  Google Scholar 

46.
Burnham, K. P. & Anderson, D. R. Model Selection and Multi-Model Inference: A Practical Information-Theoretic Approach, Second edn. (Springer, 2002).

47.
Bartón, K. MuMIn: multi-model inference R package version 1.40.0. https://CRAN.R-project.org/package=MuMIn (2019).

48.
Honaker, J., King, G. & Blackwell, M. Amelia II: a program for missing data. J. Stat. Softw. 45, 1–47 (2011).
Article  Google Scholar 

49.
Choirat, C., Honaker, J., Imai, K., King, G. & Lau, O. Zelig: everyone’s statistical software. version 5.1.6.1. http://zeligproject.org (2008).

50.
Imai, K., King, G. & Lau, O. Toward a common framework for statistical analysis and development. J. Computational Graph. Stat. 17, 892–913 (2008).
Article  Google Scholar  More

The goal of this study was to examine the fluid flows directly adjacent to propulsor surfaces in order to better understand how metachronal propulsors interact with fluids for thrust production. Based on the direct comparison of the mean contributions of pulling vs. pushing forces throughout the power stroke of replicate individual propulsors (Fig. 4, which are generated by negative vs. positive pressure fields, respectively) we suggest that the propulsors of the animals examined rely predominantly on negative pressure for generating thrust. The assertion that these propulsor level observations to apply to the movement of the whole animal requires the assumptions that, first, the propulsors we quantified are representative of the all the other propulsors contributing to swimming thrust, and second, that the thrust generated for the whole animal is due to accumulated total thrust generated by individual propulsors. Although our data addresses the first of these assumptions by replicating individual propulsors, we cannot document the second assumption that the total thrust represents the sum of all individual propulsor elements. While this second assumption is intuitively appealing, our data is confined to the small spatial and temporal scales around individual propulsor elements. Confirmation that the whole-organism thrust results from of the summation of individual contributions requires experiments at different scales than those used in the current study.
Thrust generated by a propulsor is ultimately determined by the overall pressure gradient across the propulsor. So does it matter whether that gradient is dominated by negative or positive pressure? We believe that this distinction is fundamental for understanding why animal propulsors bend in a surprisingly characteristic and narrow range. Rigid paddle designs are dominated by positive pressure pushing against a fluid, which in turn, generates thrust pushing a body forward. Bending at propulsor margins encourages vortex formation on the lee side of the propulsor (Figs. 2, 3) that differs from rigid propulsors. Counter-rotating vortices formed on the lee side of a bending propulsor accelerate fluid at the intersection of the vortices12,15. The fluid thus accelerated relative to the leading edge of the propulsor is the basis of the pressure gradient across the propulsor surface. In turn, this elevated pressure gradient generates high thrust and is the reason for the dominant contribution of suction thrust to natural bending propulsors. More generally, negative pressure fields are a fundamental feature of vortices which are universally formed around objects moving in fluids (except at the lowest Reynolds numbers). Lift, a different propulsive mode that relies on negative pressure, is a well-known example that illustrates how kinematics and morphology can enhance negative pressure for thrust2. Lift occurs when a foil separates flow traveling over and under the foil surface. With the correct foil shape and kinematics, the separation of flow can generate strong negative pressure fields above the foil leading to an upward pulling thrust on the foil. This lift relies on the negative pressure field and foil shape.
To be clear, the thrust generated by limbs and ctenes in this study is not lift because the forces generated by lift are directed perpendicular to the direction of flow and the forces we describe are oriented in the direction of flow (Fig. 3). However, like lift, we suggest that paddles must move with prescribed kinematics to generate enhanced negative pressure fields. Bending kinematics in particular have been shown to greatly enhance vorticity and along with that, negative pressure11,16,17. Rigid, non-bending paddles generate different hydrodynamic structures than we observed13,14 and do not generate strong negative pressure fields16,17,18,19. Therefore, the kinematics of bending appear to be important for generating strong negative pressure fields around moving propulsors.
Until recently, technical constraints have limited our ability to investigate the scope of the benefits of using negative pressure for thrust. However several numerical studies, and a few experimental studies, have compared rigid to flexible propulsors. These studies have demonstrated that first, bending enhances negative pressure fields, second, bending generates elevated thrust, and third, bending enhances hydrodynamic efficiency12,16,17,20,21,22,23,24. The hydrodynamic patterns around bending propulsors show that negative pressure fields associated with bends generate significantly greater flow velocities than positive pressure fields (Fig. 2e11,12,16,21). This would lead to enhanced momentum transfer and explain the enhanced thrust observed for bending propulsors. The similar bending kinematics between the limbs and ctenes in this study and the swimming and flying animals from Lucas et al.1 suggests that these small paddling swimmers may employ similar hydrodynamic features as flying birds and swimming fish. If these bending patterns are predominately used to generate negative pressure fields for thrust, it follows that there is a need for greater focus on negative pressure around bending propulsors in order to understand the extent of the benefits of animals experience by pulling rather than pushing themselves through fluids.
Despite the vast difference in scale and Reynolds number, the results of this study suggest that the small metachronal paddles of swimming invertebrates may produce some similar effects as flapping wings in birds and insects. For example, there are similarities in the degree of bending and location of bending for the paddles in this study and the spanwise flexibility of birds and insects1. Such spanwise flexibility was found to be beneficial and yielded an increase in thrust coefficient, and a small decrease in power-input requirement, resulting in higher efficiency25.
In addition to the benefits for single propulsors, negative pressure fields can facilitate the movement and coordination of multiple propulsors which have antiplectic metachronal wave kinematics. During an antiplectic metachronal wave, a leading propulsor will begin the power stroke and, after it has initiated its stroke, the propulsor immediately behind it will initiate its own power stroke. This sequential pattern will continue for all the subsequent propulsors in the antiplectic wave. The predominately negative pressure on the leeward of each propulsor can serve to facilitate the kinematics of the adjacent propulsor by reducing the hydrodynamic resistance necessary to initiate and complete its power stroke26. In addition, the negative pressure in the gap between adjacent propulsors can serve as a cue for the adjacent propulsor to initiate its power stroke. It has been suggested that the ctenes of ctenophores require such cues to coordinate the metachronal kinematics26,27,28. At lower Reynolds numbers (Re  More

1.
IPCC. Climate change 2013: the Physical Science Basis. Contribution of working group I to the fifth assessment report of the intergovernmental panel on climate change. In: Stocker TF, Qin D, Plattner G-K, Tignor M, Allen SK, Boschung J, et al. editors. Cambridge, United Kingdom and New York, USA: Cambridge University Press; 2013. p. 1535.
2.
De Azevedo TR, Junior CC, Junior AB, dos Santos Cremer M, Piatto M, Tsai DS, et al. SEEG initiative estimates of Brazilian greenhouse gas emissions from 1970 to 2015. Sci Data. 2018;5:180045.
PubMed  PubMed Central  Article  CAS  Google Scholar 

3.
West TAP, Börner J, Fearnside PM. Climatic benefits from the 2006-2017 avoided deforestation in Amazonian Brazil. Front For Glob Change. 2019;2:52.
Article  Google Scholar 

4.
Malhi Y, Roberts JT, Betts RA, Killeen TJ, Li W, Nobre CA. Climate change, deforestation, and the fate of the Amazon. Science. 2008;319:169–73.
CAS  PubMed  Article  Google Scholar 

5.
Curtis PG, Slay CM, Harris NL, Tyukavina A, Hansen MC. Classifying drivers of global forest loss. Science. 2018;361:1108–11.
CAS  PubMed  Article  Google Scholar 

6.
Paula FS, Rodrigues JLM, Zhou J, Wu L, Mueller RC, Mirza BS, et al. Land use change alters functional gene diversity, composition and abundance in Amazon forest soil microbial communities. Mol Ecol. 2014. https://doi.org/10.1111/mec.12786.

7.
Evans PN, Boyd JA, Leu AO, Woodcroft BJ, Parks DH, Hugenholtz P, et al. An evolving view of methane metabolism in the Archaea. Nat Rev Microbiol. 2019;17:219–32.
CAS  PubMed  Article  Google Scholar 

8.
Serrano-Silva N, Sarria-Guzmán Y, Dendooven L, Luna-Guido M. Methanogenesis and methanotrophy in soil: a review. Pedosphere. 2014;24:291–307.
CAS  Article  Google Scholar 

9.
Zinder SH. Physiological ecology of methanogens. In: Ferry JG, editor. Methanogenesis: ecology, physiology, biochemistry and genetics. New York, NY, USA: Chapman & Hall; 1993. p. 128–206.

10.
Schink B, Stams AJM. Syntrophism among prokaryotes. In: Rosenberg E, DeLong F, Lory S, Stackebrandt E, Thompson F, editors. The prokaryotes: prokaryotic communities and ecophysiology. Berlin: Springer; 2013. p. 471–93.

11.
Costa KC, Leigh JA. Metabolic versatility in methanogens. Curr Opin Biotechnol. 2014;29:70–5.
CAS  PubMed  Article  Google Scholar 

12.
Semrau JD, DiSpirito AA, Yoon S. Methanotrophs and copper. FEMS Microbiol Rev. 2010;34:496–531.
CAS  PubMed  Article  Google Scholar 

13.
van Teeseling MCF, Pol A, Harhangi HR, van der Zwart S, Jetten MSM, Op den Camp HJM, et al. Expanding the verrucomicrobial methanotrophic world: description of three novel species of Methylacidimicrobium gen. nov. Appl Environ Microbiol. 2014;80:6782–92.
PubMed  PubMed Central  Article  CAS  Google Scholar 

14.
Zheng Y, Hou L, Chen F, Zhou J, Liu M, Yin G, et al. Denitrifying anaerobic methane oxidation in intertidal marsh soils: Occurrence and environmental significance. Geoderma. 2020;357:113943.
CAS  Article  Google Scholar 

15.
Zheng Y, Wang H, Liu Y, Zhu B, Li J, Yang Y, et al. Methane-dependent mineral reduction by aerobic methanotrophs under hypoxia. Environ Sci Technol Lett. 2020;7:606–12.
CAS  Article  Google Scholar 

16.
Murrell JC, Dalton H. Nitrogen fixation in obligate methanotrophs. J Gen Microbiol. 1983;129:3481–6.
CAS  Google Scholar 

17.
Graham DW, Chaudhary JA, Hanson RS, Arnold RG. Factors affecting competition between type I and type II methanotrophs in two-organism, continuous-flow reactors. Microb Ecol. 1993;25:1–17.
CAS  PubMed  Article  Google Scholar 

18.
Amaral JA, Knowles R. Growth of methanotrophs in methane and oxygen counter gradients. FEMS Microbiol Lett. 1995;126:215–20.
CAS  Article  Google Scholar 

19.
Dunfield PF, Liesack W, Henckel T, Knowles R, Conrad R. High-affinity methane oxidation by a soil enrichment culture containing a type II methanotroph. Appl Environ Microbiol. 1999;65:1009–14.
CAS  PubMed  PubMed Central  Article  Google Scholar 

20.
Bull ID, Parekh NR, Hall GH, Ineson P, Evershed RP. Detection and classification of atmospheric methane oxidizing bacteria in soil. Nature. 2000;405:175–8.
CAS  PubMed  Article  Google Scholar 

21.
Meyer KM, Klein AM, Rodrigues JLM, Nüsslein K, Tringe SG, Mirza BS, et al. Conversion of Amazon rainforest to agriculture alters community traits of methane-cycling organisms. Mol Ecol. 2017;26:1547–56.
CAS  PubMed  Article  Google Scholar 

22.
Baani M, Liesack W. Two isoenzymes of particulate methane monooxygenase with different methane oxidation kinetics are found in Methylocystis sp. strain SC2. Proc Natl Acad Sci USA. 2008;105:10203–8.
CAS  PubMed  Article  Google Scholar 

23.
Sharp CE, Smirnova AV, Graham JM, Stott MB, Khadka R, Moore TR, et al. Distribution and diversity of Verrucomicrobia methanotrophs in geothermal and acidic environments. Environ Microbiol. 2014;16:1867–78.
CAS  PubMed  Article  Google Scholar 

24.
Keller M, Goreau TJ, Wofsy SC, Kaplan WA, McElroy MB. Production of nitrous oxide and consumption of methane by forest soils. Geophys Res Lett. 1986;10:1156–9.
Article  Google Scholar 

25.
Steudler PA, Melillo JM, Feigl BJ, Neill C, Piccolo MC, Cerri CC. Consequence of forest-to-pasture conversion on CH4 fluxes in the Brazilian Amazon Basin. J Geophys Res Atmos. 1996;101:18547–54.
CAS  Article  Google Scholar 

26.
Fearnside PM. Greenhouse gases from deforestation in Brazilian Amazonia: net committed emissions. Climatic Change. 1997;35:321–60.
CAS  Article  Google Scholar 

27.
Fernandes SAP, Bernoux M, Cerri CC, Feigl BJ, Piccolo MC. Seasonal variation of soil chemical properties and CO2 and CH4 fluxes in unfertilized and P-fertilized pastures in an Ultisol of the Brazilian Amazon. Geoderma. 2002;107:227–41.
CAS  Article  Google Scholar 

28.
Wick B, Veldkamp E, de Mello WZ, Keller M, Crill P. Nitrous oxide fluxes and nitrogen cycling along a pasture chronosequence in Central Amazonia, Brazil. Biogeosciences. 2005;2:175–87.
CAS  Article  Google Scholar 

29.
Tveit AT, Urich T, Svenning MM. Metatranscriptomic analysis of active peat soil microbiome. Appl Environ Microbiol. 2014;80:5761–72.
CAS  PubMed  PubMed Central  Article  Google Scholar 

30.
Singer E, Wagner M, Woyke T. Capturing the genetic makeup of the active microbiome in situ. ISME J. 2017;11:1949–63.
CAS  PubMed  PubMed Central  Article  Google Scholar 

31.
Keiblinger KM, Wilhartitz IC, Schneider T, Roschitzki B, Schmid E, Eberl L, et al. Soil metaproteomics—comparative evaluation of protein extraction protocols. Soil Biol Biochem. 2012;54:14–24.
CAS  PubMed  PubMed Central  Article  Google Scholar 

32.
Kroeger ME, Nüsslein K. Stable isotope probing—detection of active microbes in soil. In: Elsas JDV, editor. Modern soil microbiology, 3rd ed. Boca Raton: CRC Press; 2019.

33.
Esson K, Lin X, Kumaresan D, Chanton JP, Murrell JC, Kostka JE. Alpha- and gammaproteobacterial methanotrophs codominate the active methane-oxidizing communities in an acidic boreal peat bog. Appl Environ Microbiol. 2016;82:2363–71.
CAS  PubMed  PubMed Central  Article  Google Scholar 

34.
Dumont MG, Radajewski SM, Miguez CB, McDonald IR, Murrell JC. Identification of a complete methane monooxygenase operon from soil by combining stable isotope probing and metagenomic analysis. Environ Microbiol. 2006;8:1240–50.
CAS  PubMed  Article  Google Scholar 

35.
Lu Y, Leuders T, Friedrich MW, Conrad R. Detecting active methanogenic populations on rice roots using stable isotope probing. Environ Microbiol. 2005;7:326–36.
CAS  PubMed  Article  Google Scholar 

36.
Neufeld JD, Vohra J, Dumont MG, Lueders T, Manefield M, Friedrich MW, et al. DNA stable-isotope probing. Nat Protoc. 2007;2:860–6.
CAS  PubMed  Article  Google Scholar 

37.
Holmes AJ, Costello A, Lidstrom ME, Murrell JC. Evidence that particulate methane monooxygenase and ammonia monooxygenase may be evolutionarily related. FEMS Microbiol Lett. 1995;132:203–8.
CAS  PubMed  Article  Google Scholar 

38.
Costello A, Lidstrom ME. Molecular characterization of functional and phylogenetic genes from natural populations of methanotrophs in lake sediments. Appl Environ Microbiol 1999;65:5066–74.
CAS  PubMed  PubMed Central  Article  Google Scholar 

39.
Steinberg LM, Regan JM. Phylogenetic comparison of the methanogenic communities from an acidic, oligotrophic fen and an anaerobic digester treating municipal wastewater sludge. Appl Environ Microbiol. 2008;74:6663–71.
CAS  PubMed  PubMed Central  Article  Google Scholar 

40.
Meyer KM, Morris AH, Webster K, Klein AM, Kroeger ME, Meredith LK, et al. Belowground changes to community structure alter methane-cycling dynamics in Amazonia. 2020. https://www.biorxiv.org/content/10.1101/2020.03.10.984807v1. in review.

41.
Navarrete AA, Tsai SM, Mendes LW, Faust K, de Hollander M, Cassman NA, et al. Soil microbiome responses to the short-term effects of Amazonian deforestation. Mol Ecol. 2015;24:2433–48.
CAS  PubMed  Article  Google Scholar 

42.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.
CAS  PubMed  PubMed Central  Article  Google Scholar 

43.
Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome. 2018;6:90.
PubMed  PubMed Central  Article  Google Scholar 

44.
Keegan KP, Glass EM, Meyer F. MG-RAST, a metagenomics service for analysis of microbial community structure and function. Methods Mol Biol. 2016;1399:207–33.
CAS  PubMed  Article  Google Scholar 

45.
R Core Team. R: a language and environment for statistical computing. Vienna Austria: R Foundation for Statistical Computing; 2018. https://www.r-project.org/.

46.
Oksanen JF, Blanchet G, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: community ecology package. R package version. 2018;2:4–6. https://CRAN.R-project.org/package=vegan.

47.
Parks DH, Tyson GW, Hugenholtz P, Beiko RG. STAMP: statistical analysis of taxonomic and functional profiles. Bioinformatics. 2014;30:3123–4.
CAS  PubMed  PubMed Central  Article  Google Scholar 

48.
Wickham H. ggplot2: elegant graphics for data analysis. New York: Spring-Verlag; 2016.
Google Scholar 

49.
Parks DH, Tyson GW, Hugenholtz P, Beiko RG. STAMP: statistical analysis of taxonomic and functional profiles. Bioinformatics. 2014;30:3123–4.
CAS  PubMed  PubMed Central  Article  Google Scholar 

50.
Pinto AJ, Raskin L. PCR biases distort bacterial and archaeal community structure in pyrosequencing datasets. PLoS ONE. 2012;7:e43093.
CAS  PubMed  PubMed Central  Article  Google Scholar 

51.
Huse SM, Welch DM, Morrison HG, Sogin ML. Ironing out the wrinkles in the rare biosphere through improved OTU clustering. Environ Microbiol. 2010;12:1889–98.
CAS  PubMed  PubMed Central  Article  Google Scholar 

52.
Quince C, Walker AW, Simpson JT, Loman NJ, Segata N. Shotgun metagenomics, from sampling to sequencing and analysis. Nat Biotechnol. 2017;35:833–44.
CAS  PubMed  Article  Google Scholar 

53.
Coyotzi S, Pratscher J, Murrell JC, Neufeld JD. Targeted metagenomics of active microbial populations with stable-isotope probing. Curr Opin Biotechnol. 2016;41:1–8.
CAS  PubMed  Article  Google Scholar 

54.
Eyice Ö, Namura M, Chen Y, Mead A, Samavedam S, Schäfer H. SIP metagenomics identifies uncultivated Methylophilaceae as dimethylsulphide degrading bacteria in soil and lake sediment. Int Soc Microb Ecol J. 2015;9:2336.
CAS  Google Scholar 

55.
Cardoso D, Sarkinen T, Alexander S, Amorim AM, Bittrich V, Celis M, et al. Amazon plant diversity revealed by a taxonomically verified species list. Proc Natl Acad Sci USA. 2017;114:10695–700.
CAS  PubMed  Article  Google Scholar 

56.
Sgouridis F, Ullah S. Soil greenhouse gas fluxes, environmental controls, and the partitioning of N2O sources in UK natural and seminatural land use types. J Geophys Res Biogeosciences. 2017. https://doi.org/10.1002/2017JG003783.

57.
Wanyama I, Pelster DE, Butterbach-Bahl K, Verchot LV, Martius C, Rufino MC. Soil carbon dioxide and methane fluxes from forests and other land use types in an African tropical montane region. Biogeochemistry. 2019;143:171–90.
CAS  Article  Google Scholar 

58.
Lackner N, Hintersonnleitner A, Wagner AO, Illmer P. Hydrogenotrophic methanogenesis and autotrophic growth of Methanosarcina thermophila. Archaea. 2018;2018:4712608.
PubMed  PubMed Central  Article  CAS  Google Scholar 

59.
Ferry JG. Fundamentals of methanogenic pathways that are key to the biomethanation of complex biomass. Curr Opin Biotechnol. 2011;22:351–7.
CAS  PubMed  PubMed Central  Article  Google Scholar 

60.
Welander PV, Metcalf WW. Loss of the mtr operon in Methanosarcina blocks growth on methanol, but not methanogenesis, and reveals an unknown methanogenic pathway. Proc Natl Acad Sci USA. 2005;102:10664–9.
CAS  PubMed  Article  Google Scholar 

61.
Mand TD, Metcalf WW. Energy conservation and hydrogenase function in methanogenic archaea, in particular the genus Methanosarcina. Microbiol Mol Biol Rev. 2019;83:e00020–19.
CAS  PubMed  PubMed Central  Article  Google Scholar 

62.
Lammel DR, Feigl BJ, Cerri CC, Nusslein K. Specific microbial gene abundances and soil parameters contribute to C, N, and greenhouse gas process rates after land use change in Southern Amazonian Soils. Front Microbiol. 2015;6:1057.
PubMed  PubMed Central  Article  Google Scholar 

63.
Kroeger ME, Delmont TO, Meyer KM, Guo J, Khan K, Rodrigues JLM, et al. New biological insights into how deforestation in Amazonia affects soil microbial communities using metagenomics and metagenome-assembled genomes. Front Microbiol. 2018;9:1635.
PubMed  PubMed Central  Article  Google Scholar 

64.
Pedrinho A, Mendes LW, Merloti LF, de Cassia da Fonseca M, de Souza Cannavan F, Tsai SM. Forest-to-pasture conversion and recovery based on assessment of microbial communities in Eastern Amazon rainforest. FEMS Microbiol Ecol. 2019;95:fiy236.
CAS  PubMed  Article  Google Scholar 

65.
Carini P, Marsden PJ, Leff JW, Morgan EE, Strickland MS, Fierer N. Relic DNA is abundant in soil and obscures estimates of soil microbial diversity. Nat Microbiol. 2017;2:16242.
CAS  Article  Google Scholar 

66.
Lennon JT, Muscarella ME, Placella SA, Lehmkuhl BK. How, when, and where relic DNA affects microbial diversity. MBio. 2018;9:e00637–18.
PubMed  PubMed Central  Google Scholar 

67.
Sirois SH, Buckley DH. Factors governing extracellular DNA degradation dynamics in soil. Environ Microbiol Rep. 2019;11:173–84.
CAS  PubMed  Article  Google Scholar 

68.
Mancinelli RL. The regulation of methane oxidation in soil. Annu Rev Microbiol. 1995;49:581–605.
CAS  PubMed  Article  Google Scholar 

69.
Guggenheim C, Brand A, Bürgmann H, Sigg L, Wehrli B. Aerobic methane oxidation under copper scarcity in a stratified lake. Sci Rep. 2019;9:4817.
PubMed  PubMed Central  Article  CAS  Google Scholar 

70.
Knapp CW, Fowle DA, Kulczycki E, Roberts JA, Graham DW. Methane monooxygenase gene expression mediated by methanobactin in the presence of mineral copper sources. Proc Natl Acad Sci USA. 2007;104:12040–5.
CAS  PubMed  Article  Google Scholar 

71.
Rodrigues JLM, Pellizari VH, Mueller R, Baek K, da C Jesus E, Paula FS, et al. Conversion of the Amazon rainforest to agriculture results in biotic homogenization of soil bacterial communities. Proc Natl Acad Sci USA. 2013;110:988–93.
CAS  PubMed  Article  Google Scholar 

72.
Hoehler T, Gunsalus RP, McInerney MJ. Environmental Constraints that Limit Methanogenesis. In Timmis KN, editor. Handbook of Hydrocarbon and Lipid Microbiology. Berlin, Heidelberg: Springer; 2010. p. 635–54.

73.
Elshahed MS, Bhupathiraju VK, Wofford NQ, Nanny MA, McInerney MJ. Metabolism of benzoate, cyclohex-1-ene carboxylate, and cyclohexane carboxylate by “Syntrophus aciditrophicus” strain SB in syntrophic association with H2-using microorganisms. Appl Environ Microbiol. 2001;67:1728–38.
CAS  PubMed  PubMed Central  Article  Google Scholar 

74.
Zehnder AJB, Stumm W. Geochemistry and biogeochemistry of anaerobic habitats. In Zehnder AJB editor. Biology of Anaerobic Microorganisms. New York, NY, USA: John Wiley & Sons; 1988. p. 1–38.

75.
Wolfe RS. My kind of biology. Annu Rev Microbiol. 1991;45:1–35.
CAS  PubMed  Article  Google Scholar 

76.
Strittmatter AW, Liesegang H, Rabus R, Decker I, Amann J, Andres S, et al. Genome sequence of Desulfobacterium autotrophicum HRM2, a marine sulfate reducer oxidizing organic carbon completely to carbon dioxide. Environ Microbiol. 2009;11:1038–55.
CAS  PubMed  PubMed Central  Article  Google Scholar 

77.
Badziong W, Ditter B, Thauer RK. Acetate and carbon dioxide assimilation by Desulfovibrio vulgaris (Marburg), growing on hydrogen and sulfate as sole energy source. Arch Microbiol. 1979;123:301–5.
CAS  Article  Google Scholar 

78.
Watson SW, Waterbury JB. Characteristics of two marine nitrite oxidizing bacteria, Nitrospina gracilis nov. gen. nov. sp. and Nitrococcus mobilis nov. gen. nov. sp. Arch Microbiol. 1971;77:203–30.
Google Scholar 

79.
Sorokin DY, Tourova TP, Lysenko AM, Mityushina LL, Kuenen JG. Thioalkalivibrio thiocyanoxidans sp. nov. and Thioalkalivibrio paradoxus sp. nov., novel alkaliphilic, obligately autotrophic, sulfur-oxidizing bacteria capable growth on thiocyanate, from soda lakes. IJSEM. 2002;52:657–64.
CAS  PubMed  Google Scholar 

80.
Kukla J, Whitfeld T, Cajthaml T, Baldrian P, Veselá-Šimáčková H, Novotný V, et al. The effect of traditional slash-and-burn agriculture on soil organic matter, nutrient content, and microbiota in tropical ecosystems of Papua New Guinea. Land Degrad Dev. 2018;30:166–77.
Article  Google Scholar 

81.
Banat IM, Marchant R. Geobacillus activities in soil and oil contamination remediation. In: Logan NA, De Vos P, editors. Endospore-forming soil bacteria, soil biology 27, Berlin Heidelberg: Springer-Verlag; 2011. p. 259–70.

82.
Durrer A, Margenot AJ, Silva LCR, Bohannan BJM, Nüsslein K, van Haren J, et al. Beyond total carbon: long-term effect of deforestation on Amazonian soils. Biogeochemistry BIOG-D-19-00139. in review.

83.
Nazina TN, Tourova TP, Poltaraus AB, Novikova EV, Grigoryan AA, Ivanova AE, et al. Taxonomic study of aerobic thermophilic bacilli: descriptions of Geobacillus subterraneus gen. nov., sp. nov. and Geobacillus uzenensis sp. nov. from petroleum reservoirs and transfer of Bacillus thermoleovorans, Bacillus kaustophilus, Bacillus thermodenitrificans to Geobacillus as the new combinations G. stearothermophilus, G.th. IJSEM. 2001;51:433–46.
CAS  PubMed  Google Scholar 

84.
Schnürer A, Schink B, Svensson BH. Clostridium ultunense sp. nov., a mesophilic bacterium oxidizing acetate in syntrophic association with a hydrogenotrophic methanogenic bacterium. IJSEM. 1996;46:1145–52.
Google Scholar 

85.
Myhr S, Torsvik T. Denitrovibrio acetiphilus, a novel genus and species of dissimilatory nitrate-reducing bacterium isolated from an oil reservoir model column. IJSEM. 2000;50:1611–9.
CAS  PubMed  Google Scholar 

86.
Yao H, Conrad R, Wassmann R, Neue HU. Effect of soil characteristics on sequential reduction and methane production in sixteen rice paddy soils from Chine, the Philippines, and Italy. Biogeochemistry. 1999;47:269–95.
CAS  Article  Google Scholar 

87.
Achtnich C, Bak F, Conrad R. Competition for electron donors among nitrate reducers, ferric iron reducers, sulfate reducers, and methanogens in anoxic paddy soil. Biol Fertil Soils. 1995;19:65–72.
CAS  Article  Google Scholar 

88.
Verchot LV, Davidson EA, Cattânio JH, Ackerman IL. Land-use change and biogeochemical controls of methane fluxes in soils of Eastern Amazonia. Ecosystems. 2000;3:41–56.
CAS  Article  Google Scholar 

89.
Wagner R, Zona D, Oechel W, Lipson D. Microbial community structure and soil pH correspond to methane production in Arctic Alaska soils. Environ Microbiol. 2017;19:3398–410.
CAS  PubMed  Article  Google Scholar 

90.
Feigl BJ, Melillo J, Cerri CC. Changes in the origin and quality of soil organic matter after pasture introduction in Rôndonia (Brazil). Plant Soil. 1995;175:21–9.
CAS  Article  Google Scholar 

91.
Salimon CI, Davidson EA, Victoria RL, Melo AWF. CO2 flux from soil in pastures and forests in southwestern Amazonia. Glob Change Biol. 2004;10:833–43.
Article  Google Scholar 

92.
Wilhelm RC, Singh R, Eltis LD, Mohn WW. Bacterial contributions to delignification and lignocellulose degradation in forest soils with metagenomic and quantitative stable isotope probing. ISME J. 2019;13:413–8.
CAS  PubMed  Article  Google Scholar 

93.
Pepe-Ranney C, Campbell AN, Koechli CN, Berthrong S, Buckley DH. Unearthing the ecology of soil microorganisms using a high-resolution DNA-SIP approach to explore cellulose and xylose metabolism in soil. Front Microbiol. 2016;7:703.
PubMed  PubMed Central  Article  Google Scholar 

94.
Murase J, Hordijk K, Tayasu I, Bodelier PLE. Strain-specific incorporation of methanotrophic biomass into eukaryotic grazers in a rice field soil revealed by PLFA-SIP. FEMS Microbiol Ecol. 2011;75:284–90.
CAS  PubMed  Article  Google Scholar 

95.
Radajewski S, Ineson P, Parekh NR, Murrell JC. Stable-isotope probing as a tool in microbial ecology. Nature. 2000;403:646–9.
CAS  PubMed  Article  Google Scholar 

96.
Cébron A, Bodrossy L, Chen Y, Singer AC, Thompson IP, Prosser JI, et al. Identity of active methanotrophs in landfill cover soil as revealed by DNA-stable isotope probing. FEMS Microbiol Ecol. 2007;62:12–23.
PubMed  Article  CAS  Google Scholar 

97.
Mathew RP, Feng Y, Githinji L, Ankumah R, Balkcom KS. Impact of no-tillage and conventional tillage systems on soil microbial communities. Appl Environ Soil Sci. 2012;2012:548620.
Article  Google Scholar 

98.
Eilers KG, Debenport S, Anderson S, Fierer N. Digging deeper to find unique microbial communities: the strong effect of depth on the structure of bacterial and archaeal communities in soil. Soil Biol Biochem. 2012;50:58–65.
CAS  Article  Google Scholar 

99.
Zhou J, Wu L, Deng Y, Zhi X, Jiang Y-H, Tu Q, et al. Reproducibility and quantitation of amplicon sequencing-based detection. ISME J. 2011;5:1303–13.
CAS  PubMed  PubMed Central  Article  Google Scholar 

100.
Wen C, Wu L, Qin Y, Van Nostrand JD, Ning D, Sun B, et al. Evaluation of the reproducibility of amplicon sequencing with Illumina MiSeq platform. PLoS ONE. 2017. https://doi.org/10.1371/journal.pone.0176716.

101.
Berenguer E, Gardner TA, Ferreira J, Aragao LE, Camargo PB, Cerri CE, et al. Developing cost-effective field assessments of carbon stocks in human-modified tropical forests. PLoS ONE. 2015;10:e0133139.
PubMed  PubMed Central  Article  CAS  Google Scholar  More

Global food demand will rise substantially over the coming decades. Meeting this demand while decreasing the environmental footprint of agriculture is one of largest challenges of the twenty-first century1,2,3. A growing world population and changing diets are projected to double4 meat and dairy consumption between 2000 and 2050. As one of the main feed sources for livestock, grasslands play a key role in meeting this demand. With over 33 million km2, permanent grasslands account for ~ 25% of the world’s land cover. Over two thirds of this area is utilised for agriculture, making it the most dominant land use5. Sustainably increasing grassland productivity is therefore crucial to ensure future global food security6,7.
Phosphorus (P) is an essential nutrient, often limiting plant growth8. P fertilisation is therefore needed to sustain productivity in agricultural systems across the world. Because the world’s P reserves are decreasing, the importance of judicious P use will increase over the coming century. Although estimates of global P reserves vary, the costs of high quality P fertilisers will increase, as will the global demand for these fertilisers9,10,11,12. Differences in climate, geography, agricultural development, and fertilisation practices have led to great global imbalances of P in agricultural land13,14,15,16. In parts of Europe, North America, and China, historical applications of manure and fertilisers have resulted in positive P balances and increased risk of eutrophication of surface waters17. In many other regions, predominantly in tropical areas, farmers struggle to maintain soil P availability to sustain optimal rates of crop production18. Recent predictions suggest that global P inputs in grasslands will have to increase fourfold to support an 80% increase in grass yield projected for 205015, which implies an urgent need to increase use efficiency of P fertiliser sources.
The large diversity in agronomic P status of soils across the world and the projected increase in cost and demand of P fertilisers necessitate a rethink of the use of P resources: are we applying fertilisers at the right rates to the right soils? The success of fertiliser application depends on conditions created by climate and management19,20 and is strongly governed by soil properties such as pH and concentrations of metal oxides and Ca in soil that can impact P availability to plants8,21,22. However, data for these relationships are fragmentary and country- or region-specific, and global assessments are lacking23,24. Here we use a meta-analysis on a global database of 67 studies and 1227 observations with a wide range of soil properties and climatic conditions to assess the general effect of P fertilisation on grassland production across the world. Furthermore, we identify soil-related driving factors that determine the success of fertiliser applications.
Our dataset included data from field grasslands all over the world (Supplementary Fig. 1). Most studies originated from Europe and North America, but due to several studies with many observations from the Australian continent, there were almost as many observations from Oceania. We analysed our dataset using two different metrics: the response ratio (RR) as measure for the relative increase in dry matter yield as a result of P fertilisation, and P agronomic efficiency (PAE) expressing the absolute yield increase per unit of P applied.
Factors controlling the success of phosphorus fertilisation
P fertilisation increased grassland yield by 37% (95% confidence interval: 33 to 40%; Fig. 1; Supplementary Table 3) averaged over all grasslands, soil types, and fertility levels, resulting in a PAE of 32 kg kg−1 (Fig. 2; Supplementary Table 4). In other words, dry matter yields increased by 32 kg per kg of P applied on average. Yield responses to P additions increased with P application rates: rates below 25 kg P ha−1 increased yields by 40% on average, whereas applying over 100 kg P ha−1 increased grass yield by 65%. An exception to this pattern were grasslands fertilised with 25–50 kg P ha−1, which responded to a lesser extent than those in other categories. This is likely an artefact due to a relatively high average soil P status of studies included in this category (Supplementary Fig. 3), which may have led to high yields in the control treatments. The PAE, on the other hand, decreased with P application rates (Supplementary Table 4): yields increased by 53 kg per kg P applied at rates lower than 25 kg P ha−1, but only by 12 kg kg−1 P at rates higher than 100 kg P ha−1. This indicates that finding a balance between P input and yield response is crucial for optimising fertiliser effectivity, as the agronomic efficiency decreases with higher application rates.
Figure 1

Impact of phosphorus (P) fertilisation for the controlling factors crop, P rate, climate, and mean annual temperature expressed as relative yield increase per category. The 95% confidence intervals are represented by the error bars, and the number of studies and observations per category are between parentheses; *,**,***Significant controlling factor effect at an α of 0.05, 0.01 and 0.001, respectively.

Full size image

Figure 2

The Phosphorus agronomic efficiency (PAE) for different controlling factors per subgroup. The effect is expressed for crop, climate, and P status (Olsen-equivalent) × P rate (c). Low SPT: ≤ 10 mg P kg−1; high SPT:  > 10 mg P kg−1; low rate: P rate ≤ 50 kg P ha−1; high rate: P rate  > 50 kg P ha−1 The 95% confidence intervals are represented by the error bars, and the number of studies and observations per category are between parentheses; *, **,***Significant controlling factor effect at an α of 0.05, 0.01 and 0.001, respectively.

Full size image

Systems that included legumes responded more strongly to P fertilisation than systems without legumes (Fig. 1). On average, P fertiliser increased yield in grass/legume systems by 54%, but only by 25% in grassland systems without legumes. These numbers corresponded with a PAE of 46 kg kg−1 for grass/legume and 22 kg kg−1 for grass-only systems, meaning that P fertilisation was roughly twice as effective in grasslands with legumes than in those without legumes. Legumes like alfalfa and clover are regularly included in grassland mixtures, mainly because they provide extra N inputs to the plant-soil system by establishing a symbiosis with N-fixing microorganisms23. These results likely reflect that legumes generally require more P than grasses, and can acquire it less easily due to thicker roots and shorter root hairs11,25,26.
In our database, more than half (36) of the studies included more than one N treatment. Overall, the N application rate had little effect on the response of grasslands to P fertilisation. There was no significant effect of N rate on the PAE (Supplementary Table 4). Yield responses to P fertiliser at N application rates over 200 kg N ha−1 were slightly but significantly smaller than at lower N rates (Supplementary Table 3). However, if N limitation of the grasslands would have played a prominent role, a general increase in response to P fertiliser with increasing N rate would have been observed. These results suggest that differences in yield responses were mainly driven by a response to P fertilisation rather than to N fertilisation.
Geographical variation in responses
P application increased grassland yields in tropical regions (i.e. latitudes ≤ 35°) significantly more strongly than in temperate grasslands (Fig. 1, Supplementary Table 3). However, because yields of tropical grasslands were relatively low, the PAE of fertiliser application did not differ significantly between the two regions (34 and 31 kg kg−1 for tropical and temperate regions, respectively; Fig. 2, Supplementary Table 4). These results likely reflect that soils in (sub)tropical regions are often highly weathered, nutrient-poor, and have a low P availability due to high abundancy of adsorbents like Al and Fe oxides8. In contrast, decades of manure and fertiliser applications have resulted in a build-up of soil P levels well beyond crop requirements and a corresponding decrease in yield response to P fertiliser application17,27 in many temperate regions (e.g. North America, Europe, and New Zealand). The differences in response of temperate and tropical grasslands are also reflected in the results for mean annual temperature (MAT; Fig. 1, Supplementary Table 3), with grasslands in colder regions (MAT  20 °C reacting the strongest. Higher temperatures may lead to more rapid plant production and to an increase in mineralisation of organic matter. Correlation analysis of the controlling factors showed that MAT and latitude among our studies were strongly correlated (Supplementary Fig. 4; Spearman’s ρ = -0.95).
Yield responses to P fertilisation were significantly smaller in Asia, North America, and Europe (+ 15 to + 29%) than in South America, Oceania, and Africa (+ 58 to + 94%). The PAE ranged from 12 kg kg−1 for studies in Asia to 74 kg kg−1 for studies in Oceania and even 117 kg kg−1 for the one African study included in our dataset (Supplementary Table 4). The continents with grasslands that showed a strong response to P fertilisation roughly coincide with the areas that have relatively low P inputs and outputs, as modelled by Sattari et al.15. Taken together, these results imply that Africa and Oceania with low P inputs responded strongly to P fertilisation whereas grasslands in Europe, North America and Asia with relatively high P inputs over the past decades, showed a weak response to P fertilisation.
Do we apply phosphorus fertilisers to the right soils?
We used various soil parameters as controlling factors (Table 1) to identify what soil properties drive differences in yield response to P fertilisation. One of the most important parameters is the agronomic P status of the soil, which is commonly determined with a soil P test (SPT). Because soil type, climate, and crop response vary considerably across the world, each country and sometimes even region has its own SPT method and classification system28,29. Given this large variety of SPT procedures (and resulting P concentrations) in use, we applied conversion formulas published in peer-reviewed papers to express reported SPT values in our database as ‘Olsen-equivalent’ P values wherever possible (see Supplementary Methods).
Table 1 Controlling factors and categories distinguished in the meta-analysis.
Full size table

Grasslands on soils with low SPT values (≤ 5 mg P kg−1) responded strongest to P fertilisation with a yield increase of 110% on average (Fig. 3, Supplementary Table 3). Conversely, P additions to soils with SPT values  > 5 mg P kg−1 increased yields by 7–25%. Although yield response decreased dramatically with increasing SPT values, the responses at relatively high SPT values (10–25 and  > 25 mg P kg−1) were still statistically significant. Critical values (that is, SPT levels for which the yield is 95% of the maximum yield) for grass of 23–25 mg kg−1 Olsen P have been reported previously for English grasslands30, which coincides with the limited yield response for soils in the highest SPT category. A study of 25 Spanish soils also showed an average critical SPT of 24 mg Olsen P kg−1 for ryegrass, although there was a wide spread for individual soils, ranging from 11 to 46 mg kg−131. For a range of Australian grassland species, however, lower critical SPT values (between 9 and 15 mg kg−1) have been determined32. This variety of critical SPT values found in literature illustrates that the effect of P fertilisation is strongly dependent on soil, climate, and even grassland species. Therefore, our results here do not give a hard SPT limit beyond which further P applications are rendered ineffective, but do indicate a strong decrease in effectiveness at higher SPT values.
Figure 3

Effect of different soil characteristics on the impact of phosphorus (P) fertilisation. The effect is expressed for soil P status based on Olsen P-equivalent, soil pH, organic matter content, and clay content. The 95% confidence intervals are represented by the error bars, and the number of studies and observations per category are between parentheses; *, **,***Significant controlling factor effect at an α of 0.05, 0.01 and 0.001, respectively.

Full size image

The strong yield response to P fertilisation on soils with low SPT values was not merely the result of a low yield of the control treatments. PAE was also highest (75 kg kg−1) for soils with SPT ≤ 5 mg P kg−1 (Supplementary Table 4) and fertilisation on these soils was 3 to 8 times as effective as on soils with higher SPT values in terms of absolute yield increases. Without correcting for the P application rate, absolute yield responses (average yield of treated plots minus average yield of control plots) to P fertilisation varied substantially (− 2.7 to 11.3 tonnes ha−1; Supplementary Fig. 5). The largest response (on average 2.7 tonnes ha−1 increase) and variation to P fertilisation were found for soils in the lowest SPT category. The yield response decreased with higher SPT (Supplementary Fig. 5). Figure 2 shows that both relatively low (≤ 50 kg P ha−1) and relatively high ( > 50 kg P ha−1) P application rates on soils with a low P status (≤ 10 mg P kg−1 Olsen-equivalent) were more effective than any P fertilisation rate on soils with a relatively high P status ( > 10 mg P kg−1 Olsen-equivalent). The high PAE of large application rates on soils with a low P status (Fig. 2) may be the result of the binding behaviour of P in soil: in soils with a low P status (where relatively more P adsorption occurs), relatively high P inputs are required to raise the level of plant-available P, so grassland on these soils will benefit relatively more from high application rates. Conversely, applying large amounts of P to soils with a relatively high P status ( > 10 mg P kg−1) showed a low PAE.
Yield responses to P applications were highest on grasslands with a soil pH of 5–6 (60% yield increase; Fig. 3) whereas lower and higher pH levels resulted in lower (11–26%) yield responses. We observed the same pattern for PAE, where studies with a pH of 5 to 6 had a 50 kg yield increase per kg of P fertilised, whereas for soils with a pH above 7 this was only 11 kg (Supplementary Table 4). Soil pH is a crucial parameter in determining the availability of P to crops8. In acidic mineral soils, binding of P to Fe and Al (hydr)oxides is often the main factor that governs the level of plant available P. In contrast, in soils with pH values above 7, P is more likely to form poorly soluble Ca-P precipitates, decreasing plant available P. The relative availability of soil P is highest at soil pH levels of 5 to 733,34, which would imply that around this pH fertiliser P application would yield the strongest responses.
We found a positive correlation between the soil organic matter (OM) content and yield response to P fertilisation (Fig. 3). On average, P application increased yield by only 11% on soils with an OM content below 2% (PAE was 7.2 kg kg−1 on average and this effect was not statistically significant). Yield responses were much higher (41–80%) in soils with an OM content of  > 5%. The PAE was 9 times as high in soils with  > 5% OM as in soils with  More

1.
Pierce, R. H. & Henry, M. S. Harmful algal toxins of the Florida red tide (Karenia brevis): Natural chemical stressors in South Florida coastal ecosystems. Ecotoxicology 17, 623–631 (2008).
CAS  PubMed  PubMed Central  Article  Google Scholar 
2.
Fauquier, D. A. et al. Brevetoxicosis in seabirds naturally exposed to Karenia brevis blooms along the central west coast of Florida. J. Wildl. Dis. 49, 246–260 (2013).
PubMed  Article  Google Scholar 

3.
Fire, S. E. et al. Brevetoxin-associated mass mortality event of bottlenose dolphins and manatees along the east coast of Florida USA. Mar. Ecol. Prog. Ser. 526, 241–251 (2015).
ADS  Article  CAS  Google Scholar 

4.
Foley, A. M. et al. Assessing Karenia brevis red tide as a mortality factor of sea turtles in Florida USA. Dis. Aquat. Organ. 132, 109–124 (2019).
CAS  PubMed  Article  Google Scholar 

5.
Gannon, D. P. et al. Effects of Karenia brevis harmful algal blooms on nearshore fish communities in southwest Florida. Mar. Ecol. Prog. Ser. 378, 171–186 (2009).
ADS  CAS  Article  Google Scholar 

6.
Backer, L. C. et al. Occupational exposure to aerosolized brevetoxins during Florida red tide events: Effects on a healthy worker population. Environ. Health Perspect. 113, 644–649 (2005).
CAS  PubMed  PubMed Central  Article  Google Scholar 

7.
Fleming, L. E. et al. Aerosolized red-tide toxins (Brevetoxins) and asthma. Chest 131, 187–194 (2007).
CAS  PubMed  PubMed Central  Article  Google Scholar 

8.
Hoagland, P., Anderson, D. M., Kaoru, Y. & White, A. W. The economic effects of harmful algal blooms in the United States: Estimates, assessment issues, and information needs. Estuaries 25, 819–837 (2002).
Article  Google Scholar 

9.
Gilbert, P. M., Anderson, D. M., Genstien, P., Granéli, E. & Sellner, K. G. The global, complex phenomena of harmful algal blooms. Oceanography 18, 136–147 (2005).
Article  Google Scholar 

10.
Hallegraeff, G. M. Harmful algal blooms: A global overview. in Manual on Harmful Marine Microalgae. Monographs on Oceanographic Methodology 25–49 (IOC-UNE-SCO, 2003).

11.
Vargo, G. A. A brief summary of the physiology and ecology of Karenia brevis Davis (G. Hansen and Moestrup comb. nov.) red tides on the West Florida Shelf and of hypotheses posed for their initiation, growth, maintenance, and termination. Harmful Algae 8, 573–584 (2009).
CAS  Article  Google Scholar 

12.
Heil, C. A. et al. Blooms of Karenia brevis (Davis) G. Hansen & Ø. Moestrup on the West Florida shelf: Nutrient sources and potential management strategies based on a multi-year regional study. Harmful Algae 38, 127–140 (2014).
Article  Google Scholar 

13.
Stumpf, R. P. et al. Monitoring Karenia brevis blooms in the Gulf of Mexico using satellite ocean color imagery and other data. Harmful Algae 2, 147–160 (2003).
CAS  Article  Google Scholar 

14.
Flagherty, K. E. & Landsberg, J. H. Effects of a persistent red tide (Karenia brevis) bloom on community structure and species-specific relative abundance of nekton in a Gulf of Mexico estuary. Estuaries Coasts 34, 417–439 (2011).
Article  CAS  Google Scholar 

15.
McHugh, K. A., Allen, J. B., Barleycorn, A. A. & Wells, R. S. Severe Karenia brevis red tides influence juvenile bottlenose dolphin (Tursiops truncatus) behavior in Sarasota Bay Florida. Mar. Mammal Sci. 27, 622–643 (2011).
Article  Google Scholar 

16.
Deutsch, C. J. et al. Seasonal movements, migratory behavior, and site fidelity of West Indian manatees along the Atlantic coast of the United States. Wildl. Monogr. 2, 1–77 (2003).
Google Scholar 

17.
Gannon, J. G., Scolardi, K. M., Reynolds, J. E. III., Koelsch, J. K. & Kessenich, T. J. Habitat selection by manatees in Sarasota Bay Florida. Mar. Mammal Sci. 23, 133–143 (2007).
Article  Google Scholar 

18.
Garrott, R. A. et al. Trends in counts of Florida manatees at winter aggregation sites. J. Wildl. Manag. 58, 642 (1994).
Article  Google Scholar 

19.
Semeyn, D. J. et al. Aerial surveys of manatees (Trichechus manatus) in Lee County, Florida, provide insights regarding manatee abundance and real time information for managers and enforcement officers. J. Coast. Conserv. 15, 573–583 (2011).
Article  Google Scholar 

20.
Wells, R. S. et al. Bottlenose dolphins as marine ecosystem sentinels: Developing a health monitoring system. EcoHealth 1, 246–254 (2004).
Article  Google Scholar 

21.
Fire, S. E. et al. Brevetoxin exposure in bottlenose dolphins (Tursiops truncatus) associated with Karenia brevis blooms in Sarasota Bay Florida. Mar. Biol. 152, 827–834 (2007).
CAS  Article  Google Scholar 

22.
Sadchatheeswaran, S., Belanger, M. & Wittnich, C. A comparison of published brevetoxin tissue levels in West Indian manatee, bottlenose dolphin and double-crested cormorants in southwest Florida. 5, 8 (2012).

23.
Monczak, A. et al. Sound patterns of snapping shrimp, fish, and dolphins in an estuarine soundscape of the southeastern USA. Mar. Ecol. Prog. Ser. 609, 49–68 (2019).
Article  Google Scholar 

24.
Walters, S., Lowerre-Barbieri, S., Bickford, J. & Mann, D. Using a passive acoustic survey to identify spotted seatrout spawning sites and associated habitat in Tampa Bay Florida. Trans. Am. Fish. Soc. 138, 88–98 (2009).
Article  Google Scholar 

25.
Mann, D. A. Remote sensing of fish using passive acoustic monitoring. Acoust. Today 8, 8 (2012).
Article  Google Scholar 

26.
Harris, S. A., Shears, N. T. & Radford, C. A. Ecoacoustic indices as proxies for biodiversity on temperate reefs. Methods Ecol. Evol. 7, 713–724 (2016).
Article  Google Scholar 

27.
Parks, S. E., Miksis-Olds, J. L. & Denes, S. L. Assessing marine ecosystem acoustic diversity across ocean basins. Ecol. Inform. 21, 81–88 (2014).
Article  Google Scholar 

28.
Parsons, M. J. G., Salgado-Kent, C. P., Marley, S. A., Gavrilov, A. N. & McCauley, R. D. Characterizing diversity and variation in fish choruses in Darwin Harbour. ICES J. Mar. Sci. J. Cons. 73, 2058–2074 (2016).
Article  Google Scholar 

29.
Hildebrand, J. Anthropogenic and natural sources of ambient noise in the ocean. Mar. Ecol. Prog. Ser. 395, 5–20 (2009).
ADS  Article  Google Scholar 

30.
Quick, J. A. & Henderson, G. E. Effects of Gymnodinium breve red tide on fishes and birds: A preliminary report on behavior, anatomy, hematology, and histopathology. Proc. Gulf Coast Region. Symp. Dis. Aquat. Anim. 2, 85–113 (1982).
Google Scholar 

31.
Landsberg, J. H. & Steidinger, K. A. A historical review of Gymnodinium breve red tides implicated in mass mortalities of the manatee (Trichechus manatus latirostris) in Florida, USA. in In: Reguera B, Blanco J, Fernandez ML, Wyatt T (eds) (UNESCO, 1998).

32.
Walters, S., Lowerre-Barbieri, S., Bickford, J., Tustison, J. & Landsberg, J. Effects of Karenia brevis red tide on the spatial distribution of spawning aggregations of sand seatrout Cynoscion arenarius in Tampa Bay Florida. Mar. Ecol. Prog. Ser. 479, 191–202 (2013).
ADS  Article  Google Scholar 

33.
FWC. Fish Kill Database. https://public.myfwc.com/fwri/FishKillReport.

34.
Baden, D. G. & Mende, T. Toxicity of two toxins from the Florida red tide marine dinoflagellate Ptychdiscus brevis. Toxicon 20, 457–461 (1982).
CAS  PubMed  Article  Google Scholar 

35.
Steidinger, K. A., Burklew, M. & Ingle, R. The effects of Gymnodinium breve toxin on estuarine animals. In Marine pharmacognosy: action of marine toxins at the cellular level 179–202 (Academic Press, New York, 1973).
Google Scholar 

36.
Landsberg, J. H. The effects of harmful algal blooms on aquatic organisms. Rev. Fish. Sci. 10, 113–390 (2002).
Article  Google Scholar 

37.
Fire, S. E. et al. Prevalence of brevetoxins in prey fish of bottlenose dolphins in Sarasota Bay Florida. Mar. Ecol. Prog. Ser. 368, 283–294 (2008).
ADS  Article  Google Scholar 

38.
Flewelling, L. J. et al. Red tides and marine mammal mortalities: Brevetoxicosis. Nature 435, 755–756 (2005).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

39.
Naar, J. P. et al. Brevetoxins, like ciguatoxins, are potent ichthyotoxic neurotoxins that accumulate in fish. Toxicon 50, 707–723 (2007).
CAS  PubMed  PubMed Central  Article  Google Scholar 

40.
The Manatee County Water Atlas. https://manatee.wateratlas.usf.edu/ (2020).

41.
Caldwell, M. C. & Caldwell, D. K. Individualized whistle contours in bottlenosed dolphins (Tursiops truncatus). Nature 207, 434–435 (1965).
ADS  Article  Google Scholar 

42.
O’Shea, T. J. & Poché, L. B. Aspects of underwater sound communication in Florida manatees (Trichechus manatus latirostris). J. Mammal. 87, 1061–1071 (2006).
Article  Google Scholar 

43.
Fisheries, N. 2018–2019 Bottlenose dolphin unusual mortality event Southwest Florida | NOAA Fisheries. /southeast/marine-life-distress/2018–2019-bottlenose-dolphin-unusual-mortality-event-southwest (2019).

44.
Red Tide Manatee Mortalities. Florida Fish and Wildlife Conservation Commission. https://myfwc.com/research/manatee/rescue-mortality-response/statistics/mortality/red-tide/.

45.
Montie, E. W. et al. Acoustic monitoring indicates a correlation between calling and spawning in captive spotted seatrout (Cynoscion nebulosus). PeerJ 5, e2944 (2017).
PubMed  PubMed Central  Article  Google Scholar 

46.
Myrberg, A. A., Ha, S. J. & Shamblott, M. J. The sounds of bicolor damselfish (Pomacentrus partitus): Predictors of body size and a spectral basis for individual recognition and assessment. J. Acoust. Soc. Am. 94, 3067–3070 (1993).
ADS  Article  Google Scholar 

47.
Tyson, R. B. & Wells, R. S. Sarasota Bay/Little Sarasota Bay Bottlenose Dolphin Abundance Estimates: 2015. 22 http://doi.org/10.7289/V5/RD-PRBD-2016-02 (2016).

48.
Lund, J., Kipling, C. & LeCren, E. The inverted microscope method of estimating algal numbers and the statistical basis of estimations of counting. Hydrobiology 2, 143–170 (1958).
Article  Google Scholar 

49.
Sournia, A. Phytoplankton Manual. in Monographs on oceanographic methodology vol. 6 (UNESCO, 1978).

50.
Sellner, K., Doucette, G. & Kirkpatrick, G. Harmful algal blooms: Causes, impacts, and detection. J. Ind. Microbiol. Biotechnol. 30, 383–406 (2003).
CAS  PubMed  Article  Google Scholar 

51.
Sprague, M. W. & Luczkovich, J. J. Measurement of an individual silver perch Bairdiella chrysoura sound pressure level in a field recording. J. Acoust. Soc. Am. 116, 3186–3191 (2004).
ADS  PubMed  Article  Google Scholar 

52.
Breder, C. M. Jr. Seasonal and diurnal occurrences of fish sounds in a small Florida bay. Am. Mus. Nat. Hist. 138, 325–378 (1968).
Google Scholar 

53.
Fish, M. P. & Mowbray, W. H. Sounds of Western North Atlantic fishes. A reference file of biological underwater sounds (Johns Hopkins Press, Baltimore, 1970).
Google Scholar 

54.
NOAA Solar Calculator. https://www.esrl.noaa.gov/gmd/grad/solcalc/.

55.
Wilke, C. O. cowplot: Streamlined plot theme and plot annotations for ‘ggplot2’. (2019).

56.
Grolemund & Wickham. lubridate. https://lubridate.tidyverse.org/ (2011).

57.
South, A. rnaturalearth: World Map Data from Natural Earth. https://CRAN.R-project.org/package=rnaturalearth (2017).

58.
Pebesma. sf. https://github.com/r-spatial/sf/ (2018).

59.
Wickham, H. et al. Welcome to the Tidyverse. J. Open Source Softw. 4, 2 (2019).
Article  Google Scholar 

60.
Garnier. viridis. https://cran.r-project.org/package=viridis (2018).

61.
Florida Fish and Wildlife Conservation Commission GIS & Mapping Data Downloads. https://geodata.myfwc.com/.

62.
R Core Team. https://www.R-project.org/ (2019).

63.
Bauer, D. F. Constructing confidence sets using rank statistics. J. Am. Stat. Assoc. 67, 687–690 (1972).
MATH  Article  Google Scholar 

64.
Hollander, M. & Wolfe, D. A. Nonparametric Statistical Methods 27–33 (John Wiley & Sons, New York, 1973).
Google Scholar 

65.
Hollander, M. & Wolfe, D. A. Nonparametric Statistical Methods 68–75 (John Wiley & Sons, New York, 1973).
Google Scholar  More

We estimated the abundance of desert bighorn sheep in a captive facility located within the Chihuahuan Desert of New Mexico (Fig. 1). The area is arid, mountainous, with steep cliffs punctuated by ravines. The entire facility is enclosed by a high fence, preventing desert bighorn sheep ingress or egress. These animals are wild, unconditioned to humans, and used by The New Mexico Game and Fish Department (NMDGF) to establish and bolster other desert bighorn sheep populations within New Mexico16.
NMDGF leads an annual spring census to enumerate desert bighorn sheep numbers within the facility. Methodologically, the census uses a ground crew of spaced individuals walking in a line perpendicular to pen fencing, within the facility (i.e., drive count or census). Each individual keeps track of neighboring individuals to space the line, and to count any desert bighorn sheep breaking past them. Most desert bighorn sheep are herded ahead of the line. Other biologists at high topographical sites use spotting-scopes and binoculars to count and age class the moving sheep. Census counts are classified by animal ages and sex. We consider young sheep any animal ≤ 1.5 years old. Adult rams and ewes consist of males and females aged  > 1.5 years old, respectively. The group “adult” includes all sheep  > 1.5 years old.
We established 11 motion activated camera traps (Bushnell Trophy Cam) at the centers of an 800-m grid with a random geographical start. Cameras were secured to T-posts or existing vegetation when rocky areas thwarted post establishment. Most cameras were oriented north, to minimize sun exposure in the imagery. When vegetation or rock blocked the camera view, the camera orientation was rotated eastward until a clear view was obtained. Camera heights were 0.9–1.2 m with declination perpendicular to the ground. Cameras were checked at 6 month intervals and the retrieved SD cards were never full. Cameras were motion activated and set at the shortest delay possible (10 s; meaning that the camera waits at least 10 s after recording a picture before it will record another). In practice, the fastest trigger speed that cameras recorded imagery was a mean of 14.92 s (N = 6 cameras; SD = 0.94), a value we rounded to 15 s and employed in our analyses. Cameras recorded one image per trigger. We deployed cameras by 15 May 2017 and retrieved cameras no earlier than 30 April 2018. Our analyses period began on 1 October 2017 and ended on 1 May 2018. We employed a 5 month acclimation period to avoid the cameras serving as a novel attractant for desert bighorn sheep, which would violate distance sampling assumptions. Further, we experienced 1 camera failure during this acclimation period, and relocated 2 misplaced cameras. Lastly, some desert bighorn sheep gathered in shady locations near cameras during hot months (June–September) which created extreme variation in the encounter rate. The 1 October 2017–1 May 2018 period lacked all of these issues.
Imagery of desert bighorn sheep were identified and then classified by sex and age class: rams, ewes, young, adults (an adult-sized animal with sex undiscernible), and unknown [undiscernible (e.g. picture of a hoof, or animals blocking a clear view another animal)17]. To quantify distances, an individual stood at each camera and used the printed images of desert bighorn sheep to position another person at the exact locations that an animal was imaged. Distances between individuals were measured with a laser rangefinder and metric tape (Fig. 2; the authors confirm that informed consent was obtained to publish the identifying information/images in an online open-access publication). We supplemented these measurements by recording distances between the camera and several recognizable objects in the images (e.g. rocks, plants), to ensure accurate distance delineations for the imaged sheep.
The sensitivity of a trail camera’s passive infrared sensor (PIR) will decline as radial distance from the camera increases. Other factors, like vegetation and topography, also cause animal detections to decline with distance. This situation makes distance sampling an appealing approach, as its’ foundational premise is that the probability of detecting an animal declines with distance from the observation point15. The technique relies on measured distances between animals and the observation point. Therefore, we used the distance measures of desert bighorn sheep to the respective trail camera imaging them, to estimate sheep abundances using the ‘distance’ package (version 0.9.812) in program R. We analyzed data using a fall period (1 October 2017–31 January 2018) and a spring period (1 March 2018–1 May 2018). Dates for these periods were selected by season while ensuring sufficient sample sizes.
The camera trap operates like a point count sample. The sampling angle for each camera’s field of view was 50°18. Therefore, we multiplied sampling effort by this fraction (50/360) to correctly represent the amount of sampling area.
When analyzing data within a distance sampling approach, obtaining unbiased density estimates relies on correspondence between the sampling period and the period describing when the target species are active and therefore available for detection10,15. In our situation, this means aligning the sampling period of distance measurements with the period of time that desert bighorn sheep were active and able to trigger the trail camera10. In our study, we acquired very few nocturnal images of desert bighorn sheep (Fig. 4). Therefore, we censored study effort and distance data to the period defined by 1 h before sunrise through 2 h after sunset for each month, on a daily basis. Doing so removed  More

1.
Schmitz, O. J., Miller, J. R. B., Trainor, A. M. & Abrahms, B. Toward a community ecology of landscapes: predicting multiple predator–prey interactions across geographic space. Ecology 98, 2281–2292 (2017).
PubMed  Article  Google Scholar 
2.
van Noordwijk, A. J. & de Jong, G. Acquisition and allocation of resources: their influence on variation in life history. Am. Nat. 128, 137–142 (1986).
Article  Google Scholar 

3.
Gaillard, J.-M. et al. Habitat-performance relationships: finding the right metric at a given spatial scale. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 365, 2255–2265 (2010).
PubMed  PubMed Central  Article  Google Scholar 

4.
Mosser, A., Fryxell, J. M., Eberly, L. & Packer, C. Serengeti real estate: Density vs. fitness-based indicators of lion habitat quality. Ecol. Lett.12, 1050–1060 (2009).

5.
Kosterman, M. K., Squires, J. R., Holbrook, J. D., Pletscher, D. H. & Hebblewhite, M. Forest structure provides the income for reproductive success in a southern population of Canada lynx. Ecol. Appl. 28, 1032–1043 (2018).
PubMed  Article  Google Scholar 

6.
DeCesare, N. J. et al. Linking habitat selection and predation risk to spatial variation in survival. J. Anim. Ecol. 83, 343–352 (2014).
PubMed  Article  Google Scholar 

7.
Hebblewhite, M., Merrill, E. H. & McDonald, T. L. Spatial decomposition of predation risk using resource selection functions: An example in a wolf-elk predator-prey system. Oikos 111, 101–111 (2005).
Article  Google Scholar 

8.
McLoughlin, P. D., Dunford, J. S. & Boutin, S. Relating predation mortality to broad-scale habitat selection. J. Anim. Ecol. 74, 701–707 (2005).
Article  Google Scholar 

9.
Ritchie, E. G. & Johnson, C. N. Predator interactions, mesopredator release and biodiversity conservation. Ecol. Lett. 12, 982–998 (2009).
PubMed  Article  Google Scholar 

10.
Vanak, A. T. et al. Moving to stay in place: behavioral mechanisms for coexistence of African large carnivores. Ecology 94, 2619–2631 (2013).
PubMed  Article  Google Scholar 

11.
Torretta, E., Serafini, M., Puopolo, F. & Schenone, L. Spatial and temporal adjustments allowing the coexistence among carnivores in Liguria (N–W Italy). Acta Ethol. 19, 123–132 (2016).
Article  Google Scholar 

12.
Preisser, E. L., Bolnick, D. I. & Benard, M. F. Scared to death? The effects of intimidation and consumption in predator–prey interactions. Ecology 86, 501–509 (2005).
Article  Google Scholar 

13.
Swanson, A. et al. Cheetahs and wild dogs show contrasting patterns of suppression by lions. J. Anim. Ecol. 83, 1418–1427 (2014).
PubMed  Article  Google Scholar 

14.
Levi, T. & Wilmers, C. Wolves—coyotes—foxes: a cascade among carnivores. Ecology 93, 921–929 (2012).
PubMed  Article  Google Scholar 

15.
Henke, S. E. & Bryant, F. C. Effects of coyote removal on the faunal community in western Texas. J. Wildl. Manag. 63, 1066 (1999).
Article  Google Scholar 

16.
Gehrt, S. D. & Prange, S. Interference competition between coyotes and raccoons: a test of the mesopredator release hypothesis. Behav. Ecol. 18, 204–214 (2007).
Article  Google Scholar 

17.
St-Pierre, C., Ouellet, J. P. & Crête, M. Do competitive intraguild interactions affect space and habitat use by small carnivores in a forested landscape?. Ecography (Cop.) 29, 487–496 (2006).
Article  Google Scholar 

18.
Shores, C. R., Dellinger, J. A., Newkirk, E. S., Kachel, S. M. & Wirsing, A. J. Mesopredators change temporal activity in response to a recolonizing apex predator. Behav. Ecol. 30, 1324–1335 (2019).
Article  Google Scholar 

19.
Allen, M. L., Elbroch, L. M., Wilmers, C. C. & Wittmer, H. U. The comparative effects of large carnivores on the acquisition of carrion by scavengers. Am. Nat. 185, 822–833 (2015).
PubMed  Article  Google Scholar 

20.
Janssen, A., Sabelis, M. W., Magalhães, S. & Van, T. Habitat structure affects intraguild predation. Ecology 88, 2713–2719 (2007).
PubMed  Article  Google Scholar 

21.
Finke, D. L. & Denno, R. F. Intraguild predation diminished in complex-structured vegetation: Implication for prey suppression. Ecology 83, 643–652 (2002).
Article  Google Scholar 

22.
Laurenson, M. K. High juvenile mortality in cheetahs (Acinonyx jubatus) and its consequences for maternal care. J. Zool. Soc. Lond. 234, 387–408 (1994).
Article  Google Scholar 

23.
Hunter, J. S., Durant, S. M. & Caro, T. M. To flee or not to flee: predator avoidance by cheetahs at kills. Behav. Ecol. Sociobiol. 61, 1033–1042 (2007).
Article  Google Scholar 

24.
Hilborn, A. et al. Cheetahs modify their prey handling behavior depending on risks from top predators. Behav. Ecol. Sociobiol.72, 74 (2018).

25.
Swanson, A., Arnold, T., Kosmala, M., Forester, J. & Packer, C. In the absence of a “landscape of fear”: How lions, hyenas, and cheetahs coexist. Ecol. Evol. 6, 8534–8545 (2016).
PubMed  PubMed Central  Article  Google Scholar 

26.
Dröge, E., Creel, S., Becker, M. S. & M’soka, J. Spatial and temporal avoidance of risk within a large carnivore guild. Ecol. Evol.7, 189–199 (2016).

27.
Broekhuis, F., Cozzi, G., Valeix, M., Mcnutt, J. W. & Macdonald, D. W. Risk avoidance in sympatric large carnivores: Reactive or predictive?. J. Anim. Ecol. 82, 1098–1105 (2013).
PubMed  Article  Google Scholar 

28.
Mills, M. G. L. & Mills, M. E. J. Cheetah cub survival revisited: a re-evaluation of the role of predation, especially by lions, and implications for conservation. J. Zool. 292, 136–141 (2014).
Article  Google Scholar 

29.
Broekhuis, F. Natural and anthropogenic drivers of cub recruitment in a large carnivore. Ecol. Evol. 8, 6748–6755 (2018).
PubMed  PubMed Central  Article  Google Scholar 

30.
Orrock, J. L., Preisser, E. L., Grabowski, J. H. & Trussell, G. C. The cost of safety: refuges increase the impact of predation risk in aquatic systems. Ecology 94, 573–579 (2013).
PubMed  Article  Google Scholar 

31.
Donelan, S. C., Grabowski, J. H. & Trussell, G. C. Refuge quality impacts the strength of nonconsumptive effects on prey. Ecology 98, 403–411 (2016).
Article  Google Scholar 

32.
Miller, J. R. B., Ament, J. M. & Schmitz, O. J. Fear on the move: predator hunting mode predicts variation in prey mortality and plasticity in prey spatial response. J. Anim. Ecol. 83, 214–222 (2014).
PubMed  Article  Google Scholar 

33.
Michel, M. J. & Adams, M. M. Differential effects of structural complexity on predator foraging behavior. Behav. Ecol. 20, 313–317 (2009).
Article  Google Scholar 

34.
Blake, L. W. & Gese, E. M. Resource selection by cougars: influence of behavioral state and season. J. Wildl. Manag. 80, 1205–1217 (2016).
Article  Google Scholar 

35.
Hopcraft, J. G. C., Sinclair, A. R. E. & Packer, C. Planning for success: Serengeti lions seek prey accessibility rather than abundance. J. Anim. Ecol. 74, 559–566 (2005).
Article  Google Scholar 

36.
Davies, A. B., Tambling, C. J., Kerley, G. I. H. & Asner, G. P. Effects of vegetation structure on the location of lion kill sites in African thicket. PLoS ONE 11, 1–20 (2016).
Google Scholar 

37.
Camp, M. J., Rachlow, J. L., Woods, B. A., Johnson, T. R. & Shipley, L. A. When to run and when to hide: the influence of concealment, visibility, and proximity to refugia on perceptions of risk. Ethology 118, 1010–1017 (2012).
Article  Google Scholar 

38.
Wilson, A. M. et al. Locomotion dynamics of hunting in wild cheetahs. Nature 498, 185–189 (2013).
ADS  CAS  PubMed  Article  Google Scholar 

39.
Rostro-García, S., Kamler, J. F. & Hunter, L. T. B. To kill, stay or flee: the effects of lions and landscape factors on habitat and kill site selection of cheetahs in South Africa. PLoS ONE 10, e0117743 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

40.
Mills, M. G. L., Broomhall, L. S., du Toit, J. T. & Toit, J. T. Cheetah Acinonyx jubatus feeding ecology in the Kruger National Park and a comparison across African savanna habitats: is the cheetah only a successful hunter on open grassland plains?. Wildl. Biol. 10, 177–186 (2004).
Article  Google Scholar 

41.
Fuller, T. & Sievert, P. Carnivore demography and the consequences of changes in prey availability. in Carnivore conservation (eds. Gittleman, J. L., Funk, S. M., Macdonald, D. & Wayne, R. K.) 163–179 (Cambridge University Press, Cambridge 2001).

42.
Balme, G., Hunter, L. T. B. & Slotow, R. Feeding habitat selection by hunting leopards Panthera pardus in a woodland savanna: prey catchability versus abundance. Anim. Behav. 74, 589–598 (2007).
Article  Google Scholar 

43.
Gosselink, T. E., Van Deelen, T. R., Warner, R. E. & Joselyn, M. G. Temporal habitat partitioning and spatial use of coyotes and red foxes in East-Central Illinois. J. Wildl. Manag. 67, 90–103 (2003).
Article  Google Scholar 

44.
Palomares, F. & Caro, T. M. Interspecific killing among mammalian carnivores. Am. Nat. 153, 492–508 (1999).
CAS  PubMed  Article  Google Scholar 

45.
Brown, J. S. Vigilance, patch use and habitat selection: foraging under predation risk. Evol. Ecol. Res. 1, 49–71 (1999).
Google Scholar 

46.
Creel, S. The control of risk hypothesis: reactive vs. proactive antipredator responses and stress-mediated vs. food-mediated costs of response. Ecol. Lett.21, 947–956 (2018).

47.
Dudeck, B. P., Clinchy, M., Allen, M. C. & Zanette, L. Y. Fear affects parental care, which predicts juvenile survival and exacerbates the total cost of fear on demography. Ecology 99, 127–135 (2018).
PubMed  Article  Google Scholar 

48.
Gigliotti, L. C. et al. Context-dependency of top-down, bottom-up, and density-dependent influences on cheetah demography. J. Anim. Ecol. 2, 449–459 (2020).
Article  Google Scholar 

49.
Relyea, R. A. Morphological and behavioral plasticity of larval anurans in response to different predators. Ecology 82, 523–540 (2001).
Article  Google Scholar 

50.
Suraci, J. P., Clinchy, M., Dill, L. M., Roberts, D. & Zanette, L. Y. Fear of large carnivores causes a trophic cascade. Nat. Commun. 7, 1–7 (2016).
Article  CAS  Google Scholar 

51.
Valeix, M. et al. Behavioral adjustments of African herbivores to predation risk by lions: Spatiotemporal variations influence habitat use. Ecology 90, 23–30 (2009).
CAS  PubMed  Article  Google Scholar 

52.
Sheriff, M. J., Krebs, C. J. & Boonstra, R. The sensitive hare: sublethal effects of predator stress on reproduction in snowshoe hares. J. Anim. Ecol. 78, 1249–1258 (2009).
PubMed  Article  Google Scholar 

53.
Clinchy, M. et al. Multiple measures elucidate glucocorticoid responses to environmental variation in predation threat. Oecologia 166, 607–614 (2011).
ADS  PubMed  Article  Google Scholar 

54.
Travers, M., Clinchy, M., Zanette, L., Boonstra, R. & Williams, T. D. Indirect predator effects on clutch size and the cost of egg production. Ecol. Lett. 13, 980–988 (2010).
PubMed  Google Scholar 

55.
LaManna, J. A. & Martin, T. E. Costs of fear: behavioural and life-history responses to risk and their demographic consequences vary across species. Ecol. Lett. 19, 403–413 (2016).
PubMed  Article  Google Scholar 

56.
Roques, K. G., O’Connor, T. G. & Watkinson, A. R. Dynamics of shrub encroachment in an African savanna: Relative influences of fire, herbivory, rainfall and density dependence. J. Appl. Ecol. 38, 268–280 (2001).
Article  Google Scholar 

57.
Tews, J. & Jeltsch, F. Modelling the impact of climate change on woody plant population dynamics in South African savanna. BMC Ecol. 4, 1–12 (2004).
Article  Google Scholar 

58.
Joubert, D. F., Smit, G. N. & Hoffman, M. T. The role of fire in preventing transitions from a grass dominated state to a bush thickened state in arid savannas. J. Arid Environ. 87, 1–7 (2012).
ADS  Article  Google Scholar 

59.
Lohmann, D., Tietjen, B., Blaum, N., Joubert, D. F. & Jeltsch, F. Prescribed fire as a tool for managing shrub encroachment in semi-arid savanna rangelands. J. Arid Environ. 107, 49–56 (2014).
ADS  Article  Google Scholar 

60.
Durant, S. M. et al. The global decline of cheetah Acinonyx jubatus and what it means for conservation. Proc. Natl. Acad. Sci. 114, 528–533 (2017).
CAS  PubMed  Article  Google Scholar 

61.
Weise, F. J. et al. The distribution and numbers of cheetah (Acinonyx jubatus) in southern Africa. PeerJ 5, e4096 (2017).
PubMed  PubMed Central  Article  Google Scholar 

62.
Prugh, L. R. et al. Designing studies of predation risk for improved inference in carnivore-ungulate systems. Biol. Conserv. 232, 194–207 (2019).
Article  Google Scholar 

63.
Moll, R. J. et al. The many faces of fear: a synthesis of the methodological variation in characterizing predation risk. J. Anim. Ecol. 86, 749–765 (2017).
PubMed  Article  Google Scholar 

64.
Janse van Rensburg, J., McMillan, M., Giżejewska, A. & Fattebert, J. Rainfall predicts seasonal home range size variation in nyala. Afr. J. Ecol.56, 418–423 (2018).

65.
Hunter, L. T. B. The behavioural ecology of reintroduced lions and cheetahs in the Phinda Resource Reserve, Kwazulz-Natal, South Africa. Phd thesis 1–206 (1998).

66.
Caro, T. M. Cheetahs of the Serengeti Plains. (The University of Chicago Press, Chicago, 1994).

67.
Gitzen, R. A., Millspaugh, J. J. & Kernohan, B. J. Bandwidth selection for fixed-kernel analysis of animal utilization distributions. J. Wildl. Manag. 70, 1334–1344 (2006).
Article  Google Scholar 

68.
Thaker, M. et al. Minimizing predation risk in a landscape of multiple predators: effects on the spatial distribution of African ungulates. Ecology 92, 398–407 (2011).
PubMed  Article  Google Scholar 

69.
Moll, R. J., Killion, A. K., Montgomery, R. A., Tambling, C. J. & Hayward, M. W. Spatial patterns of African ungulate aggregation reveal complex but limited risk effects from reintroduced carnivores. Ecology 97, 1123–1134 (2016).
PubMed  Google Scholar 

70.
Kauffman, M. J. et al. Landscape heterogeneity shapes predation in a newly restored predator-prey system. Ecol. Lett. 10, 690–700 (2007).
PubMed  Article  Google Scholar 

71.
Tsalyuk, M., Kelly, M. & Getz, W. M. Improving the prediction of African savanna vegetation variables using time series of MODIS products. ISPRS J. Photogramm. Remote Sens. 131, 77–91 (2017).
ADS  PubMed  PubMed Central  Article  Google Scholar 

72.
Winterstein, S. R. Chi-square tests for intrabrood independence when using the Mayfield method. J. Wildl. Manage. 56, 398–402 (1992).
Article  Google Scholar 

73.
Barker, R. J., White, G. C. & McDougall, M. Movement of Paradise Shelduck between molt sites: a joint multistate-dead recovery mark-recapture model. J. Wildl. Manage. 69, 1194–1201 (2005).
Article  Google Scholar 

74.
Laake, J. L. RMark: An R Interface for Analysis of Capture-Recapture Data with MARK. Version 2.2.4. (2013).

75.
Kelly, M. J. et al. Demography of the Serengeti cheetah (Acinonyx jubatus) population: the first 25 years. J. Zool. 224, 473–488 (1998).
Article  Google Scholar 

76.
Burnham, K. P. & Anderson, D. R. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach (Springer, Berlin, 2002).
Google Scholar 

77.
Arnold, T. W. Uninformative parameters and model selection using Akaike’s Information Criterion. J. Wildl. Manag. 74, 1175–1178 (2010).
Article  Google Scholar  More

DeepHL system architecture
The DeepHL system consists of three server computers. The first one is a web server that receives a trajectory data file from a user and provides analysis results to the user (Intel Xeon E5-2620 v4, 16 cores, 32 GB RAM, Ubuntu 14.04). The second one is a storage server that stores data files and analysis results. The third one is a GPU server that analyzes data provided by the user (Intel Xeon E5-2620 v4, 32 cores, 512 GB RAM, four NVIDIA Quadro P6000, Ubuntu 14.04). Supplementary Information, Algorithm, provides a complete description of the DeepHL method. DeepHL is accessible on the Internet through http://www-mmde.ist.osaka-u.ac.jp/maekawa/deephl/. Supplementary Information, User guide to DeepHL, provides a user guide to DeepHL. In addition, Supplementary Information, Usage of Python-based Software, and Supplementary Software 1 present the Python code of DeepHL.
Preprocessing
An input trajectory is a series of timestamps and X/Y coordinates associated with a class label. To perform position- and rotation-independent analysis, we convert the series into time series of speed and relative angular speed and then standardize them (Supplementary Information, Algorithm). Note that the absolute coordinates of wild animals, which can relate to the distance from a nest or feeding location, for example, are important in understanding behavior of the animals. Hence, DeepHL allows the original coordinates to be input to DeepHL-Net along with the speed and relative angular speed. In addition, other biological time-series sensor data measured by the user can be fed into DeepHL-Net when these time-series data are included in a data file uploaded by the user. For example, a time series of the heading direction of animals obtained from digital compasses can be useful for behavior understanding. Moreover, primitive features usually used in trajectory analysis can be easily fed into DeepHL-Net. DeepHL automatically computes the travel distance from the initial position, the straight-line distance from the initial position, and the angle from the initial position (Supplementary Table 1) as primitive features. Using the web interface of DeepHL, the user can easily select primitive features and other sensor data to be fed into DeepHL-Net (Supplementary Information, User guide to DeepHL). See Supplementary Information, Effect of input features, for effects of input features on classification accuracy. Normally, the inputs of DeepHL-Net are two-dimensional time series, that is, speed and relative angular speed. When we input an additional time series (such as the original coordinates) into DeepHL-Net, the additional time series are added as additional dimensions of the inputs.
Multi-scale layer-wise attention model (DeepHL-Net)
Here, we explain DeepHL-Net shown in Fig. 2f in detail. The input of the model is a time series of primitive features, that is, an lMAX × Nf matrix, where lMAX is the maximum length of the input trajectories and Nf is the dimensionality of the time series, that is, the number of the primitive features. Because the lengths of observed trajectories are not identical to each other in many cases, we fill in missing elements in the matrix with  −1.0 and mask them when we train DeepHL-Net. In each 1D convolutional layer of the convolutional stacks, we extract features by convolving input features through the time dimension using a filter with a width (kernel size) of Ft. We use different filter widths in the four convolutional stacks (3%, 6%, 9%, and 12% of lMAX) to extract features at different levels of scale. We use a stride (step size) of one sample in terms of the time axis. We also use padding to allow the outputs of a layer to have the same length as the layer inputs. In addition, to reduce an overfitting, we employ a dropout, which is a simple regularization technique in which randomly selected neurons are dropped during training44. The dropout rate used in this study is 0.5.
In each LSTM layer of the LSTM stacks, we extract features considering the long-term dependencies of the input features. LSTM is a recurrent neural network architecture with memory cells, and it permits us to learn temporal relationships over a long time scale. LSTM learns long-term dependencies by employing memory cells that hold past information, updating the cell state using write, read, and reset operations with input, output, and forget gates (see Supplementary Information, Algorithm). In addition, we employ dropout to reduce overfitting. The attention information of each layer is computed by using Eq. (1), and then it is multiplied by the layer output. Here, the softmax and tanh functions in Eq. (1) are defined as follows:

$$,{text{softmax}},({x}_{j})=frac{exp ({x}_{j})}{{sum }_{i}exp ({x}_{i})},$$
(2)

$$tanh ({x}_{j})=frac{exp ({x}_{j})-exp (-{x}_{j})}{exp ({x}_{j})+exp (-{x}_{j})}.$$
(3)

Note that parameters in Eq. (1) for each layer, that is, Wa and ba, as well as parameters in the convolutional and LSTM layers are estimated during the network training phase. Here, we introduced the tanh activation function into Eq. (1) to smooth out the output attention values. When an outlying large value is included in WaZT + ba at time t, attention values other than time t become extremely small without using the tanh function. When we visualize a trajectory using such attention values, only a single data point is colored in red, making it difficult for a user to identify important segments.
Training and testing of DeepHL-Net
The DeepHL user can select the parameters of DeepHL-Net used in the analysis, that is, the number of convolutional/LSTM layers and the number of neurons in each layer (default: four layers with 16 neurons). Then, DeepHL-Net is trained on 80% of randomly selected trajectories to minimize the binary classification error of the training data, employing backpropagation based on Adam45 (Supplementary Information, Algorithm). (Note that each trajectory has a class label for binary classification.) Then, the trained DeepHL-Net is tested using the remaining 20% of trajectories to compute the classification accuracy, providing an indication of the degree of difference between the two classes.
Computing the score of each layer
To screen the layers in DeepHL-Net, we compute a score for each layer according to Eq. (4)

$$s({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})={s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})+{s}_{mathrm{it}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}}).$$
(4)

Here, ({A}_{i,{C}_{mathrm{A}}}) is a set of attention vectors calculated from trajectories belonging to class A using the ith layer. In addition, ({A}_{i,{C}_{mathrm{B}}}) is a set of attention vectors calculated from trajectories belonging to class B using the ith layer. As mentioned in the main text, an attention vector from a discriminator layer should have large values within limited segments. Therefore, ({s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})) in Eq. (4) calculates the averaged variance of the attention values normalized by the average length of the trajectories, as described in Eq. (5). When the layer focuses on a part of a trajectory, the variance increases

$${s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})=sqrt{frac{1}{| {A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}}| cdot l({A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}})}sum _{{bf{a}}in {A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}}}V({bf{a}})}.$$
(5)

Note that V(⋅) calculates the variance and l(⋅) calculates the average length of the trajectories. We take the square root of the average variance to derive the average standard deviation. Using (l({A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}})), which calculates the average length of ({A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}}), we normalize the computed variance. Because the softmax function in Eq. (1) ensures that all values sum to 1, resulting in a larger variance for longer trajectories, we normalize the average variance using the average length.
In addition, as mentioned in the main text, the distribution of attention values by the layer for one class should be different from that for another class. Therefore, ({s}_{mathrm{it}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})) calculates the difference between the distributions of the attention values of classes A and B as follows:

$${s}_{mathrm{it}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})=(1-,{mathrm{Intersect}},(h(A_{{i,{C}}_{mathrm{A}}}),h({{A}}_{{i,{C}}_{mathrm{B}}}))).$$
(6)

Here, h(⋅) calculates a normalized histogram of attention with 200 bins, and Intersect(⋅ , ⋅) calculates the area overlap between two histograms, and is described as follows:

$${mathrm{Intersect}},(H_{1},H_{2})=mathop{sum}limits_{i}min (H_{1}(i),H_{2}(i)),$$
(7)

where H1(i) shows the normalized frequency of the ith bin of histogram H1. As described in Eq. (4), the final score is calculated as the sum of the two scores of ({s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})) and ({s}_{mathrm{it}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})).
Here, ({s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})) in Eq. (4) is used to find a layer that focuses only on a portion of a trajectory. Owing to the term, only a small important portion of trajectories is highlighted in many cases, as shown in Figs. 3, 5, and 6, especially for the trajectories of beetles. However, substantial portions of several trajectories of the normal mice are highlighted, as shown in Fig. 4d. Because the characteristics of the normal mouse trajectories are the distance from the initial position, the segments in the trajectories far from the initial position are highlighted.
Computing the correlation between attention values and handcrafted features
To help the user understand the meaning of the highlights, DeepHL automatically computes the Pearson correlation coefficients between the attention values of each layer and handcrafted features computed by DeepHL, as shown in Supplementary Table 1. In addition, the correlation coefficients with sensor data and handcrafted features included in a trajectory data file are automatically computed. Computing the correlation with environmental sensor data can reveal the relationship between a behavior and environmental conditions. If a specific behavior is exhibited only when the temperature is high, for example, we can infer that the behavior relates to the high temperature condition. Furthermore, DeepHL automatically computes the moving average, moving variance, and derivative of each of the above features/sensor data, and then computes the correlation coefficients with the attention values, which are presented to the user (Supplementary Fig. 1).
Computing the difference between distributions of each handcrafted feature for the two classes within highlighted segments
To help the user understand the meaning of the highlights, DeepHL automatically computes the difference between distributions of each handcrafted feature for two classes within highlighted segments. The difference is computed as follows:

$${mathrm{diff}}({A}_{i,{C}_{mathrm{A}}},{F}_{j,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}},{F}_{j,{C}_{mathrm{B}}})=1-,{mathrm{Intersect}},(h(m({{A}}_{{i,{C}}_{mathrm{A}}},{{F}}_{{j,{C}}_{mathrm{A}}})),h(m({{A}}_{{i,{C}}_{mathrm{B}}},{{F}}_{{j,{C}}_{mathrm{B}}}))).$$
(8)

Here, ({F}_{j,{C}_{mathrm{A}}}) is a set of time series of the jth handcrafted feature calculated from trajectories belonging to class A. In addition, m(⋅ , ⋅) is a masking function that extracts feature values within highlighted segments. Because the softmax function in each attention layer ensures that all attention values in a sum of 1, we consider an attention value larger than c/(# time slices) as a potential attended value (c = 1.2 in our implementation).
Data acquisition of worms
Data acquisition was performed according to Yamazoe-Umemoto et al.22. In brief, several worms were placed in the center of an agar plate in a 9-cm Petri dish, 30% 2-nonanone (v/v, EtOH) was spotted on the left side of the plate, which was covered by a lid and placed on the bench upside down. Then, the images of the plate were captured with a high-resolution USB camera for 12 min at 1 Hz. Because the worms do not exhibit odor avoidance behavior during the first 2 min because of the rapid increase in odor concentration46, the data for the following 10 min (i.e., 600 s) was used. From the images, individual worms were identified and the position of the centroid was recorded by an image processing software Move-tr/2D (v. 8.31; Library Inc., Japan). The number of recorded trajectories is 325 (Supplementary Table 2). The comparison was between the naive worms (control class) and the worms after preexposure to the odor (preexposed class).
DeepHL analysis of worms
A multivariate time series of movement speed, relative angular speed, distances from the initial position, and angle from the initial position extracted from the time series of trajectories was fed into DeepHL-Net, yielding a binary classification accuracy of 93.9%, where 20% of the data are used as test data. The discriminator layer used in this investigation has the highest score of all layers. As shown in Fig. 3d, which was calculated from the moving variance of the speed within highlighted segments, we can state that the changes in the speed of preexposed worms is larger than those of control worms. Figure 3e shows spectrograms of the speed calculated from entire trajectories (Fig. 3c) with a 128-s wide sliding window shifted in 1-sample intervals. In addition, Fig. 3f shows histograms of the dominant frequency of speed calculated from entire trajectories using the 128-s wide sliding window shifted in 1-sample intervals. These results also indicate the difference in the frequency of speed between the preexposed and control worms. Our investigation revealed that the dominant frequency of speed significantly differs between the preexposed and control worms using GLMM with Gaussian distributions (t = −6.60; d.f. = 322.8; p = 1.68 × 10−10, effect size(r2) = 0.232). The p value is two sided. Individual factors were treated as random effects. The number of data points for the control class is n = 76, 784 and that for the preexposed class is n = 75, 750. We used GLMM with Gaussian distributions because the objective variable has a continuous value and we used the lmerTest package (v. 2.0–36) of R (v. 3.4.3) for the analysis.
Data acquisition of mice
We collected 52 trajectories of normal mice and unilateral 6-hydroxydopamine (OHDA) lesion mouse models of PD while they freely moved for 10 min in an open field (60 × 55 cm2, wall height = 20 cm; normal: 22, PD: 30). The trajectories were detected by the animal’s head position, which was captured by an overhead digital video camera (60 fps). Two sets of small red and green light-emitting diodes were mounted above the animal’s head so that it could be located in each frame. Custom softwares based on Matlab (R2018b, Mathworks, MA, USA) and LabVIEW (Labview 2018, National Instruments, TX, USA) were used for tracking. We then created 30-s segments by splitting each trajectory because training a DNN requires a number of trajectories. We used 966 segments in total (normal: 374, PD: 592) collected from nine C57BL/6J mice (normal: 5, PD: 4). Note that we excluded 30-s segments that contain no movements of a mouse.
DeepHL analysis of mice
Movement speed, relative angular speed, travel distances, straight-line and travel distances from the initial position, and angle from the initial position were fed into our model. The accuracy for the binary classification of normal and 6-OHDA model mice was 74.7%, where 20% of the data are used as test data. The score of the discriminator layer was the highest of all LSTM layers and the sixth highest of all layers. Our investigation revealed that the behavior of visiting locations far away from the initial position can be characteristic of normal mice.
To evaluate PD symptoms from animal behaviors, previous studies have exclusively focused on the movement speed of animals in the open-field tests (frequency and bout duration of ambulation as well as immobility or fine movement) because typical symptoms in the animal model of PD are thought to be slowness of movement and a paucity of spontaneous movements. As shown in Fig. 4e–g, we found significant differences in average movement speed during ambulation periods, average movement speed during fine movement periods, and average maximum distance within a ±60-s window in a session. These differences were derived from the findings of DeepHL using the two-sided Wilcoxon rank-sum test (W = 544, p = 3.486 × 10−5, effect size (Cliff’s delta) = −0.648; W = 511, p = 5.869 × 10−4, effect size (Cliff’s delta) = −0.548; W = 521, p = 2.666 × 10−4, effect size (Cliff’s delta) = −0.579). The 95% confidence intervals are [1.222, 3.481], [0.139, 0.468], and [13.726, 43.175], respectively. We used the exactRankTests package (v. 0.8–29) of R (v. 3.2.3). Note that these behavioral features are extracted from original 10-min trajectories.
The maximum distance, which was derived from a finding of DeepHL, is more useful for evaluating the PD symptoms than conventional measures based on the movement speed. Note that the new feature is designed based on an insight drawn from an analysis by deep learning. These results suggest that DeepHL helps find a novel measure not directly linked to the movement speed, that is, a straight-line distance within a certain time window. When the aim of an animal is to visit all locations in an area, the travel distance over a short duration commonly becomes longer. Besides, it is well known that rodents, including mice and rats, spontaneously prefer to explore an environment, particularly in novel places. Thus, DeepHL may capture the fact that the abnormal behavior of the 6-OHDA lesion model of PD hinders such spontaneous behavioral traits of normal mice. Indeed, the 6-OHDA lesion mouse model appears to remain in the same place. Although this hypothesis should be verified based on the causality between behavioral traits and neural activity patterns underlying PD symptoms using neuronal recording together with its optogenetic manipulation in the basal ganglia and motor cortex23, it is beyond the scope of this study.
Behavioral features of mice
According to Kravitz et al.23, ambulation was defined as periods when the velocity of the animal’s center point averaged >2 cm/s for at least 0.5 s. Immobility was defined as continuous periods of time during which the average change of the trajectory was More