Philippa Kaur

1.
Fox, L. R. Cannibalism in natural populations. Annu. Rev. Ecol. Syst. 6, 87–106 (1975).
Article  Google Scholar 
2.
Polis, G. A. The evolution and dynamics of intraspecific predation. Annu. Rev. Ecol. Evol. Syst. 12, 225–251 (1981).
Article  Google Scholar 

3.
Elgar, M. A. & Crespi, B. J. Ecology and evolution of cannibalism. In Cannibalism: ecology and evolution among diverse taxa (eds Elgar, M. A. & Crespi, B. J.) 1–12 (Oxford University Press, Oxford, 1992).
Google Scholar 

4.
Richardson, M. L., Mitchell, R. F., Reagel, P. F. & Hanks, L. M. Causes and consequences of cannibalism in noncarnivorous insects. Annu. Rev. Entomol. 55, 39–53 (2010).
CAS  PubMed  Article  Google Scholar 

5.
Vilaça, A. Relations between funerary cannibalism and warfare cannibalism: The question of predation. Ethnos 65, 83–106 (2000).
Article  Google Scholar 

6.
Lopez-Riquelme, G. O. & Fanjul-Moles, M. L. The funeral ways of social insects. Social strategies for corpse disposal. Trends Entomol. 9, 71–129 (2013).
Google Scholar 

7.
Walls, S. C. & Roudebush, R. E. Reduced aggression toward siblings as evidence of kin recognition in cannibalistic salamanders. Am. Nat 138, 1027–1038 (1991).
Article  Google Scholar 

8.
Pfennig, D. W. Cannibalistic tadpoles that pose the greatest threat to kin are most likely to discriminate kin. Proc. R. Soc. Lond. B 266, 57–61 (1999).
Article  Google Scholar 

9.
Bilde, T. & Lubin, Y. Kin recognition and cannibalism in a subsocial spider. J. Evolut. Biol. 14, 959–966 (2001).
Article  Google Scholar 

10.
Santana, A. F. K., Roselino, A. C., Cappelari, F. A. & Zucoloto, F. S. Cannibalism in insects. In Insect bioecology and nutrition for integrated pest management (eds Panizzi, A. R. & Parra, J. R. P.) 177–190 (CRC Press, Boca Raton, 2012).
Google Scholar 

11.
Hölldobler, B. & Wilson, E. O. The ants (The Belknap Press of Harvard University, London, 1990).
Google Scholar 

12.
Schmickl, T. & Crailsheim, K. Cannibalism and early capping: strategy of honeybee colonies in times of experimental pollen shortage. J. Comput. Physiol. A 187, 541–547 (2001).
CAS  Article  Google Scholar 

13.
Sun, Q. & Zhou, X. Corpse management in social insects. Int. J. Biol. Sci. 9, 313–321 (2013).
PubMed  PubMed Central  Article  Google Scholar 

14.
Davis, H. E., Meconcelli, S., Rudek, R. & McMahon, D. P. Termites shape their collective behavioural response based on stage of infection. Sci. Rep. 8, 14433. https://doi.org/10.1038/s41598-018-32721-7 (2018).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

15.
Mabelis, A. A. Wood ant wars: the relationship between aggression and predation in the red wood ant (Formica polyctena Först.). Neth. J. Zool. 29, 451–620 (1979).
Article  Google Scholar 

16.
Driessen, G. J. J., Van Raalte, ATh. & De Bruyn, G. Cannibalism in the red wood ant, Formica polyctena (Hymenoptera: Formicidae). Oecologia 63, 13–22 (1984).
ADS  PubMed  Article  Google Scholar 

17.
Yao, M. et al. The ancient chemistry of avoiding risks of predation and disease. Evol. Biol. 36, 267–281 (2009).
Article  Google Scholar 

18.
Visscher, P. K. The honey bee way of death: Necrophoric behaviour in Apis mellifera colonies. Anim. Behav. 31, 1070–1076 (1983).
Article  Google Scholar 

19.
Oi, D. H. & Pereira, R. M. Ant behavior and microbial pathogens (Hymenoptera: Formicidae). Florida Entomol. 76, 63–74 (1993).
Article  Google Scholar 

20.
Nazzi, F., Della Vedova, G. & D’Agaro, M. A semiochemical from brood cells infested by Varroa destructor triggers hygienic behaviour in Apis mellifera. Apidologie 35, 65–70 (2004).
CAS  Article  Google Scholar 

21.
Renucci, M., Tirrard, A. & Provost, E. Complex undertaking behavior in Temnothorax lichtensteini ant colonies: From corpse-burying behavior to necrophoric behavior. Insect. Soc. 58, 9–16 (2011).
Article  Google Scholar 

22.
Diez, L., Le Borgne, H., Lejeune, P. & Detrain, C. Who brings out the dead? Necrophoresis in the red ant Myrmica rubra. Anim. Behav. 6, 1259–1264 (2013).
Article  Google Scholar 

23.
Baracchi, D., Fadda, A. & Turillazzi, S. Evidence for antiseptic behaviour towards sick adult bees in honey bee colonies. J. Insect. Physiol. 58, 1589–1596 (2012).
CAS  PubMed  Article  Google Scholar 

24.
Pull, Ch. D. et al. Destructive disinfection of infected brood prevents systemic disease spread in ant colonies. eLife 7, e32073. https://doi.org/10.7554/eLife.32073 (2018).
Article  PubMed  PubMed Central  Google Scholar 

25.
Leclerc, J.-B. & Detrain, C. Ants detect but do not discriminate diseased workers within their nest. Sci. Nat. 103, 70. https://doi.org/10.1007/s00114-016-1394-8 (2016).
CAS  Article  Google Scholar 

26.
Williams, T. & Hernandez, O. Costs of cannibalism in the presence of an iridovirus pathogen of Spodoptera frugiperda. Ecol. Entomol. 31, 106–113 (2006).
Article  Google Scholar 

