Ecology

Maxwell, S. L., Fuller, R. A., Brooks, T. M. & Watson, J. E. M. Biodiversity: The ravages of guns, nets and bulldozers. Nature 536, 143–145 (2016).ADS 
CAS 
PubMed 

Google Scholar 
Uchida, K. & Ushimaru, A. Biodiversity declines due to abandonment and intensification of agricultural lands: Patterns and mechanisms. Ecol. Monogr. 84, 637–658 (2014).
Google Scholar 
Habel, J. C. et al. Butterfly community shifts over two centuries: Shifts in butterfly communities. Conserv. Biol. 30, 754–762 (2016).PubMed 

Google Scholar 
Seibold, S. et al. Arthropod decline in grasslands and forests is associated with landscape-level drivers. Nature 574, 671–674 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Hallmann, C. A. et al. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS One 12, e0185809 (2017).PubMed 
PubMed Central 

Google Scholar 
Wenzel, M., Schmitt, T., Weitzel, M. & Seitz, A. The severe decline of butterflies on western German calcareous grasslands during the last 30 years: A conservation problem. Biol. Cons. 128, 542–552 (2006).
Google Scholar 
Biesmeijer, J. C. et al. Parallel declines in pollinators and insect-pollinated plants in Britain and the Netherlands. Science 313, 351–354 (2006).ADS 
CAS 
PubMed 

Google Scholar 
Hallmann, C. A., Foppen, R. P. B., van Turnhout, C. A. M., de Kroon, H. & Jongejans, E. Declines in insectivorous birds are associated with high neonicotinoid concentrations. Nature 511, 341–343 (2014).ADS 
CAS 
PubMed 

Google Scholar 
Møller, A. P. Parallel declines in abundance of insects and insectivorous birds in Denmark over 22 years. Ecol. Evol. 9, 6581–6587 (2019).PubMed 
PubMed Central 

Google Scholar 
Wagner, D. L. Insect declines in the anthropocene. Annu. Rev. Entomol. 65, 457–480 (2020).CAS 
PubMed 

Google Scholar 
Habel, J. C., Samways, M. J. & Schmitt, T. Mitigating the precipitous decline of terrestrial European insects: Requirements for a new strategy. Biodivers. Conserv. 28, 1343–1360 (2019).
Google Scholar 
Uhl, B., Wölfling, M. & Fiedler, K. Understanding small-scale insect diversity patterns inside two nature reserves: The role of local and landscape factors. Biodivers. Conserv. 29, 2399–2418 (2020).
Google Scholar 
Stevens, C. J., Dise, N. B., Mountford, J. O. & Gowing, D. J. Impact of nitrogen deposition on the species richness of grasslands. Science 303, 1876–1879 (2004).ADS 
CAS 
PubMed 

Google Scholar 
Thomas, J. A. Butterfly communities under threat. Science 353, 216–218 (2016).ADS 
CAS 
PubMed 

Google Scholar 
Sanders, J. & Hess, J. Benefits of organic farming to environment and society. Thünen Report 65, 362 (2019).
Google Scholar 
Brühl, C. A. & Zaller, J. G. Biodiversity decline as a consequence of an inappropriate environmental risk assessment of pesticides. Front. Environ. Sci. 7, 177 (2019).
Google Scholar 
Brühl, C. A. et al. Direct pesticide exposure of insects in nature conservation areas in Germany. Sci. Rep. 11, 24144 (2021).ADS 
PubMed 
PubMed Central 

Google Scholar 
Wagner, D. L., Grames, E. M., Forister, M. L., Berenbaum, M. R. & Stopak, D. Insect decline in the Anthropocene: Death by a thousand cuts. Proc. Natl. Acad. Sci. USA 118, e2023989118 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
Den Boer, P. J. & van Dijk, T. S. Carabid Beetles in A Changing Environment (Agricultural Univ, 1995).
Google Scholar 
Cristescu, M. E. From barcoding single individuals to metabarcoding biological communities: Towards an integrative approach to the study of global biodiversity. Trends Ecol. Evol. 29, 566–571 (2014).PubMed 

Google Scholar 
Hausmann, A. et al. Toward a standardized quantitative and qualitative insect monitoring scheme. Ecol. Evol. 10, 4009–4020 (2020).PubMed 
PubMed Central 

Google Scholar 
Ratnasingham, S. & Hebert, P. D. N. A DNA-based registry for all animal species: The Barcode Index Number (BIN) system. PLoS One 8, e66213 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Hausmann, A. et al. Genetic patterns in european geometrid moths revealed by the Barcode Index Number (BIN) system. PLoS One 8, e84518 (2013).ADS 
PubMed 
PubMed Central 

Google Scholar 
Padial, J. M., Miralles, A., De la Riva, I. & Vences, M. The integrative future of taxonomy. Front. Zool. 7, 1–14 (2010).
Google Scholar 
Schlick-Steiner, B. C. et al. Integrative taxonomy: A multisource approach to exploring biodiversity. Ann. Rev. Entomol. 55, 421–438 (2010).CAS 

Google Scholar 
Schlick‐Steiner, B. C., Arthofer, W., & Steiner, F. M. Take up the challenge! Opportunities for evolution research from resolving conflict in integrative taxonomy (2014).Fujita, M. K., Leaché, A. D., Burbrink, F. T., McGuire, J. A. & Moritz, C. Coalescent-based species delimitation in an integrative taxonomy. Trends Ecol. Evol. 27, 480–488 (2012).PubMed 

Google Scholar 
Morinière, J. et al. A DNA barcode library for 5,200 German flies and midges (Insecta: Diptera) and its implications for metabarcoding-based biomonitoring. Mol. Ecol. Res. 19, 900–928 (2019).
Google Scholar 
Kortmann, M. et al. Arthropod dark taxa provide new insights into diversity responses to bark beetle infestations. Ecol. Appl. 32, e2516 (2022).PubMed 

Google Scholar 
Porter, T. M. & Hajibabaei, M. Automated high throughput animal CO1 metabarcode classification. Sci. Rep. 8, 1–10 (2018).
Google Scholar 
Boggs, C. L. & Inouye, D. W. A single climate driver has direct and indirect effects on insect population dynamics: Climate drivers of population dynamics. Ecol. Lett. 15, 502–508 (2012).PubMed 

Google Scholar 
Conrad, K. F., Fox, R. & Woiwod, I. P. Monitoring biodiversity: Measuring long-term changes in insect abundance. In Insect Conservation Biology (eds Stewart, A. J. A. et al.) 203–225 (CABI, 2007). https://doi.org/10.1079/9781845932541.0203.Chapter 

Google Scholar 
Flohre, A. et al. Agricultural intensification and biodiversity partitioning in European landscapes comparing plants, carabids, and birds. Ecol. Appl. Publ. Ecol. Soc. Am. 21, 1772–1781 (2011).
Google Scholar 
Emmerson, M. et al. How agricultural intensification affects biodiversity and ecosystem services. In Advances in Ecological Research, vol ***55 43–97 (Elsevier, 2016).
Google Scholar 
Segerer, A. H. & Rosenkranz, E. Das grosse Insektensterben: Was es Bedeutet und was Wir Jetzt tun Müssen (Oekom Verlag, 2019).
Google Scholar 
Batáry, et al. The former Iron Curtain still drives biodiversity-profit trade-offs in German agriculture. Nat. Ecol. Evol. 1, 1279–1284 (2017).PubMed 

Google Scholar 
Kuussaari, M. et al. Extinction debt: A challenge for biodiversity conservation. Trends Ecol. Evol. 24, 564–571 (2009).PubMed 

Google Scholar 
Birkhofer, K., Smith, H. G., Weisser, W. W., Wolters, V. & Gossner, M. M. Land-use effects on the functional distinctness of arthropod communities. Ecography 38, 889–900 (2015).
Google Scholar 
Tscharntke, T., Klein, A. M., Kruess, A., Steffan-Dewenter, I. & Thies, C. Landscape perspectives on agricultural intensification and biodiversity—ecosystem service management. Ecol. Lett. 8, 857–874 (2005).
Google Scholar 
Habel, J. C., Seibold, S., Ulrich, W. & Schmitt, T. Seasonality overrides differences in butterfly species composition between natural and anthropogenic forest habitats. Anim. Conserv. 21, 405–413 (2018).
Google Scholar 
Schmitt, T., Ulrich, W., Delic, A., Teucher, M. & Habel, J. C. Seasonality and landscape characteristics impact species community structure and temporal dynamics of East African butterflies. Sci. Rep. 11, 15103 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ssymank, A. et al. Praktische Hinweise und Empfehlungen zur Anwendung von Malaisefallen für Insekten in der Biodiversitätserfassung und im Monitoring. Entomol. Verein Krefeld 1, 1–12 (2018).
Google Scholar 
Elbrecht, V., Peinert, B. & Leese, F. Sorting things out: Assessing effects of unequal specimen biomass on DNA metabarcoding. Ecol. Evol. 7, 6918–6926 (2017).PubMed 
PubMed Central 

Google Scholar 
Elbrecht, V. & Steinke, D. Scaling up DNA metabarcoding for freshwater macrozoobenthos monitoring. Freshw. Biol. 64, 380–387 (2019).CAS 

Google Scholar 
Boetzl, F. A. et al. A multitaxa assessment of the effectiveness of agri-environmental schemes for biodiversity management. Proc. Natl. Acad. Sci. 118, 25 (2021).
Google Scholar 
Uhler, J. et al. Relationship of insect biomass and richness with land use along a climate gradient. Nat. Commun. 12, 1–9 (2021).
Google Scholar 
Leray, M. et al. A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: Application for characterizing coral reef fish gut contents. Front. Zool. 10, 34 (2013).PubMed 
PubMed Central 

Google Scholar 
Morinière, J. et al. Species identification in malaise trap samples by DNA barcoding based on NGS Technologies and a scoring matrix. PLoS One 11, e0155497 (2016).PubMed 
PubMed Central 

Google Scholar 
Rognes, T., Flouri, T., Nichols, B., Quince, C. & Mahé, F. VSEARCH: A versatile open source tool for metagenomics. PeerJ 4, e2584 (2016).PubMed 
PubMed Central 

Google Scholar 
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10 (2011).
Google Scholar 
Ondov, B. D., Bergman, N. H. & Phillippy, A. M. Interactive metagenomic visualization in a Web browser. BMC Bioinform. 12, 385 (2011).
Google Scholar  More

Goodale, E., Beauchamp, G. & Ruxton, G. D. Mixed-Species Groups of Animals: Behavior, Community Structure, and Conservation (Academic Press, 2017).
Google Scholar 
Krause, J. & Ruxton, G. D. Living in Groups (Oxford University Press, 2002).
Google Scholar 
Stensland, E., Angerbjorn, A. & Berggren, P. Mixed species groups in mammals. Mamm. Rev. 33, 205–223 (2003).
Google Scholar 
Anderson, T. M. et al. Landscape-scale analyses suggest both nutrient and antipredator advantages to Serengeti herbivore hotspots. Ecology 91, 1519–1529 (2010).PubMed 

Google Scholar 
Sinclair, A. R. E. Does interspecific competition or predation shape the African ungulate community? J. Anim. Ecol. 54, 899–918 (1985).
Google Scholar 
Kiffner, C., Kioko, J., Leweri, C. & Krause, S. Seasonal patterns of mixed species groups in large East African mammals. PLoS ONE 9, e113446 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Meise, K., Franks, D. W. & Bro-Jørgensen, J. Using social network analysis of mixed species groups in African savannah herbivores to assess how community structure responds to environmental change. Philos. Trans. R. Soc. B Biol. Sci. 374, 20190009 (2019).
Google Scholar 
de Boer, W. F. & Prins, H. H. T. Large herbivores that thrive mightily but eat and drink as friends. Oecologia 82, 264–274 (1990).ADS 
PubMed 

Google Scholar 
Beaudrot, L., Palmer, M. S., Anderson, T. M. & Packer, C. Mixed-species groups of Serengeti grazers: A test of the stress gradient hypothesis. Ecology. https://doi.org/10.1002/ecy.3163 (2020).Article 
PubMed 

Google Scholar 
He, Q., Bertness, M. D. & Altieri, A. H. Global shifts towards positive species interactions with increasing environmental stress. Ecol. Lett. 16, 695–706 (2013).PubMed 

Google Scholar 
Bertness, M. D. & Callaway, R. Positive interactions in communities. Trends Ecol. Evol. 9, 191–193 (1994).CAS 
PubMed 

Google Scholar 
Fugère, V. et al. Testing the stress-gradient hypothesis with aquatic detritivorous invertebrates: Insights for biodiversity-ecosystem functioning research. J. Anim. Ecol. 81, 1259–1267 (2012).PubMed 

Google Scholar 
Bakker, E. S., Dobrescu, I., Straile, D. & Holmgren, M. Testing the stress gradient hypothesis in herbivore communities: Facilitation peaks at intermediate nutrient levels. Ecology 94, 1776–1784 (2013).PubMed 

Google Scholar 
Hopcraft, J. G. C., Olff, H. & Sinclair, A. R. E. Herbivores, resources and risks: Alternating regulation along primary environmental gradients in savannas. Trends Ecol. Evol. 25, 119–128 (2010).PubMed 

Google Scholar 
Sih, A. Optimal behavior: Can foragers balance two conflicting demands? Science 210, 1041–1043 (1980).ADS 
CAS 
PubMed 

Google Scholar 
Creel, S. & Christianson, D. Relationships between direct predation and risk effects. Trends Ecol. Evol. 23, 194–201 (2008).PubMed 

Google Scholar 
Zollner, P. A. & Lima, S. L. Towards a behavioral ecology of ecological landscapes. Trends Ecol. Evol. 11, 131–135 (1996).PubMed 