27.
Rudolf, V. H. W. & Antonovics, J. Disease transmission by cannibalism: rare event or common occurrence?. Proc. R. Soc. Lond. B 274, 1205–1210 (2007).
Google Scholar 

28.
Sadeh, A. & Rosenheim, J. A. Cannibalism amplifies the spread of vertically transmitted pathogens. Ecology 97, 1994–2002 (2016).
PubMed  Article  Google Scholar 

29.
Claessen, D., de Roos, A. M. & Persson, L. Population dynamic theory of size-dependent cannibalism. Proc. R. Soc. Lond. B 271, 333–340 (2004).
Article  Google Scholar 

30.
Pfennig, D. W., Ho, S. G. & Hoffman, E. A. Pathogen transmission as a selective force against cannibalism. Anim. Behav. 55, 1255–1261 (1998).
CAS  PubMed  Article  Google Scholar 

31.
Loreto, R. G. & Hughes, D. P. Disease in the society: infectious cadavers result in collapse of ant sub-colonies. PLoS ONE 11, e0160820 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

32.
Hughes, W. H., Eilenberg, J. & Boomsmal, J. J. Trade-offs in group living: Transmission and disease resistance in leaf-cutting ants. Proc. R. Soc. Lond. B 269, 1811–1819 (2002).
Article  Google Scholar 

33.
Cremer, S. & Sixt, M. Analogies in the evolution of individual and social immunity. Proc. R. Soc. Lond. B 364, 129–142 (2009).
Google Scholar 

34.
Konrad, M. et al. Social transfer of pathogenic fungus promotes active immunisation in ant colonies. PLoS ONE 10, 1–15 (2012).
Google Scholar 

35.
Liu, L., Ganghua, L., Pengdong, S., Chaoliang, L. & Quiying, H. Experimental verification and molecular basis of active immunization against fungal pathogens in termites. Sci. Rep. 5, 15106. https://doi.org/10.1038/srep15106 (2015).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

36.
Marikovsky, P. I. On some features of behaviour of the ants Formica rufa L. infected with fungus disease. Insect. Soc. 2, 173–179 (1962).
Article  Google Scholar 

37.
Rutkowski, T. et al. Ants trapped for years in an old bunker; survival by cannibalism and eventual escape. J. Hymenopt. Res. 72, 177–184 (2019).
Article  Google Scholar 

38.
Seifert, B. Die Ameisen Mittel- und Nordeuropas (Lutra-Verlags-und Vertriebsgesellschaft, Görlitz, 2007).
Google Scholar 

39.
Czechowski, W., Radchenko, A., Czechowska, W. & Vepsäläinen, K. The ants of Poland with reference to the myrmecofauna of Europe. Fauna Poloniae (n.s.) 4. (Natura Optima Dux Foundation, 2012).

40.
Meyling, N. V. & Eilenberg, J. Ecology of the entomopathogenic fungi Beauveria bassiana and Metarhizium anisopliae in temperate agroecosystems: Potential for conservation biological control. Biol. Control 43, 145–155 (2007).
Article  Google Scholar 

41.
Reber, A. & Chapuisat, M. Diversity, prevalence and virulence of fungal entomopathogens in colonies of the ant Formica selysi. Insect. Soc. 59, 231–239 (2012).
Article  Google Scholar 

42.
Hajek, A. E. & St. Leger, R. J. Interactions between fungal pathogens and insect hosts. Annu. Rev. Entomol. 39, 293–322 (1994).
Article  Google Scholar 

43.
Maák, I. et al. Cues or meaningless objects? Differential responses of the ant Formica cinerea to corpses of competitors and enslavers. Anim. Behav. 91, 53–59 (2014).
Article  Google Scholar 

44.
Csata, E. & Dussutour, A. Nutrient regulation in ants (Hymenoptera: Formicidae): A review. Myrmecol. News 29, 111–124 (2019).
Google Scholar 

45.
Nonacs, P. Death in the distance: Mortality risk as information for foraging ants. Behaviour 112, 23–35 (1990).
Article  Google Scholar 

46.
Roces, F. & Núṅez, J. A. Information about food quality influences load-size selection in recruited leaf-cutting ants. Anim. Behav. 45, 135–143 (1993).
Article  Google Scholar 

47.
Song, D., Hu, X. P. & Su, N.-Y. Survivorship, cannibalism, body weight loss, necrophagy, and entombement in laboratory groups of the Formosan subterranean termite, Coptotermes formosanus under starvation (Isoptera: Rhinotermitidae). Sociobiology 47, 27–39 (2006).
Google Scholar 

48.
Heifig, I., Lima, J. T., Janei, V. & Costa-Leonardo, A. M. Effects of group size and starvation on survival of the Asian subterranean termite Coptotermes gestroi (Isoptera: Rhinotermitidae). Austral Entomol. 57, 279–284 (2017).
Article  Google Scholar 

49.
Pompilio, L., Kacelnik, A. & Behmer, S. T. State-dependent learned valuation drives choice in an invertebrate. Science 311, 1613–1615 (2006).
ADS  CAS  PubMed  Article  Google Scholar 

50.
Akino, T. & Yamaoka, R. Origin of oleic acid: Corpse recognition signal in the ant Formica japonica Motschlsky (Hymenoptera: Formicidae). Jpn. J. Appl. Entomol. Z. 40, 265–271 (1996).
CAS  Article  Google Scholar 

51.
Chouvenc, T., Robert, A., Sémon, E. & Bordereau, C. Burial behaviour by dealates of the termite Pseudacanthotermes spiniger (Termitidae, Macrotermitinae) induced by chemical signals from termite corpses. Insect. Soc. 59, 119–125 (2012).
Article  Google Scholar 

52.
Kok-Boon, N., Beng-Keok, Y., Kunio, T., Tsuyoshi, Y. & Chow-Yang, L. Do termites avoid carcasses? Behavioral responses depend on the nature of the carcasses. PLoS ONE 7, 1–11 (2012).
Google Scholar 

53.
Diez, L., Moquet, L. & Detrain, C. Post-mortem changes in chemical profile and their influence on corpse removal in ants. J. Chem. Ecol. 39, 1424–1432 (2013).
CAS  PubMed  Article  Google Scholar 

54.
Bignell, D. E., Roisin, Y. & Lo, N. Biology of Termites: A modern synthesis (Springer, Berlin, 2010).
Google Scholar 

55.
Dlusskij, G. M. Ants of the genus Formica (Hymenoptera, Formicidae, g. Formica) (Nauka, Moscow, 1967) (in Russian).
Google Scholar 

56.
Czechowski, W. Ants cemeteries. Przegląd Zoologiczny 20, 417–427 (1976) (in Polish with English summary).
Google Scholar 

57.
Czechowski, W. Around nest cemeteries of Myrmica schencki Em. (Hymenoptera: Formicidae): their origin and a possible significance. Pol. J. Ecol. 56, 359–363 (2008).
Google Scholar 

58.
Gibb, H. Experimental evidence for mediation of competition by habitat succession. Ecology 92, 1871–1878 (2011).
CAS  PubMed  Article  Google Scholar 

59.
Chouvenc, T., Su, N.-Y. & Elliott, M. L. Interaction between the subterranean termite Reticulitermes flavipes (Isoptera: Rhinotermitidae) and the entomopathogenic fungus Metarhizium anisopliae in foraging arenas. J. Econ. Entomol. 101, 885–893 (2008).
CAS  PubMed  Article  Google Scholar 

60.
Yanagawa, A., Yokohari, F. & Shimizu, S. The role of antennae in removing entomopathogenic fungi from cuticle of the termite Coptotermes formosanus. . J. Insect Sci. 9, 1–9 (2009).
Article  Google Scholar 

61.
Tranter, Ch., LeFevre, L., Evison, S. E. F. & Hughes, W. O. H. Threat detection: Contextual recognition and response to parasites by ants. Behav. Ecol. 26, 396–405 (2015).
Article  Google Scholar 

62.
Bonadies, E., Wcislo, W. T., Gálvez, D., Hughes, W. O. H. & Fernández-Marin, H. Hygiene defense behaviors used by a fungus-growing ant depend on the fungal pathogen stages. Insects 10, 130 (2019).
PubMed Central  Article  PubMed  Google Scholar 

63.
Simone-Finstrom, M. D. & Spivak, M. Increased resin collection after parasite challenge: A case of self-medication in honey bees?. PLoS ONE 7, e34601. https://doi.org/10.1371/journal.pone.0034601 (2012).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

64.
Brütsch, T. & Chapuisat, M. Wood ants protect their brood with tree resin. Anim. Behav. 93, 157–161 (2014).
Article  Google Scholar 

65.
Ormond, E. L., Thomas, A. P. M., Pell, J. K., Freeman, S. N. & Roy, H. E. Avoidance of a generalist entomopathogenic fungus by the ladybird Coccinella septempunctata. FEMS Microbiol. Ecol. 77, 229–237 (2011).
CAS  PubMed  Article  Google Scholar 

66.
Fernández-Marín, H., Zimmerman, J. K., Rehner, S. A. & Wcislo, W. T. Active use of the metapleural glands by ants in controlling fungal infection. Proc. R. Soc. Lond. B 273, 1689–1695 (2006).
Google Scholar 

67.
Tragust, S. et al. Ants disinfect fungus-exposed brood by oral uptake and spread of their poison. Curr. Biol. 23, 1–7 (2013).
Article  CAS  Google Scholar 

68.
Tragust, S., Herrmann, C., Häfner, J., Braasch, R., Tilgen, Ch., Hoock, M., Milidakis, M. A., Gross, R. & Feldhaar, H. Formicine ants swallow their highly acidic poison for gut microbial selection and control. bioRxiv preprint https://doi.org/10.1101/2020.02.13.947432 (2020).

69.
Cremer, S., Pull, Ch. D. & Fürst, M. A. Social immunity: emergence and evolution of colony-level disease protection. Annu. Rev. Entomol. 63, 105–123 (2018).
CAS  PubMed  Article  Google Scholar 

70.
Rosengaus, R. B., Jordan, C., Lefebvre, M. L. & Traniello, J. F. A. Pathogen alarm behavior in a termite: A new form of communication in social insects. Naturwissenschaften 86, 544–548 (1999).
ADS  CAS  PubMed  Article  Google Scholar 

71.
Hernandez-Lopez, J., Reissberger-Gallé, U., Crailsheim, K. & Schuehly, W. Cuticular hydrocarbon cues of immune-challenged workers elicit immune activation in honeybee queens. Mol. Ecol. 26, 3062–3073 (2017).
CAS  PubMed  Article  Google Scholar 

72.
Chouvenc, T. & Su, N.-Y. When subterranean termites challenge the rules of fungal epizootics. PLoS ONE 7, 84. https://doi.org/10.1371/journal.pone.0034484 (2012).
CAS  Article  Google Scholar 

73.
Csata, E., Erős, K. & Markó, B. Effects of the ectoparasitic fungus Rickia wasmannii on its ant host Myrmica scabrinodis: changes in host mortality and behavior. Insect. Soc. 61, 247–252 (2014).
Article  Google Scholar 

74.
Diez, L., Urbain, L., Lejeune, Ph. & Detrain, C. Emergency measures: adaptive response to pathogen intrusion in the ant nest. Behav. Process. 116, 80–86 (2015).
Article  Google Scholar 