Google Scholar 
Brown, J. S., Laundré, J. W. & Gurung, M. The ecology of fear: Optimal foraging, game theory, and trophic interactions. J. Mammal. 80, 385–399 (1999).
Google Scholar 
Gaynor, K. M., Brown, J. S., Middleton, A. D., Power, M. E. & Brashares, J. S. Landscapes of fear: Spatial patterns of risk perception and response. Trends Ecol. Evol. 34, 355–368 (2019).PubMed 

Google Scholar 
Creel, S., Schuette, P. & Christianson, D. Effects of predation risk on group size, vigilance, and foraging behavior in an African ungulate community. Behav. Ecol. 25, 773–784 (2014).
Google Scholar 
Goodale, E., Beauchamp, G., Magrath, R. D., Nieh, J. C. & Ruxton, G. D. Interspecific information transfer influences animal community structure. Trends Ecol. Evol. 25, 354–361 (2010).PubMed 

Google Scholar 
Freeberg, T. M., Eppert, S. K., Sieving, K. E. & Lucas, J. R. Diversity in mixed species groups improves success in a novel feeder test in a wild songbird community. Sci. Rep. 7, 43014 (2017).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Anderson, T. M. et al. The spatial distribution of african savannah herbivores: Species associations and habitat occupancy in a landscape context. Philos. Trans. R. Soc. B Biol. Sci. 371, 20150314 (2016).
Google Scholar 
Arsenault, R. & Owen-Smith, N. Resource partitioning by grass height among grazing ungulates does not follow body size relation. Oikos 117, 1711–1717 (2008).
Google Scholar 
Esmaeili, S. et al. Body size and digestive system shape resource selection by ungulates: A cross-taxa test of the forage maturation hypothesis. Ecol. Lett. 24, 2178–2191 (2021).PubMed 

Google Scholar 
Hopcraft, J. G. C., Anderson, T. M., Pérez-Vila, S., Mayemba, E. & Olff, H. Body size and the division of niche space: Food and predation differentially shape the distribution of Serengeti grazers. J. Anim. Ecol. 81, 201–213 (2012).PubMed 

Google Scholar 
McArthur, C., Banks, P. B., Boonstra, R. & Forbey, J. S. The dilemma of foraging herbivores: Dealing with food and fear. Oecologia 176, 677–689 (2014).ADS 
PubMed 

Google Scholar 
Gagnon, M. & Chew, A. E. Dietary preferences in extant African Bovidae. J. Mammal. 81, 490–511 (2000).
Google Scholar 
Kartzinel, T. R. et al. DNA metabarcoding illuminates dietary niche partitioning by African large herbivores. Proc. Natl. Acad. Sci. U.S.A. 112, 8019–8024 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Veldhuis, M. P. et al. Cross-boundary human impacts compromise the Serengeti-Mara ecosystem. Science 363, 1424–1428 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Kavwele, C. M. et al. Non-local effects of human activity on the spatial distribution of migratory wildlife in Serengeti National Park, Tanzania. Ecol. Solut. Evid. 3, e12159 (2022).
Google Scholar 
Bijlsma, R. & Loeschcke, V. Environmental stress, adaptation and evolution: An overview. J. Evol. Biol. 18, 744–749 (2005).CAS 
PubMed 

Google Scholar 
Schmitt, M. H., Stears, K. & Shrader, A. M. Zebra reduce predation risk in mixed-species herds by eavesdropping on cues from giraffe. Behav. Ecol. 27, 1073–1077 (2016).
Google Scholar 
Preisser, E. L., Orrock, J. L. & Schmitz, O. J. Predator hunting mode and habitat domain alter nonconsmuptive effects in predator-prey interactions. Ecology 88, 2744–2751 (2007).PubMed 

Google Scholar 
Kiffner, C. et al. Long-term persistence of wildlife populations in a pastoral area. Ecol. Evol. 10, 10000–10016 (2020).PubMed 
PubMed Central 

Google Scholar 
Hopcraft, J. G. C. et al. Competition, predation, and migration: Individual choice patterns of Serengeti migrants captured by hierarchical models. Ecol. Monogr. 84, 355–372 (2014).
Google Scholar 
Fryxell, J. M. Forage quality and aggregation by large herbivores. Am. Nat. 138, 478–498 (1991).
Google Scholar 
Fitzgibbon, C. D. Mixed-species grouping in Thomson’s and Grant’s gazelles: The antipredator benefits. Anim. Behav. 39, 1116–1126 (1990).
Google Scholar 
Brown, J. S. & Kotler, B. P. Hazardous duty pay and the foraging cost of predation. Ecol. Lett. 7, 999–1014 (2004).
Google Scholar 
Stears, K. & Shrader, A. M. Increases in food availability can tempt oribi antelope into taking greater risks at both large and small spatial scales. Anim. Behav. 108, 155–164 (2015).
Google Scholar 
Creel, S. Toward a predictive theory of risk effects: Hypotheses for prey attributes and compensatory mortality. Ecology 92, 2190–2195 (2011).PubMed 

Google Scholar 
Périquet, S. et al. Effects of lions on behaviour and endocrine stress in plains zebras. Ethology 123, 667 (2017).
Google Scholar 
Stears, K., Schmitt, M. H., Wilmers, C. C. & Shrader, A. M. Mixed-species herding levels the landscape of fear. Proc. R. Soc. B Biol. Sci. 287, 20192555 (2020).
Google Scholar 
Schmitt, M. H., Stears, K., Wilmers, C. C. & Shrader, A. M. Determining the relative importance of dilution and detection for zebra foraging in mixed-species herds. Anim. Behav. 96, 151–158 (2014).
Google Scholar 
Meise, K., Franks, D. W. & Bro-Jørgensen, J. Alarm communication networks as a driver of community structure in African savannah herbivores. Ecol. Lett. 23, 293–304 (2020).PubMed 

Google Scholar 
Codron, D., Hofmann, R. R. & Clauss, M. Morphological and physiological adaptations for browsing and grazing. In The Ecology of Browsing and Grazing II (eds Gordon, I. J. & Prins, H. H. T.) 81–125 (Springer, 2019).
Google Scholar 
Odadi, W. O., Karachi, M. K., Abdulrazak, S. A. & Young, T. P. African wild ungulates compete with or facilitate cattle depending on season. Science 333, 1753–1755 (2011).ADS 
CAS 
PubMed 

Google Scholar 
Maestre, F. T., Callaway, R. M., Valladares, F. & Lortie, C. J. Refining the stress-gradient hypothesis for competition and facilitation in plant communities. J. Ecol. 97, 199–205 (2009).
Google Scholar 
de Jonge, M. M. J. et al. Conditional love? Co-occurrence patterns of drought-sensitive species in European grasslands are consistent with the stress-gradient hypothesis. Glob. Ecol. Biogeogr. 30, 1609–1620 (2021).PubMed 
PubMed Central 

Google Scholar 
Franks, D. W., Weiss, M. N., Silk, M. J., Perryman, R. J. Y. & Croft, D. P. Calculating effect sizes in animal social network analysis. Methods Ecol. Evol. 12, 33–41 (2021).
Google Scholar 
Estes, J. A. et al. Trophic downgrading of planet earth. Science 333, 301–306 (2011).ADS 
CAS 
PubMed 

Google Scholar 
Meise, K., Franks, D. W. & Bro-Jørgensen, J. Multiple adaptive and non-adaptive processes determine responsiveness to heterospecific alarm calls in African savannah herbivores. Proc. R. Soc. B Biol. Sci. 285, 20172676 (2018).
Google Scholar 
Blumstein, D. T., Bitton, A. & DaVeiga, J. How does the presence of predators influence the persistence of antipredator behavior? J. Theor. Biol. 239, 460–468 (2006).ADS 
MathSciNet 
PubMed 
MATH 

Google Scholar 
Riggio, J. et al. Lion populations may be declining in Africa but not as Bauer et al. suggest. Proc. Natl. Acad. Sci. 113, 201521506 (2015).
Google Scholar 
Bauer, H. et al. Lion (Panthera leo) populations are declining rapidly across Africa, except in intensively managed areas. Proc. Natl. Acad. Sci. 112, 14894–14899 (2015).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Pettorelli, N., Bro-Jørgensen, J., Durant, S. M., Blackburn, T. & Carbone, C. Energy availability and density estimates in African ungulates. Am. Nat. 173, 698–704 (2009).PubMed 

Google Scholar 
Haile, G. G. et al. Projected impacts of climate change on drought patterns over East Africa. Earth’s Future 8, 1–23 (2020).
Google Scholar 
Devine, A. P., McDonald, R. A., Quaife, T. & Maclean, I. M. D. Determinants of woody encroachment and cover in African savannas. Oecologia 183, 939–951 (2017).ADS 
PubMed 
PubMed Central 

Google Scholar 
Kiffner, C. et al. Long-term population dynamics in a multi-species assemblage of large herbivores in East Africa. Ecosphere 8, e02027 (2017).
Google Scholar 
Prins, H. H. T. & Loth, P. E. Rainfall patterns as background to plant phenology in northern Tanzania. J. Biogeogr. 15, 451–463 (1988).
Google Scholar 
Beattie, K., Olson, E. R., Kissui, B., Kirschbaum, A. & Kiffner, C. Predicting livestock depredation risk by African lions (Panthera leo) in a multi-use area of northern Tanzania. Eur. J. Wildl. Res. 66, 11 (2020).
Google Scholar 
Kasozi, H. & Montgomery, R. A. Variability in the estimation of ungulate group sizes complicates ecological inference. Ecol. Evol. 10, 6881–6889 (2020).PubMed 
PubMed Central 

Google Scholar 
USGS. MOD13Q1 v006 MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid. 10.5067/MODIS/MOD13Q1.006 (2020).R Core Team. R: A Language and Environment for Statistical Computing. http://www.r-project.org/. Accessed January 02, 2022 (2021).Dice, L. R. Measures of the amount of ecologic association between species. Ecology 26, 297–302 (1945).
Google Scholar 
Croft, D. P., James, R. & Krause, J. Exploring Animal Social Networks (Princeton University Press, 2008).
Google Scholar 
Besag, J. & Clifford, P. Generalized Monte Carlo significance tests. Biometrika 76, 633–642 (1989).MathSciNet 
MATH 

Google Scholar 
Hayward, M. W. & Kerley, G. I. H. Prey preferences of the lion (Panthera leo). J. Zool. 267, 309–322 (2005).
Google Scholar 
Codron, D. et al. Diets of savanna ungulates from stable carbon isotope composition of faeces. J. Zool. 273, 21–29 (2007).
Google Scholar 
Kartzinel, T. R. & Pringle, R. M. Multiple dimensions of dietary diversity in large mammalian herbivores. J. Anim. Ecol. 89, 1482–1496 (2020).PubMed 

Google Scholar 
Prins, H. H. T. & Douglas-Hamilton, I. Stability in a multi-species assemblage of large herbivores in East Africa. Oecologia 83, 392–400 (1990).ADS 
CAS 
PubMed 

Google Scholar 
Tournier, E. et al. Differences in diet between six neighbouring groups of vervet monkeys. Ethology 120, 471–482 (2014).
Google Scholar 
Humphries, B. D., Ramesh, T. & Downs, C. T. Diet of black-backed jackals (Canis mesomelas) on farmlands in the KwaZulu-Natal Midlands, South Africa. Mammalia 80, 405–412 (2016).
Google Scholar  More