75.
Qui, H.-L. et al. Differential necrophoric behaviour of the ant Solenopsis invicta towards fungal infected corpses of workers and pupae. Bull. Entomol. Res. 105, 607–614 (2015).
Article  CAS  Google Scholar 

76.
Pereira, H. & Detrain, C. Pathogen avoidance and prey discrimination in ants. R. Soc. Open Sci. 7, 191705 (2020).
ADS  PubMed  PubMed Central  Article  Google Scholar 

77.
Cremer, S., Armitage, S. A. O. & Schmid-Hempel, P. Social immunity. Curr. Biol. 17, 693–702 (2007).
Article  CAS  Google Scholar 

78.
Pull, Ch. D. & Cremer, S. Co-founding ant queens prevent disease by performing prophylactic undertaking behaviour. BMC Evol. Biol. 219, 17. https://doi.org/10.1186/s12862-017-1062-4 (2017).
Article  Google Scholar 

79.
Kramm, K. R., West, D. F. & Rockenbach, P. G. Pathogens of termites: transfer of the entomopathogen Metarhizium anisopliae between the termites of Reticulitermes sp.. J. Invertebr. Pathol. 40, 1–6 (1982).
Article  Google Scholar 

80.
Kesäniemi, J., Koskimäki, J. J. & Jurvansuu, J. Corpse management of the invasive Argentine ant inhibits growth of pathogenic fungi. Sci. Rep. 9, 7593. https://doi.org/10.1038/s41598-019-44144-z (2019).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

81.
Greenwald, E. E., Baltiansky, L. & Feinerman, O. Individual crop loads provide local control for collective food intake in ant colonies. eLife 7, e31730 (2018).
PubMed  PubMed Central  Article  Google Scholar 

82.
Horstmann, K. Untersuchungen über den Nahrungserwerb der Waldameisen (Formica polyctena Foerster) im Eichenwald. Oecologia 5, 138–157 (1970).
ADS  PubMed  Article  Google Scholar 

83.
Bhatkar, A. & Whitcomb, W. H. Artificial diet for rearing various species of ants. Florida Entomol. 53, 229–232 (1970).
Article  Google Scholar 

84.
Choe, D. H. & Rust, M. K. Horizontal transfer of insecticides in laboratory colonies of the argentine ant (Hymenoptera: Formicidae). J. Econ. Entomol. 101, 1397–1405 (2008).
CAS  PubMed  Article  Google Scholar 

85.
Pereira, R. M. & Stimac, J. L. Transmission of Beauveria bassiana within nests of Solenopsis invicta (Hymenoptera: Formicidae) in the laboratory. Environ. Entomol. 21, 1427–1432 (1992).
Article  Google Scholar 

86.
Liu, H., Skinner, M., Parker, B. L. & Brownbridge, M. Pathogenicity of Beauveria bassiana, Metarhizium anisopliae (Deuteromycotina: Hyphomycetes), and other entomopathogenic fungi against Lygus lineolaris (Hemiptera: Miridae). J. Econ. Entomol. 95, 675–681 (2002).
PubMed  Article  Google Scholar 

87.
Loreto, R. G. & Hughes, D. P. Disease dynamics in ants. Adv. Genet. 94, 287–306. https://doi.org/10.1016/bs.adgen.2015.12.005 (2016).
CAS  Article  PubMed  Google Scholar 

88.
R Core Team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing, 2017). https://www.R-project.org/.

89.
Therneau, T. coxme: Mixed Effects Cox Models. R package version 2.2-5. https://CRAN.R-project.org/package=coxme (2015).

90.
Bates, D., Maechler, M., Bolker, B. & Walker, S. lme4: Linear mixed-effects models using Eigen and S4. R package version 1.0-5. https://CRAN.R-project.org/package=lme4 (2013).

91.
Bartoń, K. MuMIn: Multi-model inference. R package version 1.9.13. https://CRAN.R-project.org/package=MuMIn (2013).

92.
Grueber, C. E., Nakagawa, S., Laws, R. J. & Jamieson, I. G. Multimodel inference in ecology and evolution: Challenges and solutions. J. Evol. Biol. 24, 699–711 (2011).
CAS  PubMed  Article  Google Scholar  More

The nature of the slug movement data collected in our field experiment (i.e. position on or beneath the soil surface observed at discrete moments of time) makes it possible to analyse the slug locomotory track in terms of the discrete-time random movement framework12,16,35,36. Within this framework, a curvilinear movement path is approximated by a broken line (see Fig. 2) and the movement of an individual slug is parameterized by the following frequency distributions:
1.
The distribution of the step sizes along the movement path (i.e. the distance between sequential pairs of recorded positions; Fig. 2) or the corresponding average speed

2.
The distribution of turning angle (the angle between the straight lines drawn between sequential pairs of recorded positions; Fig. 2).

Figure 2

A sketch of animal movement path and its discretization (adapted from36). (a) The original movement path is normally curvilinear. (b) Due to the limitations of the radio-tracking technique, position of the animal is only known at certain discrete moments of time; correspondingly, the curve is approximated by a broken line. (c) The movement path as a broken line is fully described by the sequence of the step sizes (lengths) along the path, i.e. the distances travelled between any two sequential recorded positions, and the sequence of the corresponding turning angles.

Full size image

Once all the information is available, it is possible to calculate the mean squared displacement as a function of time12,37. Additionally, in case the movement consists of alternating periods of active movement and immobility (periods with no recorded displacement resulting from feeding or inactivity, hereafter referred to as“resting time”), one should also consider the distribution of the corresponding periods.
Speed, squared displacements and the straightness index
It is apparent from the data that slug movement is intermittent, with periods of locomotion interspersed between periods in which they remain motionless. Tables 1 and 2 show, for the sparse and dense releases respectively, the number of ‘active’ time intervals when the slugs were moving. Periods during which slugs were motionless are marked by the zeros in Tables 1 and 2, but all these individuals resumed their movement during the following hours, confirming that they were alive throughout the assessment period. We therefore retain the zeros in the data for the subsequent analysis.
Table 1 Slug mean speed (averaged over the whole movement path), the mean SSD (see Eqs. (3) and (5), respectively) and the straightness index in the case of sparse release for each of 17 slugs used in the experiment. Here the straightness index is calculated using Eq. (4) where the values of the step size are immediately available from our field data.
Full size table

Table 2 Slug mean speed (averaged over the whole movement path), the mean SSD (see Eqs. (3) and (5), respectively) and the straightness index in the case of dense release for each of 11 slugs used in the experiment. Here the straightness index is calculated using Eq. (4) where the values of the step size are immediately available from our field data.
Full size table

The baseline discrete-time framework considers animal position at equidistant moments of time. However, in the field experiment (as described in the previous section), time taken to locate slugs at each assessment resulted in the time interval varying between measurements (sparse release treatment: 27–87 mins; dense release treatment: 20–103 mins). The step size, i.e. the displacement during one time interval, depends in part on the duration of that interval, hence risking bias in the results. We address this issue by scaling the step size by the duration of the corresponding time interval, i.e. by considering the average speed during the step:

$$begin{aligned} v_k(i)=, & {} frac{|Delta {mathbf{r}|_k(i)}}{Delta {t}_k(i)}, quad i=1,2,ldots ,N, end{aligned}$$
(1)

where

$$begin{aligned} |Delta {mathbf{r}|_k(i)}=, & {} |mathbf{r}_k(t_i)-mathbf{r}_k(t_{i-1})|, end{aligned}$$
(2)

is the displacement of the kth slug during the ith time interval, i.e. the distance between the two sequential positions in the field. Here N is the total number of steps made by the given slug during the full period of the experiment (in our field data, for all slugs (N=10)).
For each individual slug, we then calculate the mean speed over all steps along the movement path:

$$begin{aligned} _k=, & {} frac{1}{N} sum ^{N}_{i=1} v_k(i). end{aligned}$$
(3)

The results for the sparse and dense releases are shown in Tables 1 and 2, respectively; see also Fig. 3a.
The mean speed of slug movement, although being an important factor for slug dispersal, does not provide enough information about the rate at which the slug increases its linear distance from the point of release, because it does not provide information on the frequency of turning or the turning angle. In order to take that into account, we calculate the straightness index35, i.e. the ratio of the total displacement (distance between the point of release and the final position at the end of the experiment) to the total distance travelled along the path:

$$begin{aligned} s_k= & {} |mathbf{r}_k(t_N)-mathbf{r}_k(t_0)|/left( sum ^{N}_{i=1}|Delta {mathbf{r}}|_k(i) right) , end{aligned}$$
(4)

where (t_0) is the time of slug release and (t_N) is the time of the final observation. The actual distance travelled is approximated by the length of the corresponding broken line (see the dark solid line in Fig. 2).
Figure 3

(a) Slug mean spead and (b) slug mean SSD, black diamonds for the sparse release and red circles for the dense release.

Full size image

The straightness index quantifies the amount of turning (a combination of the frequency and angles of turns) along the whole movement path, i.e. over the whole observation time, but it says nothing about the rate of turning on the shorter time scale of a single ‘step’ along the movement path. To account for this, along with the mean speed we calculate the mean scaled squared displacement (SSD):

$$begin{aligned} langle sigma ^2 rangle _k=, & {} frac{1}{N} sum ^{N}_{i=1} sigma ^2_k(i) qquad text{ where }qquad sigma ^2_k(i)~=~frac{|Delta {mathbf{r}|^2_k(i)}}{Delta {t}_k(i)}, end{aligned}$$
(5)

see Tables 1 and 2 and Fig. 3b. For the same value of mean speed, a larger value of the SSD corresponds to a straighter movement on the timescale of a single step, with a smaller turning rate.
An immediate observation from visual analysis of the data shown in Fig. 3 is that both slug speed and the SSD are smaller in the case of dense release than in the sparse release. Therefore, a preliminary conclusion can be drawn that average slug movement is slower in the dense release compared to the sparse release treatment.
Turning angles
Figure 4

Frequency distribution of the turning angle in the case of (a) sparse and (b) dense releases of slugs. In calculating the turning angle, the periods of no movement were disregarded. The red curve shows the best-fitting of the data with the exponential function; see details in the text.

Full size image

We now proceed to analyse the distribution of turning angles. The histogram of different values of the angle is shown in Fig. 4. Let us consider first the case of sparse release (see Fig. 4a). We readily observe that the distribution is roughly symmetrical and has a clear maximum at (theta _T=0). The latter indicates that, on this timescale, slug movement is better described as the CRW than the standard diffusion1,16. Indeed, the standard diffusion (also known as the simple random walk) assumes that there is no bias in the movement direction, in particular there is no correlation in the movement direction in the intervals before and after the recorded position, which means that the turning angle is uniformly distributed over the whole circle. On the contrary, in the case where a correlation between the movement directions exists (hence resulting in the CRW), the distribution of the turning angle becomes hump-shaped. This is in agreement with the results of previous studies on animal movement (in particular, invertebrates12,38) as well as a general theoretical argument13.
In order to provide a more quantitative insight, we look for a functional description of the turning angle distribution using several distributions that are commonly used in movement ecology. The results are shown in Table 3. We establish that the turning angle data are best described by the exponential distribution. Somewhat unexpectedly, it outperforms the Von Mises distribution, although the latter is often regarded as a benchmark and its use has some theoretical justification1. However, the exponential distribution of the turning angle has previously been observed in movement data on some other species, e.g. on swimming invertebrates39.
Table 3 The (r^2) values for the turning angle movement data (in case of sparsely released slugs) described by different standard frequency distributions. The corresponding data are shown in Fig. 4a.
Full size table