Thermogravimetric/differential thermogravimetry analyses of corn stoverThermogravimetric/Differential Thermogravimetry (TG/DTG) curves are shown in Fig. 2. The pyrolysis process of corn stover could be divided into three stages. The first stage was the dehydration stage, which occurred at approximately 55–125 °C, and the weight loss was mainly accounted for by water19. The second stage was the pyrolysis stage, which occurred at approximately 200–400 °C and mainly involved the decomposition of cellulose, hemicellulose and a small amount of lignin. This process involved the generation of CO and CO2 and the breaking of carbonaceous polymer bonds20. In addition, a shoulder peak in the range of 265 to 300 °C in the DTG diagram could be caused by side chain decomposition and glycosidic bond cleavage of xylan during the pyrolysis of corn stover21. The third stage was the carbonization stage, which occurred above 400 °C; this stage mainly involved the decomposition of lignin22,23. The carbonization process was relatively slow after 600 °C; this process was called the passive pyrolysis stage24. In general, the TG loss in the pyrolysis process of corn stover was mainly from the moisture in the biomass sample in the first stage. Hemicellulose and cellulose decomposition occurred in the second stage, and lignin decomposition occurred in the third stage25. In this experiment, the minimum pyrolysis temperature for the preparation of biochar was 400 °C. Therefore, the pyrolysis of biochar was relatively complete.Figure 2TG/DTG curves of corn stover.Full size imageCharacterization of biocharYield and specific surface area analysesThe yield and SBET are presented in Table 2. BC, BC-H and BC-OH represent the origin, acid-modified, and base-modified biochar, respectively. The yield of corn stover biochar exhibited a negative correlation with the temperature and decreased from 39.65 to 28.26% when the pyrolysis temperature increased from 400 to 700 °C. This phenomenon could have occurred due to the loss of more volatile substances and the thermal degradation of lignocellulose with increasing temperature, thus reducing the yield of biochar26,27. The SBET of the original biochar showed little difference below 700 °C but increased significantly at 700 °C. Combined with the SEM analysis (Fig. 3), at low temperatures, more ashes on the surface of biochar could block its pores so that the change in SBET was not obvious. At 700 °C, because the ash content significantly reduced and the pyrolysis was more sufficient, the pores of the biochar were more developed, and the SEBT significantly increased. The SBET of the acid/base-modified biochar increased with increasing temperature. The SBET of biochar was larger than that of the original biochar after acid and base modification at 400–600 °C. This phenomenon occurred because the porous structure of biochar was enhanced by acid and base modification28. Moreover, pickling removed most of the inorganic substances in biochar and reduced ash content, while alkali washing removed the tar on the surface of biochar to a certain extent29. However, at 700 °C, the SBET of biochar after acid/base modification was lower than that of the original biochar. Combined with the SEM (Fig. 4), the acid/base modification caused the nanopores of biochar to collapse into mesopores or macropores30. Therefore, the well-developed pore structure of the biochar prepared at 700 °C was destroyed by acid/base modification, resulting in a significant decrease in SBET.Table 2 Yield and SBET of different biochars.Full size tableFigure 3SEM (ZEISS) images of biochar at different pyrolysis temperatures: (a) C1, (b) C8, (c) C12, and (d) C16.Full size imageFigure 4SEM (OPTON) images of C16 biochar and its acid/base modification: (a) C16, (b) C16-H, and (c) C–OH.Full size imageScanning electron microscopy analysisThe C1, C8, C12 and C16 biochars had the highest Cd2+ removal rates at 400, 500, 600 and 700 °C, respectively. Therefore, these BCs were selected for SEM analysis. Figure 3 clearly showed that as the pyrolysis temperature increased from 400 to 700 °C, the pore structure of biochar became more developed, with a smaller pore size and more pores. Although there were numerous pores at 500 °C, the pores were not fully developed and were blocked inside. At 700 °C, the skeleton structure appeared, and the particle size of ash blocked in the pores decreased.By taking C16 biochar with the highest removal rate of Cd2+ as the research object, the changes in the biochar surface before and after modification were compared. C16-H and C16-OH represent acid-modified and base-modified biochar, respectively. After acid/base modification, the ash content on the surface of the biochar decreased, and the pore size increased (Fig. 4). Therefore, some skeleton structures could collapse after corrosion, which was consistent with the previous SBET results. Sun et al. discovered that citric acid-modified biochar would lead to micropore wall collapse and micropore loss, resulting in a reduction in SBET31. This finding was in agreement with the results of our study.Fourier transform infrared spectroscopy analysisThe FTIR spectra of biochar at different pyrolysis temperatures are presented in Fig. 5a.Figure 5FTIR spectra of corn stover biochar: (a) different pyrolysis temperatures and (b) different modification treatments.Full size imageAs the pyrolysis temperature increased from 300 to 700 °C, the absorption peak intensity showed a downwards trend. There was a remarkable decrease in features associated with stretch O–H (3400 cm−1)32. The vibration peaks of C–H (2924 cm−1) and C=O (1610 cm−1) decreased with increasing temperature, which could be due to the reduction in –CH2 and –CH3 groups of small molecules and the pyrolysis of C=O into gas or liquid byproducts at high temperatures33. In addition, the peak at 1435 cm−1 was identified as the vibration of C=C bonds belonging to the aromatic skeleton of biochar. A decrease in the absorbance peaks was found at 1115 cm−1, which corresponded to C–O–C bonds. The ratio of intensities for C=C/C=O (1550–1650 cm−1) and C–O–C (1115 cm−1) to the shoulder (1100–1200 cm−1) gradually decreased, and the loss of –OH at 3444 cm−1 indicated that the oxygen content in biochar reduced. The cellulose and wood components were dehydrated, and the degree of biochar condensation increased at higher temperatures. The bending vibration peaks of Ar–H at 856 and 877 cm−1 changed little at different temperatures, which showed that the aromatic rings were relatively stable below 700 °C34. Combined with the above analysis the condensation degree of biochar increased gradually above 400 °C35,36. In summary, as the pyrolysis temperature increased, the degree of aromatization of biochar improved, and the numbers of oxygen-containing functional groups decreased continuously.Figure 5b showed that after acid/base modification, the absorbance peaks at 3444 cm−1, 1610 cm−1 and 1115 cm−1 increased, indicating that the number of oxygen-containing functional groups increased. However, the stretching vibration peak of aromatic ring skeleton C=C (1435 cm−1) and the bending vibration peaks of Ar–H (856–877 cm−1) changed little. The number of functional groups of acid-modified biochar increased more than that of alkali-modified biochar. Mahdi et al. found that acid modification increased the number of functional groups in a study of biochar modification37. After acid/base modification, the number of oxygen-containing functional groups, such as hydroxyl and carboxyl groups, increased.Optimization of biocharFigure 6 illustrates that the removal rates of Cd2+ by corn stover biochar (original, acid-modified, and base-modified biochars) consistently increased with increasing pyrolysis temperature. The highest removal rate reached 95.79% at 700 °C. The removal rate decreased after modification, especially after pickling. The results showed that C16 biochar had the best removal effect on Cd2+.Figure 6Cd2+ removal rate of different biochars (BC: original biochar, BC-OH: alkali-modified biochar, and BC-H: acid-modified biochar).Full size imageIntuitive and variance analyses were employed to explore the influences of biochar preparation conditions on the removal rate of Cd2+.

1.

Intuitive analysis
The intuitive analysis of the orthogonal experiment is shown in Table 3 and Fig. 7. The pyrolysis temperature had the most significant influence on the removal of Cd2+, followed by the retention time and finally the heating rate. Therefore, the optimal conditions for biochar preparation were a pyrolysis temperature of 700 °C, a retention time of 2.5 h, and a heating rate of 5 °C/min.

2.

Variance analysis
Variance analysis showed that the effect of pyrolysis temperature on the removal rate of Cd2+ was very significant (Table 4). The effects of retention time and heating rate were not significant. This phenomenon was consistent with the conclusions obtained in the intuitive analysis.

Table 3 Intuitive analyses of influencing factors of biochar preparation.Full size tableFigure 7Intuitive analysis diagram of influencing factors for biochar preparation.Full size imageTable 4 Variance analysis.Full size tableAnalysis of adsorption mechanismThe SBET of the unmodified biochar did not change significantly with temperature, which indicated that SBET could potentially not be a critical factor for Cd2+ adsorption. Qi et al. obtained a similar conclusion when studying the adsorption of Cd2+ in water by chicken litter biochar38. In addition to SBET, the four primary mechanisms involved in the removal of heavy metal ions by biochar were as follows: (1) Ion exchange: the alkali or alkaline earth metals in biochar (K+, Ca2 +, Na+, and Mg2+) were the dominant cations in ion exchange39. (2) The complexation of oxygen-containing functional groups mainly included hydroxyl and carboxyl groups40. (3) Mineral precipitation: Cd2+ was precipitated by minerals on the surface of biochar to form Cd3(PO4)2 and CdCO341. Soluble cadmium precipitated with some anions released by biochar, such as CO32−, PO43− and OH−42,43. (4) π electron interaction: Cd2+ coordinated with the π electrons of C=C or C=O at low pyrolysis temperatures43,44. Biochar contains more aromatic structures at high pyrolysis temperatures, which could provide more π electrons. Therefore, the π electron interaction in adsorption of Cd2+ was effectively enhanced45.C1, C8, C12 and C16 were selected to study the adsorption mechanism. Related physicochemical properties are given in Table 5.Table 5 Physicochemical properties of biochar at different pyrolysis temperatures.Full size tableThe CEC of biochar gradually increased as the pyrolysis temperature increased, reaching a maximum at 600 °C and slightly decreasing at 700 °C. This phenomenon could have occurred because the crystalline minerals under high pyrolysis temperatures inhibited the exchange of cations on the surface of biochar with Cd2+ in aqueous solution46. Nevertheless, CEC did not change significantly with temperature; thus, CEC was not the main adsorption mechanism. With increasing pyrolysis temperature, the number of acidic functional groups decreased gradually, while the number of alkaline functional groups increased. The main functional groups used to remove Cd2+ were generally considered acidic oxygen-containing functional groups. However, the number of these functional groups decreased with increasing pyrolysis temperature, which weakened the complexation on the surface of the biochar. However, this result was contradictory to the results of Cd2+ adsorption. Therefore, the functional groups were not the main adsorption mechanism.To further explore the adsorption mechanism of Cd2+, the biochar before and after the adsorption of Cd2+ was characterized by XRD. As shown in Fig. 7a, C16-100Cd and C16-200Cd represented the biochar after Cd2+ adsorption when the concentrations of cadmium solution were 100 mg/l and 200 mg/l, respectively. The results showed that new peaks appeared at 30.275° and 36.546° after adsorption, corresponding to CdCO3. The spike at 29.454° was due to Cd(OH)2. Additionally, the intensity of the CdCO3 peak increased significantly from C16-100Cd to C16-200Cd, indicating that mineral precipitation occurred in adsorption. Liu et al. found similar results in a study on removing Cd2+ from water by blue algae biochar12. However, as the concentration of Cd2+ increased from 0 to 200 mg/L, the diffraction peak at 2θ = 29.454° first increased and then decreased. This because the peak position of CaCO3 at 2θ = 29.369° was very close to Cd(OH)2 at 2θ = 29.454°. At low concentrations, the production of Cd(OH)2 was greater than that of CdCO3. When the initial concentration of Cd2+ increased, more CO32− released by CaCO3 combined with Cd2+ to form CdCO3, resulting in a reduction in the diffraction peak.As presented in Fig. 8b, the peak intensities of CdCO3 and Cd(OH)2 gradually increase with increasing pyrolysis temperature. On the one hand, this phenomenon could be ascribed to the increase in the mineral content of biochar with increasing pyrolysis temperature. On the other hand, the pH value of biochar increased with increasing pyrolysis temperature. In this way, more OH− was released, thus forming more Cd(OH)2. Wang et al. obtained similar results42. Moreover, the peak intensity of KCl at 2θ = 28.347° decreased after adsorption, as shown in Fig. 8a, which indicated that ion exchange took part in adsorption.Figure 8XRD images: (a) before and after adsorption of Cd2+ on C16 biochar and (b) Cd2+ adsorption by biochar at different pyrolysis temperatures.Full size imageIn addition, the FTIR spectra showed that the number of functional groups, such as C=C and C=O, in biochar decreased with increasing pyrolysis temperature, leading to the weakening of cation–π interactions between Cd2+ and C=C and C=O. In contrast, due to the enhanced aromatization of functional groups on the surface of biochar, many lone pair electrons existed in the electron-rich domains of the graphene-like structure, which in turn enhanced the cation–π interactions. Harvey et al., based on the study of Cd2+ adsorption by plant biochar, concluded that the electron-rich domain bonding mechanism between Cd2+ and the graphene-like structure on the surface of biochar played a more significant role in biochar with a high degree of carbonization45. Therefore, π-electron interactions could play a dominant role in Cd2+ adsorption on high-temperature pyrolysis biochar. Moreover, the results showed that the number of alkaline functional groups increased while acidic functional groups decreased with the increase in pyrolysis temperature. It is generally believed that acidic functional groups could withdraw electrons, and basic functional groups could donate electrons47,48. The biochar with higher pyrolysis temperature contained more alkaline functional groups, which improved the electron donating ability of biochar and enhanced the cation–π electron effect.In summary, mineral precipitation and π electron coordination were the main mechanisms of removing Cd2+ from water by corn stover biochar. This phenomenon explained why the Cd2+ removal rate of acid/base–modified biochar decreased. After modification, the functional groups on the surface of biochar increased, but the inorganic minerals were removed. Pickling resulted in the loss of soluble minerals and alkaline functional groups on the surface of biochar, which was not conducive to adsorption49. After alkaline washing, more PO43−, CO32− and HCO3− were released, thereby reducing the mineral precipitation50,51. Since NaOH had a weaker destructive effect than HCl and introduced some OH−, alkaline washing had little effect on the removal rate of Cd2+.Adsorption isotherm and adsorption kineticsAdsorption isothermThe adsorption isotherms were fitted with Langmuir (Eq. 3) and Freundlich (Eq. 4) models, as shown in Fig. 9, and the fitting parameters are listed in Table 6.Figure 9Adsorption isotherm.Full size imageTable 6 Fitting parameters of the adsorption isotherm model.Full size tableThe Langmuir model (R2  > 0.963) was more suitable than the Freundlich model (R2  > 0.919), indicating that the adsorption sites of biochar were evenly distributed, and adsorption was mainly monolayer. Parameter KL reflected the difficulty of adsorption and was generally divided into four types: unfavourable (KL  > 1), favourable (0  More

Buchmann, K. & Mehrdana, F. Effects of anisakid nematodes Anisakis simplex (s.l.), Pseudoterranova decipiens (s.l.) and Contracaecum osculatum (s.l.) on fish and consumer health. Food Waterborne Parasitol. 4, 13–22. https://doi.org/10.1016/j.fawpar.2016.07.003 (2016).Article 

Google Scholar 
Mattiucci, S., Cipriani, P., Levsen, A., Paoletti, M. & Nascetti, G. Molecular epidemiology of Anisakis and anisakiasis: An ecological and evolutionary road map. Adv. Parasitol. https://doi.org/10.1016/bs.apar.2017.12.001 (2018).Article 
PubMed 

Google Scholar 
Mattiucci, S., Cipriani, P., Paoletti, M., Levsen, A. & Nascetti, G. Reviewing biodiversity and epidemiological aspects of anisakid nematodes from the North-east Atlantic Ocean. J. Helminthol. https://doi.org/10.1017/S0022149X1700027X (2017).Article 
PubMed 

Google Scholar 
Moravec, F. & Justine, J.-L. Erection of Euterranova n. gen. and Neoterranova n. gen. (Nematoda, Anisakidae), with the description of E. dentiduplicata n. sp. and new records of two other anisakid nematodes from sharks off New Caledonia. Parasite 27, 58. https://doi.org/10.1051/parasite/2020053 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Shamsi, S. & Suthar, J. Occurrence of Terranova larval types (nematoda: Anisakidae) in Australian marine fish with comments on their specific identities. Peer J. 4, e1722. https://doi.org/10.7717/peerj.1722 (2016).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Timi, J. T. et al. Molecular identification, morphological characterization and new insights into the ecology of larval Pseudoterranova cattani in fishes from the Argentine coast with its differentiation from the Antarctic species, P. decipiens sp. E (Nematoda: Anisakidae). Vet. Parasitol. 199, 59–72 (2014).CAS 
PubMed 

Google Scholar 
Deardorff, T. L. Redescription of Pulchrascaris chiloscyllii (Johnston and Mawson, 1951) (Nematoda: Anisakidae), with comments on species in Pulchrascaris and Terranova. Proc. Helminthol. Soc. Wash. 54, 28–39 (1987).
Google Scholar 
Cannon, L. R. G. Some larval ascaridoids from south-eastern queensland marine fishes. Int. J. Parasitol. 7, 233–243 (1977).CAS 
PubMed 

Google Scholar 
Levsen, A. & Lunestad, B. T. Anisakis simplex third stage larvae in Norwegian spring spawning herring (Clupea harengus L.), with emphasis on larval distribution in the flesh. Vet. Parasitol. 171, 247–253 (2010).PubMed 

Google Scholar 
Berland, B., (1989) Identification of fish larval nematodes from fish. In: Möller H, editor. Nematode problems in North Atlantic fish. Report from a workshop in Kiel, 3 4 16–22.Zhu, X., D’Amelio, S., Paggi, L. & Gasser, R. B. Assessing sequence variation in the internal transcribed spacers of ribosomal DNA within and among members of the Contracaecum osculatum complex (nematoda: Ascaridoidea: Anisakidae). Parasitol. Res. 86, 677–683 (2000).CAS 
PubMed 

Google Scholar 
Nadler, S. A. & Hudspeth, D. S. S. Phylogeny of the ascaridoidea (Nematoda: Ascaridida) based on three genes and morphology hypotheses of structural and sequence evolution. J. Parasitol. 86, 380–393. https://doi.org/10.1645/0022-3395(2000)086[0380:POTANA]2.0.CO;2 (2000).Article 
CAS 
PubMed 

Google Scholar 
Mattiucci, S. et al. Genetic and morphological approaches distinguish the three sibling species of the Anisakis simplex species complex, with a species designation as Anisakis berlandi n. sp. for A simplex sp. C (Nematoda: Anisakidae). J. Parasitol. 100, 199–214. https://doi.org/10.1645/12-120.1 (2014).Article 
CAS 
PubMed 

Google Scholar 
Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).CAS 
PubMed 