The distribution of turning angle obtained in the case of dense release exhibit different features; see Fig. 4b. However, in this case, the distribution is not symmetric and has a clear bias towards positive values: the mean turning angle corresponding to the data shown in Fig. 4b is (langle theta _T rangle = 0.772approx pi /4). Since the slugs used in the dense release are from the same cohort as those used in the sparse release, we consider this bias as an effect of the slug density: the movement pattern of an individual slug is affected by the presence of con-specifics. We discuss possible specific mechanisms for the responsiveness to this factor in the Discussion.
An attempt to describe the turning angle data from the dense release by a symmetric distribution returns low values of (r^2) (see Suppl. Appendix A.1). However, the accuracy of data fitting comparable with the sparse release can be achieved by using an asymmetric distribution, i.e. where the corresponding function has different parameters for the positive and negative values of the angle. The results are shown in Table 4.
Table 4 The (r^2) values for the turning angle movement data (in case of densely released slugs) described by asymmetric frequency distributions. In calculating the turning angle, the periods of no movement were disregarded; the corresponding data are shown in Fig. 4b.
Full size table

The turning angle data shown in Fig. 4 were obtained using all active steps along the movement paths. However, since periods of slug movement alternate with periods of resting, it may raise the question of the relevance of the turning angle at the locations where slugs remained motionless for some time. In order to check the robustness of our results, we now repeat the analysis to calculate the turning angle differently by omitting the segments adjoined with the rest position. The results are shown in Fig. 5. In this case, a reliable fit may not be possible due to there being insufficient data. However, a visual inspection of the corresponding histograms suggests that the main properties of the turning angle distribution agree with those observed above for the bigger data set. Namely, in both cases the distribution has a clear maximum at (theta _T=0) (this is seen particularly well in the case of sparse release). In the case of sparse release the distribution is approximately symmetric, while in the case of dense release there is a clear bias towards positive values. We therefore conclude that the properties of the turning angle distribution are robust with regard to the details of its definition.
Figure 5

Frequency distribution of the turning angle in case of (a) sparse release, (b) dense release. The turning angle is only calculated for consecutive movements, i.e. if a slug does not move during a time step then its previous angle of movement is not used.

Full size image

Movement and resting times
Figure 6

Distribution of the proportion of the total time spent in movement in case of (a) sparse release and (b) dense release. The red curve shows the best-fit of the data with the normal distribution; see details in the text.

Full size image

Our field data shows that, while foraging, slugs do not move continuously but alternate periods of movement and rest; see the second column in Tables 1 and 2. Such behaviour is typical of many animal species38,40. In this section, we analyse the proportion of time that slugs spend moving, in particular to reveal the differences, if any, between the sparse and dense release.
Figure 6 shows the corresponding data where for the convenience of analysis the slugs are renumbered in a hierarchical order, so that slug 1 spends the highest proportion of time moving, slug 2 has the second highest, etc. We readily observe that the sparse release slugs tend to move more frequently than those from the dense release treatment: slugs that move for more than half of the total observation time constitute about 50% of the group in the case of sparse release but less than 30% in the case of dense release.
Figure 7

Distribution of the movement frequencies in case of (a) sparse release and (b) dense release.

Full size image

In order to make a more quantitative insight, we endeavour to describe the data using several standard distributions; see Tables 5 and 6. We find that the normal distribution performs better than others both in sparse and dense release treatments. Importantly, however, the parameters of the distribution are significantly different between the two cases; in particular, the standard deviation appears to be approximately twice as large in the case of sparse release. Arguably, it confirms the above conclusion that slugs move more frequently or for longer in the case of sparse release. Slugs released as a group tend to spend considerably more time at rest compared to the slugs released individually.
Table 5 The (r^2) values for the proportion of movement time described by different standard frequency distributions in the case of sparse release.
Full size table

To avoid a possible bias due to the different group size (17 slugs in the sparse release and 11 in the dense release), we now rearrange the data in terms of the proportion of the group that moves with a given frequency. The results are shown in Fig. 7. Although the amount of data in this case does not allow us to describe them using a particular function, the two cases clearly exhibit distributions with different properties. In particular, the average movement frequency is 0.467 for the sparse release and 0.264 for the dense release, and the corresponding variances are 0.090 and 0.065, respectively.
Table 6 The (r^2) values for the proportion of movement time described by different standard frequency distributions in the case of dense release.
Full size table

To further quantify the differences, Fig. 8 shows the number of slugs moving in each observation interval. Once again, we observe that the graph exhibits essentially different properties between the two releases. In particular, over the first interval, the majority of slugs (14 out of 17) move in the case of sparse release but none of the slugs move in the case of dense release. In the second half of the observation time (intervals 6–10) on average about 50% of slugs (8 out of 17) move in the case of sparse release but only about 25% of slugs (2–3 out of 11) move in the case of dense release.
Based on the differences between the two releases, we conclude that the presence of con-specifics is the factor that affects the distribution of slug movement time. Thus, along with the results of the previous sections, it suggests that slug movement is density dependent.
Figure 8

The number of moving slugs at each observation moment in case of (a) sparse release and (b) dense release.

Full size image More

Definitions of terms
GPP is defined here as ‘the sum of gross carbon fixation (carboxylation minus photorespiration) by autotrophic carbon-fixing tissues per unit area and time54. GPP is expressed as mass of organic carbon produced per unit area and time, over at least one year. NPP consists of all organic carbon that is fixed, but not respired over a given time period54:

$${mathrm{NPP}} = {mathrm{GPP}}-R_{mathrm{a}} = {Delta}B + L + F + H + O = {mathrm{BP}} + O$$
(3)