Google Scholar 
Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35, 1547–1549 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Nagy, L. G. et al. Re-mind the gap! Insertion – deletion data reveal neglected phylogenetic potential of the nuclear ribosomal internal transcribed spacer (ITS) of fungi. PLoS ONE 7, 1–9 (2012).
Google Scholar 
Suchard, M. A. et al. Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 4, 1–5 (2018).
Google Scholar 
Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior summarization in bayesian phylogenetics using tracer 1.7. Syst. Biol. 67, 901–904 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Cipriani, P. et al. Anisakid nematodes in Trichiurus lepturus and Saurida undosquamis (Teleostea) from the South-West Indian Ocean : Genetic evidence for the existence of sister species within Anisakis typica (s.l.), and food-safety considerations. Food Waterborne Parasitol. 28, e00177. https://doi.org/10.1016/j.fawpar.2022.e00177 (2022).Article 
PubMed 
PubMed Central 

Google Scholar 
Safonova, A. E. First report on molecular identification of Anisakis simplex in Oncorhynchus nerka from the fish market, with taxonomical issues within Anisakidae. J. Nematol. 53(1), 10. https://doi.org/10.21307/jofnem-2021-023 (2021).Article 
CAS 

Google Scholar 
Takano, T. & Sata, N. Multigene phylogenetic analysis reveals non-monophyly of Anisakis s.l. and Pseudoterranova (Nematoda: Anisakidae). Parasitol. Int. 91, 102631. https://doi.org/10.1016/j.parint.2022.102631 (2022).Article 
CAS 
PubMed 

Google Scholar 
Leiper, R. T. & Atkinson, E. L. Parasitic worms, with a note on a free-living nematode. British Museum (Natural History). Bristish Antarctic (“Terra Nova”) expedition, 1910. Natural History Report. Zool 2(3), 19–60 (1915).
Google Scholar 
Leiper, R. T. & Atkinson, E. L. Helminthes of the British Antarctic expedition 1910–1913. Proc. Zool. Soc. London, 222–226 (1914).Myers, B. J. Phocanema, a new genus for the anisakid nematode of seals. Can. J. Zool. 37, 459–465 (1959).
Google Scholar 
Mattiucci, S., Paoletti, M., Webb, S. C. & Nascetti, G. Pseudoterranova and Contracaecum. In Molecular detection of human parasitic pathogens (ed. Liu, D.) 645–656 (CRC Press, 2012).
Google Scholar 
Mozgovoĭ, A.A., (1953) Ascaridata of animals and man, and the diseases caused by them. In: Osnovy nematodologii. Vol. II. Izd. AN SSSR, Moskva (In Russian)Johnston, T.H., Mawson, P.M., (1939) Internal parasites of the pigmy sperm whale. Rec. South Aust Museum.6. http://www.biodiversitylibrary.org/item/126147.Gibson, D. I. The systematics of ascaridoid nematodes-a current assessment. In Stone A (eds Platt, H. & Khalil, L.) 321–338 (Academic Press, 1983).
Google Scholar 
Shamsi, S., Barton, D. P. & Zhu, X. Description and characterisation of Terranova pectinolabiata n. sp. (Nematoda: Anisakidae) in great hammerhead shark, Sphyrna mokarran (Rüppell, 1837), in Australia. Parasitol. Res. 118, 2159–2168. https://doi.org/10.1007/s00436-019-06360-4 (2019).Article 
PubMed 

Google Scholar 
Shamsi, S., Barton, D. P. & Zhu, X. Description and genetic characterisation of Pulchrascaris australis n. sp. in the scalloped hammerhead shark, Sphyrna lewini (Griffin & Smith) in Australian waters. Parasitol. Res. https://doi.org/10.1007/s00436-020-06672-w (2020).Article 
PubMed 

Google Scholar 
González-Solís, D. et al. Parasitic nematodes of marine fishes from Palmyra Atoll, East Indo-Pacific, including a new species of Spinitectus (Nematoda, Cystidicolidae). Zookeys. 2019, 1–26 (2019).
Google Scholar 
Jabbar, A. et al. Larval anisakid nematodes in teleost fishes from Lizard Island, northern great barrier reef Australia. Mar. Freshw. Res. 63, 1283. https://doi.org/10.1071/MF12211 (2012).Article 

Google Scholar 
ICES. (2012) Pseudoterranova larvae (“codworm”; Nematoda) in fish. Revised and updated by Matt Longshaw. ICES Identification Leaflets for diseases and parasites of fish and shellfish. Leaflet No. 7. 4 pp.Arai, H. P. & Smith, J. W. Guide to the parasites of fishes of Canada part V: Nematoda. Zootaxa 4185, 1. https://doi.org/10.11646/zootaxa.4185.1.1 (2016).Article 

Google Scholar 
Hurst, H. J. Identification and description of larval Anisakis simplex and Pseudoterranova decipiens (anisakidae: Nematoda) from New Zealand waters. New Zeal J. Mar. Freshw. Res. 18, 177–186 (1984).
Google Scholar 
Hernández-Orts, J. S. et al. Description, microhabitat selection and infection patterns of sealworm larvae (Pseudoterranova decipiens species complex, Nematoda: Ascaridoidea) in fishes from Patagonia Argentina. Parasite Vector. 6, 1–15 (2013).
Google Scholar 
Shiraki, T. Larval nematodes of family anisakidae (Nematoda) in the northern sea of Japan as a causative agent of eosinophilic phlegmone of granuloma in the human gastro-intestinal tract. Acta Med. Biol. 22, 57–98 (1974).
Google Scholar 
Berland, B. Nematodes from some Norwegian marine fishes. Sarsia 2, 1–50. https://doi.org/10.1080/00364827.1961.10410245 (1961).Article 

Google Scholar 
George-Nascimento, M. & Llanos, A. Micro-evolutionary implications of allozymic and morphometric variations in sealworms Pseudoterranova sp. (Ascaridoidea: Anisakidae) among sympatric hosts from the Southeastern Pacific Ocean. Int. J. Parasitol. 25, 1163–1171 (1995).CAS 
PubMed 

Google Scholar 
Deardorff, T. L., Kliks, M. M., Rosenfeld, M. E., Rychlinski, R. A. & Desowitz, R. S. Larval, ascaridoid nematodes from fishes near the Hawaiian Islands, with commonents on pathogenicity experiments. Pacific Sci. 36, 187–201 (1982).
Google Scholar 
Deardorff, T. L., Kliks, M. M. & Desowitz, R. S. Histopathology induced by larval Terranova (Type HA) (nematoda: Anisakinae) in experimentally infected rats. J. Parasitol. 69, 191–195 (1983).CAS 
PubMed 

Google Scholar 
Kuramochi, T. et al. Stomach nematodes of the family anisakidae collected from the cetaceans stranded on or incidentally caught off the coasts of the Kanto districts and adjoining areas. Mem. Nat. Museum. Nat. Sci. 37, 177–192 (2001).
Google Scholar 
Deardorff, T. L., Raybourne, R. B. & Desowitz, R. S. Description of a third-stage larva, Terranova type Hawaii A (nematoda: Anisakinae), from Hawaiian fishes. J. Parasitol. 70, 829–831 (1984).CAS 
PubMed 

Google Scholar 
González-Solís, D., Vidal-Martínez, V. M., Antochiw-Alonso, D. M. & Ortega-Argueta, A. Anisakid nematodes from stranded pygmy sperm whales, Kogia breviceps (Kogiidae), in three localities of the Yucatan peninsula. Mexico. J. Parasitol. 92, 1120–1122 (2006).
Google Scholar 
Santos, C. P. & Lodi, L. Occurrence of Anisakis physeteris Baylis, 1923 and Pseudoterranova sp. (Nematoda) in pygmy sperm whale Kogia breviceps (De Blainvillei, 1838) (Physeteridae) in northeastern coast of Brazil. Mem. Inst. Oswaldo Cruz. 93, 187–188. https://doi.org/10.1590/s0074-02761998000200009 (1998).Article 
CAS 
PubMed 

Google Scholar 
Bloodworth, B. E. & Odell, D. K. Kogia breviceps (cetacea: Kogiidae). Mam. Species. 819, 1–12. https://doi.org/10.1644/819.1 (2008).Article 

Google Scholar 
Deardorff, T. L. & Overstreet, R. M. Terranova ceticola n. sp. (Nematoda: Anisakidae) from the dwarf sperm whale; Kogia simus (Owen), in the Gulf of Mexico. Syst. Parasitol. 3, 25–28 (1981).
Google Scholar 
Abollo, E., Santiago, P., (2002) SEM study of Anisakis brevispiculata Dollfus, 1966 and Pseudoterranova ceticola (Deardoff and Overstreet, 1981) (Nematoda: Anisakidae), parasites of the pygmy sperm whale Kogia breviceps. Sci. Mar. 66 3 49 255Di Deco, M. A., Orecchia, P., Paggi, L. & Petrarca, V. Morphometric stepwise discriminant analysis of three genetically identified species within Pseudoterranova decipiens (Krabbe, 1878) (Nematoda: Ascaridida). Syst. Parasitol. 29, 81–88 (1994).
Google Scholar 
George-Nascimento, M. & Urrutia, X. Pseudoterranova cattani sp. nov. (Ascaridoidea: Anisakidae), a parasite of the South American sea lion Otaria byronia De Blainville from Chile. Rev. Chil. Hist. Nat. 73, 93–98. https://doi.org/10.4067/s0716-078×2000000100010 (2000).Article 

Google Scholar 
Mattiucci, S. et al. Allozyme and morphological identification of Anisakis, Contracaecum and Pseudoterranova from Japanese waters (Nematoda, Ascaridoidea). Syst Parasitol. 40, 81–92 (1998).
Google Scholar 
Paggi, L. et al. Pseudoterranova decipiens species A and B (Nematoda, Ascaridoidea): Nomenclatural designation, morphological diagnostic characters and genetic markers. Syst. Parasitol. 45, 185–197. https://doi.org/10.1023/A:1006296316222 (2000).Article 
CAS 
PubMed 

Google Scholar 
Valentini, A. et al. Genetic relationships among Anisakis species (Nematoda: Anisakidae) inferred from mitochondrial cox2 sequences, and comparison with allozyme data. J. Parasitol. 92, 156–166 (2006).CAS 
PubMed 

Google Scholar 
Colón-Llavina, M. M. et al. Additional records of metazoan parasites from Caribbean marine mammals, including genetically identified anisakid nematodes. Parasitol Res. 105, 1239–1252 (2009).PubMed 

Google Scholar 
Cavallero, S., Nadler, S. A., Paggi, L., Barros, N. B. & D’Amelio, S. Molecular characterization and phylogeny of anisakid nematodes from cetaceans from southeastern Atlantic coasts of USA, Gulf of Mexico, and Caribbean Sea. Parasitol. Res. 108, 781–792 (2011).PubMed 

Google Scholar 
Kijewska, A., Dzido, J., Shukhgalter, O. & Rokicki, J. Anisakid parasites of fishes caught on the African shelf. J. Parasitol. 95, 639–645 (2009).PubMed 