with all terms expressed in unit of mass of carbon per unit area and time. Ra is autotrophic respiration (composed of growth and maintenance respiration components); ΔB is the annual change in standing biomass carbon; litter production (roots, leaves and woody debris) is L; fruit production is F; the loss to herbivores is H, which was not accounted here because of the very limited number of observations available. BP is biomass production4. Symbol O represents occult, carbon flows, i.e. all other allocations of assimilated carbon, including changes in the nonstructural carbohydrate pool, root exudates, carbon subsidies to symbiotic fungi (mycorrhizae) or bacteria (e.g. nitrogen fixers), and BVOCs emissions (Supplementary Fig. 1). These ‘occult’ components are often ignored or unaccounted when estimating NPP, hence this bias is necessarily propagated into the Ra estimate when Ra is calculated as the difference between GPP and NPP55.
Estimation methods
We grouped the ‘methods’ into four categories:

biometric: direct tree stock measurements, or proxy data together with biomass expansion factors, allometric equations and the stock change as a BP component. If not otherwise stated, we assumed that the values included both above- and below-ground plant parts (n = 13 for GPP; n = 200 for NPP or BP).

micrometeorological: micrometeorological flux measurements using the eddy-covariance technique to measure CO2 flux and partitioning methods to estimate ecosystem respiration and GPP (n = 98 for GPP; n = 4 for NPP or BP).

model: model applications ranging from single mathematical equations (for canopy photosynthesis and whole-tree respiration) to more complex mechanistic process-based models to estimate GPP and Ra, with NPP as the net difference between them (n = 53 for GPP; n = 24 for NPP or BP).

scaling: upscaling of chamber-based measurements of assimilation and respiration (GPP and Ra) fluxes at the organ scale, or the entire stand (n = 73 for GPP; n = 9 for NPP or BP).