Google Scholar 
Quiazon, K. M. A., Santos, M. D. & Yoshinaga, T. Anisakis species (nematoda: Anisakidae) of dwarf sperm whale kogia sima (Owen, 1866) stranded off the pacific coast of southern Philippine archipelago. Vet. Parasitol. https://doi.org/10.1016/J.VETPAR.2013.05.019 (2013).Article 
PubMed 

Google Scholar 
Zhang, L., Du, X., An, R., Li, L. & Gasser, R. B. Identification and genetic characterization of Anisakis larvae from marine fishes in the South China Sea using an electrophoretic-guided approach. Electrophoresis 34, 888–894 (2013).CAS 
PubMed 

Google Scholar 
Luo, H.-Y., Chen, H.-Y., Chen, H.-G. & Shih, H.-H. Scavenging hagfish as a transport host of anisakid nematodes. Vet. Parasitol. 218, 15–21. https://doi.org/10.1016/j.vetpar.2016.01.005 (2016).Article 
PubMed 

Google Scholar 
Kuhn, T., Hailer, F., Palm, H. W. & Klimpel, S. Global assessment of molecularly identified Anisakis dujardin, 1845 (nematoda: Anisakidae) in their teleost intermediate hosts. Folia Parasitol. (Praha). 60, 123–134. https://doi.org/10.14411/fp.2013.013 (2013).Article 
CAS 
PubMed 

Google Scholar 
Grainger, J. N. R. The Identity of the larval nematodes found in the body muscles of the cod (Gadus callarias L.). Parasitology 49, 121–131 (1959).CAS 
PubMed 

Google Scholar 
Costa, G., Chada, T., Melo-Moreira, E., Cavallero, S. & D’Amelio, S. Endohelminth parasites of the leafscale gulper shark, Centrophorus squamosus (Bonnaterre, 1788) (Squaliformes: Centrophoridae) off Madeira archipelago. Acta Parasitol. 59, 316–322. https://doi.org/10.2478/s11686-014-0247-x (2014).Article 
PubMed 

Google Scholar 
Hermida, M. et al. Infection levels and diversity of anisakid nematodes in blackspot seabream, Pagellus bogaraveo, from Portuguese waters. Parasitol. Res. 110, 1919–1928 (2012).PubMed 

Google Scholar 
Sequeira, V. et al. Macroparasites as biological tags for stock identification of the bluemouth, Helicolenus dactylopterus (Delaroche, 1809) in Portuguese waters. Fish Res. 106, 321–328. https://doi.org/10.1016/j.fishres.2010.08.014 (2010).Article 

Google Scholar 
Shamsi, S., Spröhnle-Barrera, C. & Shafaet, H. M. Occurrence of Anisakis spp. (Nematoda: Anisakidae) in a pygmy sperm whale Kogia breviceps (Cetacea: Kogiidae) in Australian waters. Dis. Aquat. Organ. 134, 65–74. https://doi.org/10.3354/dao03360 (2019).Article 
CAS 
PubMed 

Google Scholar 
Mcalpine, D. F., Murison, L. D. & Hoberg, E. P. New records for the pygmy sperm whale, Kogia breviceps (physeteridae) from Atlantic Canada with notes on diet and parasites. Mar. Mammal. Sci. 13, 701–704. https://doi.org/10.1111/j.1748-7692.1997.tb00093.x (1997).Article 

Google Scholar 
Gunter, G. & Overstreet, R. Cetacean notes. I. Sei and rorqual whales on the Mississippi coast, a correction. II. A dwarf sperm whale in Mississippi sound and its helminth parasites. Gulf Res. Rep. 4, 479–481 (1974).
Google Scholar 
Mignucci-Giannoni, A. A., Hoberg, E. P., Siegel-Causey, D. & Williams, E. H. Metazoan parasites and other symbionts of cetaceans in the Caribbean. J. Parasitol. 84, 939–946 (1998).CAS 
PubMed 

Google Scholar 
Vidal, O., Findley, L. T., Turk, P. J. & Boyer, R. E. Recent records of pygmy sperm whales in the Gulf of California. Mexico. Mar. Mammal. Sci. 3, 354–356. https://doi.org/10.1111/J.1748-7692.1987.TB00323.X (1987).Article 

Google Scholar 
Dollfus, R. P. Helminthofaune de Kogia breviceps (Blainxille, 1938) cetace odontocete. Recoltes du Dr R. Duguy. Ann. Sci. Natl. Charente-Maritime 4, 3–6 (1966).
Google Scholar 
MCAlpine, D.F., (2018) Pygmy and dwarf sperm whales. In: Encyclopedia of Marine Mammals. Elsevier p. 786–8.Fernández, R., Santos, M. B., Carrillo, M., Tejedor, M. & Pierce, G. J. Stomach contents of cetaceans stranded in the canary Islands 1996–2006. J. Mar. Biol. Assoc. United Kingdom. 89, 873–883 (2009).

Google Scholar 
Berrow, S., López Suárez, P., Jann, B., Ryan, C., Varela, J., Hazevoet, C.J., (2015) Recent and noteworthy records of Cetacea from the Cape Verde Islands. www.scvz.org. Accessed 1 Mar 2021.Mattiucci, S., Nascetti, G., (2008) Chapter 2 advances and trends in the molecular systematics of anisakid nematodes, with implications for their evolutionary ecology and host-parasite co-evolutionary processes. Adv. Parasitol. 66 47 148Measures, L.N., (2014) Anisakiosis and pseudoterranovosis. Reston, Virginia; https://doi.org/10.3133/cir1393McClelland, G. The trouble with sealworms (Pseudoterranova decipiens species complex, nematoda): A review. Parasitology 2002(124 Suppl), S183-203 (2009).
Google Scholar 
Alt, K. G., Cunze, S., Kochmann, J. & Klimpel, S. Parasites of three closely related Antarctic fish species (teleostei: Nototheniinae) from Elephant Island. Acta Parasitol. https://doi.org/10.1007/s11686-021-00455-8 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
McClelland, G. Phocanema decipiens (Nematoda: Anisakinae): Experimental infections in marine copepods. Can. J. Zool. 60, 502–509. https://doi.org/10.1139/z82-075 (1982).Article 

Google Scholar 
Marcogliese, D. J. Review of experimental and natural invertebrate hosts of sealworm (Pseudoterranova decipiens) and its distribution and abundance in macroinvertebrates in eastern Canada. NAMMCO Sci. Publ. 3, 27–37 (2001).
Google Scholar 
West, K. L. et al. Diet of pygmy sperm whales (Kogia breviceps) in the Hawaiian Archipelago. Mar. Mammal. Sci. 25, 931–943. https://doi.org/10.1111/j.1748-7692.2009.00295.x (2009).Article 

Google Scholar 
Kleinertz, S., Damriyasa, I. M., Hagen, W., Theisen, S. & Palm, H. W. An environmental assessment of the parasite fauna of the reef-associated grouper Epinephelus areolatus from Indonesian waters. J. Helminthol. 88, 50–63 (2014).CAS 
PubMed 

Google Scholar 
Nadler, S. A. et al. Molecular phylogenetics and diagnosis of Anisakis, Pseudoterranova, and Contracaecum from northern pacific marine mammals. J. Parasitol. 91, 1413–1429 (2005).CAS 
PubMed 

Google Scholar 
Weitzel, T. et al. Human infections with Pseudoterranova cattani nematodes. Chile. Emerg. Infect. Dis. 21, 1874–1875 (2015).CAS 
PubMed 

Google Scholar 
Arizono, N., Miura, T., Yamada, M., Tegoshi, T. & Onishi, K. Human infection with Pseudoterranova azarasi roundworm. Emerg. Infect. Dis. 17, 555–556. https://doi.org/10.3201/eid1703.101350 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Kleinertz, S. et al. Gastrointestinal parasites of free-living Indo-Pacific bottlenose dolphins (Tursiops aduncus) in the Northern Red Sea. Egypt. Parasitol Res. 113, 1405–1415. https://doi.org/10.1007/s00436-014-3781-4 (2014).Article 
CAS 
PubMed 

Google Scholar 
Aco Alburqueque, R., Palomba, M., Santoro, M. & Mattiucci, S. Molecular identification of zoonotic parasites of the genus Anisakis (Nematoda: Anisakidae) from fish of the southeastern Pacific Ocean (off Peru coast). Pathogens. 9, 910. https://doi.org/10.3390/pathogens9110910 (2020).Article 
CAS 
PubMed Central 

Google Scholar 
Di Azevedo, M. I. N., Carvalho, V. L. & Iñiguez, A. M. Integrative taxonomy of anisakid nematodes in stranded cetaceans from Brazilian waters: An update on parasite’s hosts and geographical records. Parasitol. Res. 116, 3105–3116. https://doi.org/10.1007/s00436-017-5622-8 (2017).Article 
PubMed 

Google Scholar 
Quiazon, K. M. A., Santos, M. D., Blatchley, D. D., Aguila, R. D. & Yoshinaga, T. Molecular and morphological identifications of Anisakis dujardin, 1845 (Nematoda: Anisakidae) from a rare deraniyagala’s beaked whale (Mesoplodon hotaula deraniyagala, 1963) and blainville’s beaked whale (Mesoplodon densirostris blainville, 1817) stranded. Philipp. J. Sci. 150, 823–835 (2021).
Google Scholar 
Bao, M. et al. Air-dried stockfish of Northeast Arctic cod do not carry viable anisakid nematodes. Food Cont. 116, 107322. https://doi.org/10.1016/j.foodcont.2020.107322 (2020).Article 
CAS 

Google Scholar 
Liu, G. H. et al. Mitochondrial phylogenomics yields strongly supported hypotheses for ascaridomorph nematodes. Sci. Rep. https://doi.org/10.1038/srep39248 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Hrabar, J. et al. Phylogeny and pathology of anisakids parasitizing stranded California sea lions (Zalophus californianus) in Southern California. Front Mar. Sci. https://doi.org/10.3389/fmars.2021.636626 (2021).Article 

Google Scholar  More

United Nations, Department of Economic and Social Affairs & Population Division. World Urbanization Prospects: The 2018 Revision (ST/ESA/SER.A/420). (United Nations, 2019).Wei, Y. D. & Ewing, R. Urban expansion, sprawl and inequality. Landsc. Urban Plan. 177, 259–265. https://doi.org/10.1016/j.landurbplan.2018.05.021 (2018).Article 

Google Scholar 
Ayers, A. C. & Rehan, S. M. Supporting bees in cities: How bees are influenced by local and landscape features. Insects 12 (2021).Grimm, N. B. et al. Global change and the ecology of cities. Science 319, 756–760. https://doi.org/10.1126/science.1150195 (2008).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Fahrig, L. Effects of habitat fragmentation on biodiversity. Annu. Rev. Ecol. Evol. Syst. 34, 487–515. https://doi.org/10.1146/annurev.ecolsys.34.011802.132419 (2003).Article 

Google Scholar 
Shochat, E. et al. Invasion, competition, and biodiversity loss in urban ecosystems. Bioscience 60, 199–208. https://doi.org/10.1525/bio.2010.60.3.6 (2010).Article 

Google Scholar 
Sánchez-Bayo, F. & Wyckhuys, K. A. G. Worldwide decline of the entomofauna: A review of its drivers. Biol. Cons. 232, 8–27. https://doi.org/10.1016/j.biocon.2019.01.020 (2019).Article 

Google Scholar 
Seibold, S. et al. Arthropod decline in grasslands and forests is associated with landscape-level drivers. Nature 574, 671–674. https://doi.org/10.1038/s41586-019-1684-3 (2019).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Wagner, D. L. Insect declines in the anthropocene. Annu. Rev. Entomol. 65, 457–480. https://doi.org/10.1146/annurev-ento-011019-025151 (2020).Article 
CAS 
PubMed 

Google Scholar 
Brown, M. J. & Paxton, R. J. The conservation of bees: A global perspective. Apidologie 40, 410–416 (2009).
Google Scholar 
Hallmann, C. A. et al. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS ONE 12, e0185809 (2017).PubMed 
PubMed Central 

Google Scholar 
Kennedy, C. M. et al. A global quantitative synthesis of local and landscape effects on wild bee pollinators in agroecosystems. Ecol. Lett. 16, 584–599 (2013).PubMed 

Google Scholar 
Potts, S. G. et al. Summary for policymakers of the assessment report of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services on pollinators, pollination and food production. (2016).Winfree, R., Aguilar, R., Vázquez, D. P., LeBuhn, G. & Aizen, M. A. A meta-analysis of bees’ responses to anthropogenic disturbance. Ecology 90, 2068–2076 (2009).PubMed 

Google Scholar 
Millard, J. et al. Global effects of land-use intensity on local pollinator biodiversity. Nat. Commun. 12, 1–11 (2021).ADS 

Google Scholar 
Baldock, K. C. et al. A systems approach reveals urban pollinator hotspots and conservation opportunities. Nat. Ecol. Evolut. 3, 363–373 (2019).
Google Scholar 
Banaszak-Cibicka, W., Twerd, L., Fliszkiewicz, M., Giejdasz, K. & Langowska, A. City parks vs. natural areas—Is it possible to preserve a natural level of bee richness and abundance in a city park?. Urban Ecosyst. 21, 599–613 (2018).
Google Scholar 
Hall, D. M. et al. The city as a refuge for insect pollinators. Conserv. Biol. 31, 24–29 (2017).PubMed 

Google Scholar 
Theodorou, P. et al. Urban areas as hotspots for bees and pollination but not a panacea for all insects. Nat. Commun. 11, 1–13 (2020).
Google Scholar 
Wilson, C. J. & Jamieson, M. A. The effects of urbanization on bee communities depends on floral resource availability and bee functional traits. PLoS ONE 14, e0225852 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Samuelson, A. E., Schürch, R. & Leadbeater, E. Dancing bees evaluate central urban forage resources as superior to agricultural land. J. Appl. Ecol. 59, 79–88 (2022).
Google Scholar 
Fortel, L. et al. Decreasing abundance, increasing diversity and changing structure of the wild bee community (Hymenoptera: Anthophila) along an urbanization gradient. PLoS ONE 9, e104679 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Roffet-Salque, M. et al. Widespread exploitation of the honeybee by early Neolithic farmers. Nature 527, 226–230 (2015).ADS 
CAS 
PubMed 