The difference between ‘scaling’ and ‘modelling’ lies in the data used. In the case of ‘scaling’ the data were derived from measurements at the site. ‘Model’ means that a dynamic process-based model was used, but with parameters calibrated and optimized at the site, based on either biometric or micrometeorological measurements.
Data selection
The data were obtained from more than 300 peer-reviewed articles (see also ref. 5), adding, merging and extending published works worldwide on CUE or BPE4,9,11,23,25,56,57. Data were extracted from the text, Tables or directly from Figures using the Unix software g3data (version 1.5.2, Jonas Frantz). In most studies, NPP, BP and GPP were estimated for the tree stand only. However, GPP estimated from CO2 flux by micrometeorological methods applies to the entire stand including ground vegetation. We therefore included only those micrometeorological studies where the forest stand was the dominant primary producer. The database contains 244 records (197 for BPE and 47 for CUE) from >100 forest sites (including planted, managed, recently burned, N-fertilized, irrigated and artificially CO2-fertilized forests; Supplementary Information, Supplementary Fig. 3 and online Materials; https://doi.org/10.5281/zenodo.3953478), representing 89 different tree species. Globally, 170 records out of the total data are from temperate sites, 51 from boreal, and 23 for tropical sites, corresponding to 79 deciduous broad-leaf (DBF), 14 evergreen broad-leaf (EBF), 132 evergreen needle-leaf (ENF) and 19 mixed-forests records (MX). The majority of the data (∼93%) cover the time-span from 1995 to 2015. We assume that when productivity data came from biometric measurements the reported NPP would have to be considered as BP because ‘occult’, nonstructural and secondary carbon compounds (e.g. BVOCs or exudates) are not included. In some cases, multiple datasets from the same site were included, covering different years or published by different authors. We considered only those values where either NPP (or BP) and GPP referred to the same year. From studies where data were available from more than 1 year, mean values across years were calculated. When the same reference for data was found in different papers or collected in different databases, where possible, we used data from the original source. When different authors described the same values for the same site, one single reference (and value) was used (in principle the oldest one). By using only commonly available environmental drivers to analyse the spatial variability in CUE and BPE, we were able to include almost all of the data that we found in the literature. We examined as potential predictors site-level effects of: average stand age (n = 204; range from 5 to ∼500 years), mean annual temperature (MAT; n = 230; range −6.5 to 27.1 °C) and total annual precipitation (TAP; n = 232; range from ∼125 to ∼3500 mm yr−1), method of determination (n = 237), geographic location (latitude and longitude; n = 241, 64°07′N to −42°52′S and 155°70′W to −173°28′E), elevation (n = 217; 5–2800 m, above sea level), leaf area index (LAI, n = 117; range from 0.4 to 13 m2 m−2), treatment (e.g.: ambient or artificially increased atmospheric CO2 concentration; n = 34), disturbance type (e.g.: fire n = 6; management n = 55), and the International Geosphere-Biosphere Programme (IGBP) vegetation classification and biomes (n = 244), as reported in the published articles (online Materials). The methods by which GPP, NPP, BP (and Ra) were determined were included as random effects in a number of possible mixed-effects linear regression models (Supplementary Table 4).
We excluded from statistical analysis all data where GPP and NPP were determined based on assumptions (e.g. data obtained using fixed fractions of NPP or Ra of GPP). In just one case GPP was estimated as the sum of upscaled Ra and NPP58; however, this study was excluded from the statistical analysis. NPP or Ra estimates obtained by process-based models (n = 23) were also not included in the statistical analysis. No information was available on prior natural disturbance events (biotic and abiotic, e.g. insect herbivore and pathogen outbreaks, and drought) that could in principle modify production efficiency, apart from fire. The occurrence of fire was reported by only a few studies59,60,61. These data were included in the database but fire, as an explanatory factor, was not considered due to the small number of samples in which it was reported (n = 6).
Data uncertainty
Uncertainties of GPP, NPP and BP data were all computed following the method based on expert judgment as described in Luyssaert et al.55. First, ‘gross’ uncertainty in GPP (gC m−2 yr−1) was calculated as 500 + 7.1 × (70−|lat|) gC m−2 yr−1 and gross uncertainties in NPP and BP (gC m−2 yr−1) were calculated as 350 + 2.9 × (70−|lat|). The absolute value of uncertainty thus decreases linearly with increasing latitude for GPP and for NPP and BP, because we assumed that the uncertainty is relative to the magnitude of the flux, which also decreases with increasing |lat|. Subsequently, as in Luyssaert et al.55, uncertainty was further reduced considering the methodology used to obtain each variable, by a method-specific factor (from 0 to 1, final uncertainty (δ) = gross uncertainty × method-specific factor). Luyssaert et al.55 reported for GPP-Micromet a method-specific factor of 0.3 (i.e. gross uncertainty is reduced by 70% for micrometeorological measurements); and for GPP-Model, 0.6. GPP-Scaling and GPP-Biometric were not explicitly considered in ref. 55 for GPP. We we used values of 0.8 and 0.3, respectively. For BP-Biometric and NPP-Micromet we used a reduction factor of 0.3; for NPP-Model, 0.6; and for NPP-Scaling (as obtained from chamber-based Ra measurements), 0.8. When GPP and/or NPP or BP methods were not known (n = 7), a factor of 1 (i.e. no reduction of uncertainty for methods used, hence maximum uncertainty) was used. The absolute uncertainties on CUE (δCUE) and BPE (δBPE) were considered as the weighted means62 by error propagation of each single variable (δNPP or δBP and δGPP) as follows:

$$delta {mathrm{CUE}} = sqrt {left( {frac{{delta {mathrm{NPP}}}}{{{mathrm{GPP}}}}} right)^2 + left( {delta {mathrm{GPP}}frac{{{mathrm{NPP}}}}{{{mathrm{GPP}}^2}}} right)^2}$$
(4)

and similarly for δBPE, by substituting NPP with BP and CUE with BPE.
Data and model selection
The CUE and BPE data were combined into a single variable, as sites for which both types of estimates existed did not show any significant differences between these entities (Supplementary Fig. 2). CUE values based on modelling were excluded (in our database we do not have BPE data from modelling). Tests showed that the CUE value was systematically higher when GPP was estimated with micrometeorological methods, compared to values based on biometric or scaling methods. Only data with complete information on CUE, MAT, age, TAP, and latitude were used. Altogether, 142 observations were selected.
In order to use the most complete information possible, a full additive model was constructed first (Eq. (1)). The method used for estimation of GPP (GPPmeth) was specified as a random effect on the intercept, as visual inspection suggested that CUE values were smaller where ‘scaling’ was used to estimate GPP compared to cases where ‘micromet’ was used to estimate GPP.
In Eq. (1) the variable ‘age’ represents the development status of the vegetation, i.e. either average age of the canopy forming trees or the period since the last major disturbance. The other three parameters represent different aspects of the climate. The absolute latitude, |lat|, was chosen as a proxy of radiation climate, i.e. day length and the seasonality of daily radiation. The term ηZGPPMeth represents the random effect on the intercept due to the different methods of estimating GPP.
These variables were not independent (Supplementary Table 1). If the different driver variables contain information that is not included in any of the other driver variables, multiple linear regression is nonetheless able to separate the individual effects. If, on the contrary, two variables exert essentially the same effect on the response variable (CUE) this can be seen in an ANOVA based model comparison. These considerations led us to the selection procedure in which we started with the full model (Eq. (1)) and compared it with all possible reduced models (Supplementary Table 2). The result of this analysis is the model with the smallest number of parameters that does not significantly differ from the full model.
We also examined, whether there were any significant interactions of predictor variables. There were not.
We used the R function lmer from the R-package lme463 to fit the mixed and ordinary multiple linear models to the data. We checked for potential problems of multicollinearity using the variance inflation factor (VIF)64. All predictors had VIF  More

In this study we analysed data from monthly ASP episodes over a 25 year period to determine if there was a seasonal pattern of presentation. We found a seasonal variation in scrotal torsion events (TT and TA), with higher frequency in the colder months, and an inverse correlation between monthly frequency and ambient temperature. There was increased frequency of EO in March, May and October but no correlation with temperature.
There has been interest in the seasonality of TT for many years6. Decreasing temperature causes increased contractility of the cremasteric muscles4,5, which may lead to an increase in the frequency of TT in colder months. Previous studies have yielded conflicting results3. Most of these studies reported case series of a relatively small number of patients (n = 39 to n = 2876). A large study from Brazil3 analysing 21,289 episodes of TT found seasonality of presentation with higher incidence in colder months, which was more significant in the more temperate regions than tropical regions of Brazil. A previous report from Dundee in Scotland showed an increased frequency of TT during the colder months from a series of 173 patients1. The present report is a larger study involving 33,855 episodes, of which 7882 had torsion events, and provides more robust evidence of seasonality of torsion.
Seasonal variation in the frequency of EO has not previously been reported to our knowledge. We were unable to explain the increased frequency of EO in March, May and October within our dataset. Further epidemiological study will be required to elucidate the reasons. Possibilities to consider include sexual behaviour patterns of the male population.
Limitations of this study include the use of data from a large public database with well reported advantages and disadvantages7, and the ecological fallacy, meaning that it may not be appropriate to apply these generalised population-based findings to individual patient care.
We do not suggest, based on our findings, that the threshold for surgical exploration be raised for patients with ASP presenting during warmer months. Public health measures could be considered, for example encouraging the wearing of warm clothing and undergarments by young males during colder months may reduce the frequency of TT and TA, as the style of clothing could have a direct effect on scrotal temperature8.
In conclusion, the findings of this large ecological study provide further robust evidence of seasonality of ASP, with the frequency of torsion events correlating negatively with ambient temperature. Further study is required to explain monthly variations in presentation of EO. 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