Google Scholar 
Crane, E. Recent research on the world history of beekeeping. Bee World 80, 174–186 (1999).
Google Scholar 
Dietemann, V., Pirk, C. W. W. & Crewe, R. Is there a need for conservation of honeybees in Africa?. Apidologie 40, 285–295 (2009).
Google Scholar 
Jaffe, R. et al. Estimating the density of honeybee colonies across their natural range to fill the gap in pollinator decline censuses. Conserv. Biol. 24, 583–593 (2010).PubMed 

Google Scholar 
Browne, K. A. et al. Investigation of free-living honey bee colonies in Ireland. J. Apic. Res. 60, 229–240. https://doi.org/10.1080/00218839.2020.1837530 (2021).Article 

Google Scholar 
Kohl, P. L. & Rutschmann, B. The neglected bee trees: European beech forests as a home for feral honey bee colonies. PeerJ 6, e4602 (2018).PubMed 
PubMed Central 

Google Scholar 
Oleksa, A., Gawroński, R. & Tofilski, A. Rural avenues as a refuge for feral honey bee population. J. Insect Conserv. 17, 465–472 (2013).
Google Scholar 
Rutschmann, B., Kohl, P. L., Machado, A. & Steffan-Dewenter, I. Semi-natural habitats promote winter survival of wild-living honeybees in an agricultural landscape. Biol. Cons. 266, 109450 (2022).
Google Scholar 
Thompson, C. E., Biesmeijer, J. C., Allnutt, T. R., Pietravalle, S. & Budge, G. E. Parasite pressures on feral honey bees (Apis mellifera sp). PLoS ONE 9, e105164 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Bila Dubaić, J. et al. Unprecedented density and persistence of feral honey bees in urban environments of a large SE-European City (Belgrade, Serbia). Insects 12, 1127 (2021).PubMed 
PubMed Central 

Google Scholar 
Alaux, C., Le Conte, Y. & Decourtye, A. Pitting wild bees against managed honey bees in their native range, a losing strategy for the conservation of honey bee biodiversity. Front. Ecol. Evol. 7, 60 (2019).
Google Scholar 
Requier, F. et al. The conservation of native honey bees is crucial. Trends Ecol. Evol. 34, 789–798 (2019).PubMed 

Google Scholar 
Mladenović, S. et al. Environment in Belgrade in 2018. (in Serbian: Kvalitet životne sredine u Beogradu u 2018. godini). (The City Administration, Secretariat for Environmental Protection, 2019).Statistical Office of the Republic of Serbia. https://data.stat.gov.rs/Home/Result/130202010207?languageCode=en-US.
(“The Official Gazette of the Republic of Serbia”, Nos. 41/2009, 93/2012 and 14/2106 [In Serbian], 2009).Johnson, M. T. & Munshi-South, J. Evolution of life in urban environments. Science 358, eaam8327 (2017).PubMed 

Google Scholar 
Jara, L. et al. Stable genetic diversity despite parasite and pathogen spread in honey bee colonies. Sci. Nat. 102, 1–8 (2015).
Google Scholar 
Tanasković, M. et al. MtDNA analysis indicates human-induced temporal changes of serbian honey bees diversity. Insects 12, 767 (2021).PubMed 
PubMed Central 

Google Scholar 
Wang, J. COANCESTRY: A program for simulating, estimating and analysing relatedness and inbreeding coefficients. Mol. Ecol. Resour. 11, 141–145. https://doi.org/10.1111/j.1755-0998.2010.02885.x (2011).Article 
PubMed 

Google Scholar 
Wang, J. Triadic IBD coefficients and applications to estimating pairwise relatedness. Genet. Res. 89, 135–153. https://doi.org/10.1017/s0016672307008798 (2007).Article 
CAS 
PubMed 

Google Scholar 
Jacobson, S. Locally adapted, varroa resistant honey bees: ideas from several key studies. Am. Bee J. (2010).McNeely, J. A., Miller, K. R., Reid, W. V., Mettermeier, R. A. & Werner, T. B. Conserving the world’s biological diversity. (UICN, Morges (Suiza) WRI, Washington DC (EUA) CI, Washington DC (EUA) WWF …, 1990).Hoban, S. M. et al. Bringing genetic diversity to the forefront of conservation policy and management. Conserv. Genet. Resour. 5, 593–598 (2013).
Google Scholar 
Hohenlohe, P. A., Funk, W. C. & Rajora, O. P. Population genomics for wildlife conservation and management. Mol. Ecol. 30, 62–82 (2021).PubMed 

Google Scholar 
Shafer, A. B. et al. Genomics and the challenging translation into conservation practice. Trends Ecol. Evol. 30, 78–87 (2015).PubMed 

Google Scholar 
Mattila, H. R. & Seeley, T. D. Genetic diversity in honey bee colonies enhances productivity and fitness. Science 317, 362–364 (2007).ADS 
CAS 
PubMed 

Google Scholar 
Oddie, M. A. & Dahle, B. Insights from Norway: Using natural adaptation to breed Varroa-resistant honey bees. Bee World 98, 38–43 (2021).
Google Scholar 
Oddie, M. A., Dahle, B. & Neumann, P. Norwegian honey bees surviving Varroa destructor mite infestations by means of natural selection. PeerJ 5, e3956 (2017).PubMed 
PubMed Central 

Google Scholar 
Oldroyd, B. P. & Fewell, J. H. Genetic diversity promotes homeostasis in insect colonies. Trends Ecol. Evol. 22, 408–413 (2007).PubMed 

Google Scholar 
Tarpy, D. R. Genetic diversity within honeybee colonies prevents severe infections and promotes colony growth. Proc. R Soc. Lond. Series B Biol. Sci. 270, 99–103 (2003).
Google Scholar 
van Baalen, M. & Beekman, M. The costs and benefits of genetic heterogeneity in resistance against parasites in social insects. Am. Nat. 167, 568–577 (2006).PubMed 

Google Scholar 
Eckholm, B. J., Anderson, K. E., Weiss, M. & DeGrandi-Hoffman, G. Intracolonial genetic diversity in honeybee (Apis mellifera) colonies increases pollen foraging efficiency. Behav. Ecol. Sociobiol. 65, 1037–1044 (2011).
Google Scholar 
Graham, S., Myerscough, M., Jones, J. & Oldroyd, B. Modelling the role of intracolonial genetic diversity on regulation of brood temperature in honey bee (Apis mellifera L.) colonies. Insectes Soc. 53, 226–232 (2006).
Google Scholar 
Jones, J. C., Myerscough, M. R., Graham, S. & Oldroyd, B. P. Honey bee nest thermoregulation: Diversity promotes stability. Science 305, 402–404 (2004).ADS 
CAS 
PubMed 

Google Scholar 
Tanasković, M. et al. Further evidence of population admixture in the Serbian honey bee population. Insects 13, 180 (2022).PubMed 
PubMed Central 

Google Scholar 
Nedić, N. et al. Detecting population admixture in honey bees of Serbia. J. Apic. Res. 53, 303–313. https://doi.org/10.3896/IBRA.1.53.2.12 (2014).Article 

Google Scholar 
Nedić, N., Stanisavljević, L., Mladenović, M. & Stanisavljević, J. Molecular characterization of the honeybee Apis mellifera carnica in Serbia. Arch. Biol. Sci. 61, 587–598 (2009).
Google Scholar 
Kükrer, M., Kence, M. & Kence, A. Honey bee diversity is swayed by migratory beekeeping and trade despite conservation practices: Genetic evidence for the impact of anthropogenic factors on population structure. Front. Ecol. Evolut. 9 (2021).Bouga, M., Harizanis, P. C., Kilias, G. & Alahiotis, S. Genetic divergence and phylogenetic relationships of honey bee Apis mellifera (Hymenoptera: Apidae) populations from Greece and Cyprus using PCR–RFLP analysis of three mtDNA segments. Apidologie 36, 335–344 (2005).CAS 

Google Scholar 
Dall’Olio, R., Marino, A., Lodesani, M. & Moritz, R. F. Genetic characterization of Italian honeybees, Apis mellifera ligustica, based on microsatellite DNA polymorphisms. Apidologie 38, 207–217 (2007).CAS 

Google Scholar 
Neumann, P. & Blacquière, T. The Darwin cure for apiculture? Natural selection and managed honeybee health. Evol. Appl. 10, 226–230 (2017).PubMed 

Google Scholar 
Kulinčević, J., Rinderer, T., Mladjan, V. & Buco, S. Five years of bi-directional genetic selection for honey bees resistant and susceptible to Varroa jacobsoni. Apidologie 23, 443–452 (1992).
Google Scholar 
2011 Census of Population, Households and Dwellings in the Republic of Serbia: Comparative Overview of the Number of Population in 1948, 1953, 1961, 1971, 1981, 1991, 2002 and 2011. (Statistical Office of the Republic of Serbia, 2014).Techer, M. A. et al. Large-scale mitochondrial DNA analysis of native honey bee Apis mellifera populations reveals a new African subgroup private to the South West Indian Ocean islands. BMC Genet. 18, 1–21 (2017).
Google Scholar 
Garnery, L., Cornuet, J. M. & Solignac, M. Evolutionary history of the honey bee Apis mellifera inferred from mitochondrial DNA analysis. Mol. Ecol. 1, 145–154 (1992).CAS 
PubMed 

Google Scholar 
Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35, 1547 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Excoffier, L. & Lischer, H. E. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resources. 10, 564–567 (2010).
Google Scholar 
Kalinowski, S. T. hp-rare 1.0: A computer program for performing rarefaction on measures of allelic richness. Mol. Ecol. Notes. 5, 187–189 (2005).CAS 

Google Scholar 
Stoneking, M., Hedgecock, D., Higuchi, R. G., Vigilant, L. & Erlich, H. A. Population variation of human mtDNA control region sequences detected by enzymatic amplification and sequence-specific oligonucleotide probes. Am. J. Hum. Genet. 48, 370 (1991).CAS 
PubMed 
PubMed Central 

Google Scholar 
Hammer, Ø., Harper, D. A. & Ryan, P. D. PAST: Paleontological statistics software package for education and data analysis. Palaeontol. Electron. 4, 9 (2001).
Google Scholar 
Crozier, R. & Crozier, Y. The mitochondrial genome of the honeybee Apis mellifera: Complete sequence and genome organization. Genetics 133, 97–117 (1993).CAS 
PubMed 
PubMed Central 

Google Scholar 
Falush, D., Stephens, M. & Pritchard, J. K. Inference of population structure using multilocus genotype data: Linked loci and correlated allele frequencies. Genetics 164, 1567–1587. https://doi.org/10.1093/genetics/164.4.1567 (2003).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Falush, D., Stephens, M. & Pritchard, J. K. Inference of population structure using multilocus genotype data: Dominant markers and null alleles. Mol. Ecol. Notes 7, 574–578 (2007).CAS 
PubMed 
PubMed Central 

Google Scholar 
Hubisz, M. J., Falush, D., Stephens, M. & Pritchard, J. K. Inferring weak population structure with the assistance of sample group information. Mol. Ecol. Resour. 9, 1322–1332. https://doi.org/10.1111/j.1755-0998.2009.02591.x (2009).Article 
PubMed 
PubMed Central 

Google Scholar 
Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 155, 945–959 (2000).CAS 
PubMed 
PubMed Central 

Google Scholar 
Earl, D. A. & VonHoldt, B. M. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361 (2012).
Google Scholar 
Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 14, 2611–2620 (2005).CAS 
PubMed 

Google Scholar 
Jombart, T., Devillard, S. & Balloux, F. Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. BMC Genet. 11, 94. https://doi.org/10.1186/1471-2156-11-94 (2010).Article 
PubMed 
PubMed Central 

Google Scholar 
Ward, J. H. Hierarchical grouping to optimize an objective function. J. Am. Stat. Assoc. 58, 236–244. https://doi.org/10.1080/01621459.1963.10500845 (1963).Article 
MathSciNet 

Google Scholar 
Wang, J. An estimator for pairwise relatedness using molecular markers. Genetics 160, 1203–1215. https://doi.org/10.1093/genetics/160.3.1203 (2002).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Li, C. C., Weeks, D. E. & Chakravarti, A. Similarity of DNA fingerprints due to chance and relatedness. Hum. Hered. 43, 45–52. https://doi.org/10.1159/000154113 (1993).Article 
CAS 
PubMed 

Google Scholar 
Lynch, M. Estimation of relatedness by DNA fingerprinting. Mol. Biol. Evol. 5, 584–599. https://doi.org/10.1093/oxfordjournals.molbev.a040518 (1988).Article 
CAS 
PubMed 

Google Scholar 
Lynch, M. & Ritland, K. Estimation of pairwise relatedness with molecular markers. Genetics 152, 1753–1766. https://doi.org/10.1093/genetics/152.4.1753 (1999).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Ritland, K. Estimators for pairwise relatedness and individual inbreeding coefficients. Genet. Res. 67, 175–185 (1996).
Google Scholar 
Queller, D. C. & Goodnight, K. F. Estimating relatedness using genetic markers. Evolution 43, 258–275. https://doi.org/10.1111/j.1558-5646.1989.tb04226.x (1989).Article 
PubMed 

Google Scholar 
Milligan, B. G. Maximum-likelihood estimation of relatedness. Genetics 163, 1153–1167 (2003).PubMed 
PubMed Central 

Google Scholar 
del Felipe, P. et al. Genetic diversity and structure of the commercially important native fish pacu (Piaractus mesopotamicus) from cultured and wild fish populations: Relevance for broodstock management. Aquacult. Int. 29, 289–305. https://doi.org/10.1007/s10499-020-00626-w (2021).Article 

Google Scholar  More

Wolf, S. et al. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-022-01904-x (2022).Article 

Google Scholar 
Kattge, J. et al. Glob. Change Biol. 26, 119–188 (2020).Article 

Google Scholar 
Sabatini, F. M. et al. Glob. Ecol. Biogeogr. 30, 1740–1764 (2021).Article 

Google Scholar 
Łopucki, R., Kiersztyn, A., Pitucha, G. & Kitowski, I. Ecol. Modell. 468, 109964 (2022).Article 

Google Scholar 
Kosmala, M., Wiggins, A., Swanson, A. & Simmons, B. Front. Ecol. Environ. 14, 551–560 (2016).Article 

Google Scholar 
White, C. R. et al. Funct. Ecol. 35, 1572–1578 (2021).Article 

Google Scholar 
Xirocostas, Z. A., Debono, S. A., Slavich, E. & Moles, A. T. Methods Ecol. Evol. 13, 596–602 (2022).Article 

Google Scholar 
Callaghan, C. T. et al. Bioscience 71, 55–63 (2020).
Google Scholar  More

Data collectionIn phase 1 (2017), we first circulated a comprehensive multi-language questionnaire and associated guidelines on protected area personnel numbers to major national protected area agencies, focusing on the 50 countries listed in the WDPA as having the most protected areas. The questionnaire requested information on personnel numbers, type of employers and management levels (from executive to skilled practical workers). Protected area personnel were defined as those spending at least 50% of their work time on protected area-related tasks. The questionnaire also requested information about job titles used for personnel equivalent to rangers. This phase produced usable data for 28 countries/territories.In phase 2 (2018 onwards), we conducted online searches for published data on protected area personnel numbers in the countries/territories not included in the questionnaire survey or where questionnaire responses were incomplete or unclear. The resulting information came from official organizational reports (10 countries/territories), published external studies, project documents and journal papers (35 countries/territories) and websites of protected area organizations or individual sites (9 countries/territories).In phase 3 (2018–2021), we directly requested personal contacts to locate or supply information from official sources both for the remaining countries/territories and to improve or verify data from phases 1 and 2. The minimum data requested were the overall number of protected area personnel, the number of those personnel that could be categorized as rangers, the terrestrial area of protected areas managed by the listed personnel and the source of the information. This phase contributed usable data for 68 countries and territories. Data for a further 17 countries/territories were assembled from multiple sources.The final dataset covered 176 countries/territories: 167 surveyed countries/territories and a further 9 countries/territories that have no WDPA-listed protected areas (Supplementary Table 1), with contributions from more than 150 individuals.Initial data processingTo assess and, where necessary, improve the reliability of data obtained in a wide range of formats and levels of detail and from multiple sources, we scored the data for each country/territory from 0 to 5 for each of four criteria—detail, accuracy, source and age of the data—with a maximum score of 20 (Supplementary Table 1 and Supplementary Fig. 1). For all low-scoring records (a score of less than 15), we sought more-reliable sources in later phases of the study, rejecting any final scores of less than 10.On reviewing the data, we excluded from the analysis protected areas identified in the WDPA as predominantly or entirely marine, Antarctica and countries/territories categorized in the WDPA as polar (Greenland, French Southern Territories, Bouvet Island, Heard Island and McDonald Islands, South Georgia and the South Sandwich Islands). These large, remote and/or largely uninhabited areas are likely to have quite different management models and scales of staffing from terrestrial protected areas (although marine protected areas are also widely understaffed11). For example, in 2012 the 972,000 km2 of Northeast Greenland Protected Area (categorized by the WDPA as polar) was only periodically visited by six two-person teams of naval personnel47, and the 2008 management plan of the 1.51 million km2 Papahānaumokuākea Marine National Monument (Hawai’i, USA) specifies just nine personnel, working in conjunction with several other agencies48. Data for one country were supplied by officials on the agreement that the country was not specifically identified in publications (the country is given the three-letter code ZZZ in relevant tables and figures).Because the format, completeness and level of detail of the data varied widely, from comprehensive personnel lists to single figures, we restricted our raw dataset to six variables that could be consistently extracted from data obtained for each country/territory:

1.

Total number of non-ranger personnel (if known)

2.

Total number of rangers (if known)

3.

Total number of protected area personnel (either the sum of 1 and 2 or provided as an undifferentiated total)

4.

Terrestrial area of protected areas covered by surveyed personnel (km2)

5.

Total terrestrial area of protected areas of the country/territory (km2)

6.

Year of the data

We used the WDPA, official publications and websites to determine (or verify) the area of terrestrial protected areas covered by the personnel listed for each country/territory, using WDPA data if there were discrepancies. Total national terrestrial protected area coverage was taken from the WDPA, with the exception of Turkey, where the area officially reported to the WDPA is significantly less than the nationally published area.The raw data from the survey are shown in Supplementary Table 1.Candidate predictorsTo predict the number of rangers and non-rangers in countries and territories for which we had no data (Statistical analysis), we collected information on the following set of variables, hereafter referred to as candidate predictors:Location dataThe WGS84 latitude and longitude of the centroid of the largest land mass associated with each country/ territory (to obtain the polygons defining the land masses, we used the R package rnaturalearth version 0.1.0; https://github.com/ropensci/rnaturalearth)2020 data from the World Bank (https://data.worldbank.org/indicator)

Area of the country/territory

Population density: the mid-year population divided by land area

Gross domestic product (GDP) in US dollars

GDP per capita in US dollars (GDP divided by mid-year population)

Growth rate of GDP

The proportion of rural inhabitants

The proportion of unemployed inhabitants

The forested proportion of the country/territory

2020 data for each country/territory from the WDPA (https://www.protectedplanet.net/)

The total terrestrial area of WDPA-listed protected areas

The proportion of the terrestrial area of all IUCN-categorized protected areas (Categories I–VI) that falls within protected areas in Category I or II

The proportion of the terrestrial area of all IUCN-categorized protected areas (Categories I–VI) that falls within protected areas in Categories I–IV

2020 data from the Yale Center for Environmental Law and Policy Environmental Performance Index (https://epi.yale.edu/)

Environmental Performance Index (EPI): a composite index using 32 performance indicators across 11 categories

Ecosystem Vitality Index (EVI): an indicator of how well countries preserve, protect and enhance ecosystems and the services they provide

Species Protection Index (SPI): an indicator of the species-level ecological representativeness of each country’s/territory’s protected area network

Not all this information was available for all countries/territories. Most of the missing data were for small territories that account for only a very small proportion of the total area of protected areas worldwide (Supplementary Table 2c).Statistical analysisOur primary objective was to estimate the total number of all personnel engaged in managing all the world’s WDPA-listed terrestrial protected areas and the number categorized as rangers. Our raw data collection yielded full, partial or no information on total personnel and ranger numbers for each country/territory (Supplementary Table 1 shows the completeness of all the data collected). Our first task, therefore, was (1) to impute the information for unsurveyed protected areas on the basis of information from surveyed protected areas within the same countries/territories and (2) to predict those numbers for countries/territories where no information was available on overall personnel numbers and/or ranger numbers on the basis of relationships we could establish between available information and candidate predictors in other countries/territories (Supplementary Table 7). A brief description of these two approaches follows, and full details on the analysis are provided in Supplementary Information.Data imputationFor countries/territories where we had obtained information about numbers of personnel and/or rangers for only some protected areas, our strategy was to populate the unsurveyed protected areas in proportion to the densities of personnel or rangers from the surveyed protected areas of the same countries/territories. For example, for Spain we obtained evidence that there are 619 rangers responsible for protected areas covering 44,328 km2, out of a national total protected area system covering 142,573 km2. To impute the number of rangers for the remaining 98,245 km2, we used the density of rangers in the surveyed area (one ranger per 44,328/619 = 71.6 km2) and applied that to the unsurveyed area, giving a total of 1,991 rangers (619 + (98,245/71.6)). This imputation assumes that unsurveyed areas are staffed at the same density as surveyed areas, whereas in reality the relative densities are likely to vary in unknown ways within different countries/territories. To study the sensitivity of our results to the assumed proportion, we repeated our analysis using the following proportions of the observed densities: 0, 0.25, 0.50, 0.75 and 1.00. This provided a range of personnel numbers from a minimum (based on a proportion of 0) to a presumed maximum (based on a proportion of 1.00). From the data obtained, it was not possible to calculate the actual proportions, but based on the experience of the practitioners in the author team, the unsurveyed areas are highly unlikely to be staffed at higher densities than surveyed areas and, on average, are very likely to be staffed at lower densities. After all, most survey respondents were national or subnational agencies responsible for protected areas subject to stronger formal requirements for protection and management and therefore likely to have larger workforces. Unsurveyed protected areas are more likely to be managed by local entities, with fewer resources, less-stringent management obligations and therefore fewer personnel. The range of proportions we considered to populate unsurveyed areas should therefore yield predictions encompassing the actual (unknown) numbers of rangers and non-rangers with a conservative margin of error. In the main text, we have reported the results of imputation assuming a proportion of 1, which is probably the most optimistic assessment of the current workforce in protected areas within the proportions of the observed densities considered. Results using lower proportions are shown in Extended Data Fig. 2 and Supplementary Tables 4 and 5.Data predictionOur imputation approach was not possible for countries/territories where (1) zero ranger or personnel data had been obtained and (2) specific data had not been obtained that allowed imputation either for rangers or for total personnel (where only total personnel numbers or only ranger numbers had been obtained). To predict the missing information, we used two different statistical approaches: linear mixed models (LMMs)49 and a general implementation of random forests, which we term RF/ETs because it encompasses both random forests sensu stricto (RFs)50 and a variant called extremely randomized trees (ETs)51. LMMs and RFs have been extensively discussed and reviewed in the literature49,52,53. We adopted these approaches because both have proved successful in producing accurate predictions for a wide range of applications and because both are well suited to our data since they both produce predictions from a set of predictors and allow for the consideration of spatial effects54,55. Furthermore, comparing predictions generated through very different methods informs us about the robustness of our results with respect to key statistical assumptions. LMMs come from the ‘data modelling culture’56 and belong to parametric statistics; RF/ETs come from the ‘algorithmic modelling culture’ and belong to non-parametric statistics.We followed the same workflow for both statistical approaches, comprising eight steps: (1) general data preparation; (2) preparation of initial training datasets; (3) selection of predictor variables and of the method used for handling spatial autocorrelation; (4) preparation of final training datasets; (5) fine tuning; (6) final training; (7) preparation of datasets for predictions and simulations; and (8) predictions and simulations (see Supplementary Information for details).Both approaches yielded very similar results with our data. We chose to present the LMM results in the main text, but we provide and compare the results obtained by both approaches in Supplementary Information.SoftwareWe performed all the data analyses using the free open-source statistical software R version 4.157. We used the R package spaMM version 3.9.13 to implement LMMs58 and the R package ranger version 0.13.1 to implement RF/ETs59. To reformat and plot the data, we used the Tidyverse suite of packages60. Details are provided in an R package we specifically developed so that findings presented in this paper can readily be reproduced (see Code availability). Using a workstation with an AMD Ryzen Threadripper 3990 × 64-core processor and 256 GB of RAM, our complete workflow ran in ~3,000 CPU hours.Estimation of required numbers and densities of personnelTo estimate the numbers of personnel and rangers required for effective management of existing protected areas, we referred to ref. 25. This estimates that the minimum budget needed to adequately manage the existing protected area system is US$67.6 billion per year and that current annual expenditure is US$24.3 billion. From these figures, we can calculate that resources invested in the current global system of protected areas are approximately 36% of what is required. We consulted data from https://ourworldindata.org to determine that the proportion of global public expenditure on employee compensation has remained between 21.01% and 23.33% in the years from 2006 to 2019. We obtained these figures from the ‘Government Spending’ section of the site, consulting the chart ‘Share of employee compensation in public spending, 2002 to 2019’ and selecting data for ‘World’. On the basis of this broadly constant proportion and the assumption that total employee compensation is an indicator of total employee numbers, we inferred that current numbers of protected area employees are also around 36% of what is required. We therefore multiplied our estimations of personnel and ranger numbers by 1/0.36 and recalculated the densities on this basis (current requirement = 1/0.36 × current estimate).To estimate staffing requirements for 30% global coverage of protected areas—the global target intended to be reached by 2030—we used the mean personnel and ranger densities calculated as being required at present to ‘populate’ a global area of terrestrial protected areas if increased from the percentage at the time of our study (15.7%) to 30% (current requirement × (0.300/0.157)).Economic calculationsWe based our calculations on published data from 202025, which estimate that expanding the protected areas to 30% would generate higher overall output (revenues) than non-expansion (an extra US$64–454 billion per year by 2050). This figure is only an indicative, partial estimate, generated for the purposes of comparison and to illustrate the substantial return on investment that protected area staff investments imply. Using these figures and our estimates of personnel requirements to ensure effective management of 30% coverage, we calculated the range of sums that each additional protected area staff member has the potential to generate (Supplementary Table 8). For clarity, we rounded these figures to the nearest hundred US dollars in the main text.Our estimates of the gross value added per worker in forestry and agriculture (sectors responsible for similar proportions of the world as protected areas) are included to provide a point of comparison for the figures showing the economic benefit generated per protected area personnel member (see the preceding). The data for the gross annual value of world agricultural production (US$3,550,231,736,000) and the number of workers employed in agriculture (343,527,711) come from the Food and Agriculture Organization of the United Nations30, providing an average gross value of annual agricultural production per worker of US$10,335. We adjusted these 2018 data to 2020 price levels using a deflator based on the US consumer price index (CPI) from the World Economic Outlook database61 (Supplementary Table 9). This ensures that all the economic value data we present are directly comparable for protected area, agricultural and forestry workers. We calculated the gross value of forest production per worker on the basis of direct contribution of forestry of more than US$539 billion to world GDP in 201162 and total forest-sector employment of 11.881 million full-time-equivalent jobs in 201032. These were the most up-to-date global estimates we could locate from credible sources that presented comparable estimates of forest-sector employment and contribution to GDP. This gives an average gross value of forest production per worker of US$45,367 per year. We used the same method as for agriculture to bring these figures to 2020 price levels (Supplementary Table 9). These figures are rounded to the nearest hundred US dollars in the main text. More

Study areaThe present study was conducted in Uttarkashi district, Uttarakhand, located between 38° 28′ to 31°28′ N latitude and 77°49′ to 79°25′ E longitude with an area of about 8016 km2, covering primarily hilly terrain with an altitudinal range of 715–6717 m (Fig. 3). The terrain is mountainous, consisting of undulating hill ranges and narrow valleys with temperate climatic conditions. The district lies in the upper catchment of two major rivers of India, viz., the Ganges (Bhagirathi towards upstream) and the Yamuna. The major vegetation types of the study area are Himalayan moist temperate forest, sub-alpine forest and alpine scrub59. The Uttarkashi district forests are managed under three Forest Divisions viz., (i) Uttarkashi Forest Division (ii) Upper Yamuna Badkot Forest Division and (iii) Tons Forest Division) with two Protected Areas (PAs) (i) Gangotri National Park and (ii) Govind Pashu Vihar National Park. The forested habitats of the study landscape are home to top conservation priority species, including Asiatic Black bear (Ursus thibetanus), Musk deer (Moschus spp.), Common leopard (Panthera pardus), Himalayan brown bear (Ursus arctos isabellinus) and Western Tragopan (Tragopan melanocephalus), Himalayan monal (Lophophorus impejanus). The study was conducted after a study permit issued by the Chief Wildlife Warden, Forest Department, Uttarakhand government, vide letter no. 848/5-6 dated 31/08/2019, we have not handled the species for doing research. Instead, remote camera traps have been used for collecting the data with the permission of the Chief Wildlife Warden, Government of Uttarakhand. Further, informed consent was taken before interviewing the local communities. The data was collected according to the institutional guidelines and approved by the Research Advisory and Monitoring Committee of the Zoological Survey of India.Figure 3Map of the study area Uttarkashi, Uttarakhand. ArcGIS 10.6 (ESRI, Redlands, CA) was used to create the map. (Map created using ArcGIS 10.6; http://www.esri.com).Full size imageSampling protocolThe basic sampling protocol and assumptions for multi-species occupancy modelling are identical to the single-species case7. Briefly, a set of 62 intensive sites, were randomly selected, and each site i was surveyed j times. During each survey, detection/non-detection of S focal species was recorded. Additionally, direct or indirect evidences of species presence from the different areas were also recorded.Data collectionThe complete study area was divided into 10 × 10 km grids, consisting of n = 60 grids. Based on the reconnaissance survey, out of these 60 grids, we selected 25 girds that were accessible to conduct the survey and have the species presence. Further, these grids were divided into 2 × 2 km grids to maximize our effort so that all logistically accessible grids could be covered, and we conducted intensive sampling in N = 62 grids after excluding the grids with human settlements. T The field surveys were conducted during 2018–2019, and a team of researchers systematically visited selected grids to collect data on the detection/non-detection of these ungulates. A total of 62 camera traps were deployed in selected grids, and 650 km were traversed, accounting for N = 54 trails in these sampled grids. These camera traps were visited once in every fifteen days for replacing the batteries as well as documenting the presence of the species through the sign surveys. The ultra-compact SPYPOINT FORCE-11D trail camera (SPYPOINT, GG Telecom, Canada, QC) and Browning trail camera (Defender 850, 20 MP, Prometheus Group, LLC Birmingham, Alabama, https://browningtrailcameras.com) camera traps were used to detect the presence/absence of ungulate species. The cameras were mounted 40–60 cm above ground on natural trails without lures.Data explorationWhile deploying camera traps, we also noted habitat variables through on-site observation such as distance to the village and human disturbance. We tested site covariates for collinearity and discarded one of a pair if the Pearson’s correlation was greater than 0.760. Hence, we assumed each of the site covariates could influence the occupancy and detectability of these ungulates.CovariatesWe hypothesized that habitat variables may influence these ungulates’ occupancy and detection probability. A total of 21 variables were extracted either from the field or using the ArcGIS v. 10.6 software (ESRI, Redlands, CA), and only 14 were retained after collinearity testing60 (Table 3). These covariates were classified into the following categories (Topographic variables, Habitat variables and anthropogenic variables). The topographic variables (elevation, slope and aspect) were generated using 30× resolution SRTM (Shuttle Radar Topography Mission) image downloaded from EarthExplorer (https://earthexplorer.usgs.gov/). The habitat/ land cover classification was carried out using Landsat 8 satellite imagery (Spatial resolution = 30 m) downloaded from Global Land Cover Facility by following the methodology suggested by61 using the ArcGIS v. 10.6 software (ESRI, Redlands, CA). The study area was classified into nine Land use/land cover (LULC) classes viz., West Himalayan Sub-alpine birch/fir Forest (FT 188), West Himalayan upper oak/fir forest (FT 162), West Himalayan Dry juniper forest (FT 180), Ban oak forest (FT 152), Moist Deodar Forest (FT 155), Western mixed coniferous forest (FT 156), Moist temperate Deciduous Forest (FT 157) which were used for further analysis considering their importance to species ecology and behavior60. The values for all the covariates were extracted at 30 m resolution, and a single value per site was obtained by averaging all the pixel values within each sampling site (camera trap locations).Table 3 Habitat variables used for multi species occupancy analysis of three ungulate species in Uttarkashi, Uttarakhand.Full size tableOccupancy modelling frameworkWe used multi-species occupancy modelling62 of barking deer, goral and sambar to estimate the probability of the species (s) occurred within the area (i) sampled during our survey period (j), for accounting the imperfect detection of the species8. Distinguishing the true presence/absence of a species from detection/non-detection (i.e., species present and captured or species present but not captured) requires spatially or temporally replicated data. We used camera stations to record the presence/absence of species along with sign survey in all the studied grids. The camera traps were placed along the trail/transects in the studied grids hence each grid needs to be visited once in every fifteen days to check the camera traps as well as to document the presence of the studied species. Therefore, we treated 15 trap nights as one sampling occasion at a particular camera station resulting in ~ 7 sampling occasions per camera station.Our aim was to record the presence/ absence of the species at a particular gird hence we incorporated sign survey data if the species was not detected in camera station but recorded through sign survey. We pooled the presence/absence data in a single sheet of each species following6 and fitted occupancy and detectability models using programme Mark63,64. We model the species (s) presence (ysij = 1) and absence (ysij = 0) at site i during survey j, and the sampling protocol was identical to single species case65, where the Bernoulli random variable was conditional on the presence of species s (Zs = 1) following6$${text{y}}_{sij} sim {text{ Bernoulli}}left( {{text{p}}_{sij} {text{z}}_{si} } right),$$
where Psij represents the probability of detecting species S during replicate survey j at site i and Zsi = presence or absence of species s at site i.Furthermore, we model the latent occupancy state of species s at site i as a multivariate Bernoulli random variable:$${text{Z}}_{i} sim {text{MVB}}left( {uppsi _{i} } right)$$
where Zi = {Z1i, Z2i….., ZSi} is an S-dimensional vector of 1’s and 0’s denoting the latent occupancy state of all S species and (ψi) is a 2S-dimensional vector denoting the probability of all possible sequences of 1’s and 0’s Zi can attain such that ∑ ψi = 1 with corresponding probability mass function (PMF) adopted from6,64.$$fleft( {{text{Z}}_{i} } right) = {text{ exp}}left( {left( {{text{Z}}_{i} {text{log}}(uppsi_{{text{i}}} {1}/uppsi_{{text{i}}} 0} right) , + {text{ log}}left( {uppsi_{{text{i}}} 0} right)} right).$$The quantity f = log (ψi1/ψi0), is the log odds species S occupies a site often referred to as a ‘natural parameter’.Since we are modeling three ungulate species (S = 3), 2S = 23 the possible encounter histories included in the dataset were eight, if neither of the two species were detected the value of ‘00’ was assigned; similarly ‘01’ indicates detection of species 1; ‘02’ indicates detection of species 2; ‘03’ indicates detection of both the species; ‘04’ indicates detection of species 3; ‘05’ indicates detection of species 1 and species 3; ‘06’ indicates detection of species 2 and species 3 and ‘07’ indicates detection of all the three species. We modelled constant occupancy and detection probability for each of the three species. Hence, we specified 6 f and p parameters, an intercept (β) for each of one-way f parameter and detection parameter p following64.$$f_{{1}}=upbeta_{{{1},}} ;;{text{p}}=upbeta_{{4}}$$$$f_{{2}} = upbeta_{{{2},}} ;{text{p }} = , upbeta 5$$$$f_{{3}} = , upbeta_{{{3},}}; {text{p }} = , upbeta_{{6}}$$We fit a set of models including the detection probability as a constant, p(.), and variable function to occupancy ψ(covariate) for site-specific covariates and models include occupancy as constant ψ(.) and variable function of the detection p(covariates) for the respective site covariates.As we have assumed the independence among all three species, the model shows marginal occupancy probabilities of species 1, species 2 and species 3 varies as a function of environmental variables. We incorporated site-level characteristics affecting species-specific occurrence (f1: occupancy of species 1, f2: occupancy of species 2, & f3: occupancy of species 3) and detection probabilities using a generalized linear modelling approach42. This requires 9 parameters: an intercept (β1, β3, β5) and slope (β2, β4, β6) coefficient for each 1-way f parameter f1, f2, f3 and an intercept parameter for each detection parameter (β7, β8, β9). Below mentioned is the model for 1-way f parameters.$$f_{{1}} = , upbeta_{{{1 } + }} upbeta_{{2}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{7}}$$$$f_{{2}} = , upbeta_{{{3 } + }} upbeta_{{4}} left( {{text{Covariate}}} right),;;{text{ p}} = , upbeta_{{8}}$$$$f_{{3}} = , upbeta_{{5}} + , upbeta_{{6}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{9}} .$$All covariates were standardized before model fitting. We fitted the most complex model to each species and considered all possible combinations of covariates using the logit link function. Our rationale for including these variables in the occupancy and detectability component of the model was that we expected these variables to influence the occupancy and detectability of the study species.Since multi-species occupancy simultaneously model environmental variables, & interspecific interactions. Further it also allows to understand the influence of environmental variables on one species occupancy, in the presence or absence of other sympatric species64. Hence, we also modeled two species occur together as a function of covariates. We examined how the variables of each camera site influenced the pair-wise interaction of the three ungulate species. This model assumes that the conditional probability of one species varies in the presence or absence of other species. We assumed f123: co-occurrence of species 1, species 2 & species 3 = 0, hence we did not include higher-order interactions in any of our models, we assumed the conditional probability of 3 species occurred together was purely a function of species-specific (f1, f2, f3) and pair-wise interaction (f12: co-occurrence of species1 & species 2, f13: co-occurrence of species 1 & species 3, f23: co-occurrence of species 2 & species 3) parameters. We modeled pair-wise interaction of species varies as a function of environmental variables keeping detection probability constant. Hence, we specified 15 f and p parameters, an intercept and slope coefficient for each of the one-way (f1, f2, f3) and the two-way f parameters (f12, f13, and f23); as well as an intercept parameter for each of the detection models. The model equation below implies for 2-way f parameters:$$f_{{{12}}} = , upbeta_{{{7 } + }} upbeta_{{8}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{{13}}}$$$$f_{{{13}}} = , upbeta_{{{9 } + }} upbeta_{{{1}0}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{{14}}}$$$$f_{{{23}}} = , upbeta_{{{11 } + }} upbeta_{{{12}}} left( {{text{Covariate}}} right),;;{text{ p }} = , upbeta_{{{15}}} .$$We also fitted models including co-occurrence and detection probability of a species varies as a function of environmental variables. Hence, we specified 18 f and p parameters, an intercept and slope coefficient for each of one-way (f1, f2, f3) and two-way f parameters (f12, f13, f23); and an intercept as well as the slope parameters for each of the detection models. The model equation below implies for 2-way f parameters:$$f_{{{12}}} = , upbeta_{{{7 } + }} upbeta_{{8}} left( {{text{Covariate}}} right),{text{ p }} = , upbeta_{{{13 } + }} upbeta_{{{14}}} left( {{text{covariate}}} right)$$$$f_{{{13}}} = , upbeta_{{{9 } + }} upbeta_{{{1}0}} left( {{text{Covariate}}} right),{text{ p }} = , upbeta_{{{15}}} + , upbeta_{{{16}}} left( {{text{covariate}}} right)$$$$f_{{{23}}} = , upbeta_{{{11 } + }} upbeta_{{{12}}} left( {{text{Covariate}}} right),{text{ p }} = , upbeta_{{{17}}} + , upbeta_{{{18}}} left( {{text{covariate}}} right)$$A total of 38 models were run to test the influence of environmental variables on occupancy and detection probability of species-specific (f1, f2, f3) and pair-wise interaction of the three ungulate species. The best-supported model was identified by selecting the model with the lowest AICc value and highest model weights66, where higher model weights indicate a better fit of the model to the data. Second-Order Information Criterion (AICc)67 values were used to rank the occupancy models, and all the models whose ΔAICc  More