Philippa Kaur

Wu, A. Social buffering of stress – Physiological and ethological perspectives. Appl. Anim. Behav. Sci. 239, 105325 (2021).
Google Scholar 
Hennessy, M. B., Kaiser, S. & Sachser, N. Social buffering of the stress response: diversity, mechanisms, and functions. Front. Neuroendocrinol. 30, 470–482 (2009).CAS 
PubMed 

Google Scholar 
Young, C., Majolo, B., Heistermann, M., Schülke, O. & Ostner, J. Responses to social and environmental stress are attenuated by strong male bonds in wild macaques. Proc. Natl Acad. Sci. USA 111, 18195–18200 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
Stanton, M. E., Patterson, J. M. & Levine, S. Social influences on conditioned cortisol secretion in the squirrel monkey. Psychoneuroendocrinology 10, 125–134 (1985).CAS 
PubMed 

Google Scholar 
Caldji, C., Diorio, J. & Meaney, M. J. Variations in maternal care in infancy regulate the development of stress reactivity. Biol. Psychiatry 48, 1164–1174 (2000).CAS 
PubMed 

Google Scholar 
Novak, M. A., Hamel, A. F., Kelly, B. J., Dettmer, A. M. & Meyer, J. S. Stress, the HPA axis, and nonhuman primate well-being: a review. Appl. Anim. Behav. Sci. 143, 135–149 (2013).PubMed 
PubMed Central 

Google Scholar 
Sapolsky, R. M., Romero, L. M. & Munck, A. U. How do glucocorticoids influence stress responses? Integrating permissive, suppressive, stimulatory, and preparative actions. Endocr. Rev. 21, 55–89 (2000).CAS 
PubMed 

Google Scholar 
Liu, D. et al. Maternal Care, hippocampal glucocorticoid receptors, and hypothalamic-pituitary-adrenal responses to stress. Sci. Ment. Heal. Stress Brain 9, 75–78 (1997).
Google Scholar 
Gjerstad, J. K., Lightman, S. L. & Spiga, F. Role of glucocorticoid negative feedback in the regulation of HPA axis pulsatility. Stress 21, 403–416 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Spiga, F., Walker, J. J., Terry, J. R. & Lightman, S. L. HPA axis-rhythms. Compr. Physiol. 4, 1273–1298 (2014).PubMed 

Google Scholar 
Sapolsky, R. M. Why Zebras Don’t Get Ulcers (Henry Holt and Company, LLC, 2004).Campos, F. A. et al. Glucocorticoid exposure predicts survival in female baboons. Sci. Adv. 7, 1–10 (2021).
Google Scholar 
Banerjee, S. B., Arterbery, A. S., Fergus, D. J. & Adkins-Regan, E. Deprivation of maternal care has long-lasting consequences for the hypothalamic-pituitary-adrenal axis of zebra finches. Proc. R. Soc. B Biol. Sci. 279, 759–766 (2012).
Google Scholar 
Hennessy, M. B., Nigh, C. K., Sims, M. L. & Long, S. J. Plasma cortisol and vocalization responses of postweaning age guinea pigs to maternal and sibling separation: evidence for filial attachment after weaning. Dev. Psychobiol. 28, 103–115 (1995).CAS 
PubMed 

Google Scholar 
Hennessy, M. B., O’Leary, S. K., Hawke, J. L. & Wilson, S. E. Social influences on cortisol and behavioral responses of preweaning, periadolescent, and adult guinea pigs. Physiol. Behav. 76, 305–314 (2002).CAS 
PubMed 

Google Scholar 
Wiener, S. G., Johnson, D. F. & Levine, S. Influence of postnatal rearing conditions on the response of squirrel monkey infants to brief perturbations in mother-infant relationships. Physiol. Behav. 39, 21–26 (1987).CAS 
PubMed 

Google Scholar 
Girard-Buttoz, C. et al. Early maternal loss leads to short-but not long-term effects on diurnal cortisol slopes in wild chimpanzees. Elife 10, e64134 (2021).Rosenbaum, S. et al. Social bonds do not mediate the relationship between early adversity and adult glucocorticoids in wild baboons. Proc. Natl Acad. Sci. USA 117, 20052–20062 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Moss, C. Elephant Memories: Thirteen Years in the Life of an Elephant Family (Univ. Chicago Press, 1988).Douglas-Hamilton, I., Bhalla, S., Wittemyer, G. & Vollrath, F. Behavioural reactions of elephants towards a dying and deceased matriarch. Appl. Anim. Behav. Sci. 100, 87–102 (2006).
Google Scholar 
Shoshani, J., Kupsky, W. J. & Marchant, G. H. Elephant brain. Part I: gross morphology, functions, comparative anatomy, and evolution. Brain Res. Bull. 70, 124–157 (2006).PubMed 

Google Scholar 
Goldenberg, S. Z. & Wittemyer, G. Orphaned female elephant social bonds reflect lack of access to mature adults. Sci. Rep. 7, 14408 (2017).PubMed 
PubMed Central 

Google Scholar 
Goldenberg, S. Z. & Wittemyer, G. Orphaning and natal group dispersal are associated with social costs in female elephants. Anim. Behav. 143, 1–8 (2018).
Google Scholar 
Lee, P. C. Allomothering among African elephants. Anim. Behav. 35, 278–291 (1987).
Google Scholar 
Parker, J. M. et al. Poaching of African elephants indirectly decreases population growth through lowered orphan survival. Curr. Biol. 31, 4156–4162.e5 (2021).Wittemyer, G. et al. Where sociality and relatedness diverge: the genetic basis for hierarchical social organization in African elephants. Proc. R. Soc. B Biol. Sci. 276, 3513–3521 (2009).
Google Scholar 
Goldenberg, S. Z., Douglas-Hamilton, I. & Wittemyer, G. Vertical transmission of social roles drives resilience to poaching in elephant metworks. Curr. Biol. 26, 75–79 (2016).CAS 
PubMed 

Google Scholar 
Gobush, K. S., Mutayoba, B. M. & Wasser, S. K. Long-term impacts of poaching on relatedness, stress physiology, and reproductive output of adult female African elephants. Conserv. Biol. 22, 1590–1599 (2008).CAS 
PubMed 

Google Scholar 
Gobush, K. S. et al. Loxodonta africana (African Savanna Elephant). Loxodonta africana: the IUCN red list of threatened species 2021 e.T181008073A181022663 https://doi.org/10.2305/IUCN.UK.2021-1.RLTS.T181008073A181022663.en (2021).Wittemyer, G. et al. Illegal killing for ivory drives global decline in African elephants. Proc. Natl Acad. Sci. USA 111, 13117–13121 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
Wittemyer, G., Daballen, D. & Douglas-Hamilton, I. Comparative Demography of an At-Risk African Elephant Population. PLoS ONE 8, e53726 (2013).McCormick, S. D. & Romero, L. M. Conservation endocrinology. Bioscience 67, 429–442 (2017).
Google Scholar 
Wittemyer, G. The elephant population of Samburu and Buffalo Springs National Reserves, Kenya. Afr. J. Ecol. 39, 357–369 (2001).
Google Scholar 
Cockrem, J. F. Individual variation in glucocorticoid stress responses in animals. Gen. Comp. Endocrinol. 181, 45–58 (2013).CAS 
PubMed 

Google Scholar 
Taff, C. C., Schoenle, L. A. & Vitousek, M. N. The repeatability of glucocorticoids: a review and meta-analysis. Gen. Comp. Endocrinol. 260, 136–145 (2018).CAS 
PubMed 

Google Scholar 
Hooten, M. B. & Hobbs, N. T. A guide to Bayesian model selection for ecologists. Ecol. Monogr. 85, 3–28 (2015).
Google Scholar 
Wittemyer, G. & Getz, W. M. Hierarchical dominance structure and social organization in African elephants, Loxodonta africana. Anim. Behav. 73, 671–681 (2007).
Google Scholar 
Heim, C., Ehlert, U. & Hellhammer, D. H. The potential role of hypocortisolism in the pathophysiology of stress-related bodily disorders. Psychoneuroendocrinology 25, 1–35 (2000).CAS 
PubMed 

Google Scholar 
Dickens, M. J. & Romero, L. M. A consensus endocrine profile for chronically stressed wild animals does not exist. Gen. Comp. Endocrinol. 191, 177–189 (2013).CAS 
PubMed 

Google Scholar 
Ma, D., Serbin, L. A. & Stack, D. M. How children’s anxiety symptoms impact the functioning of the hypothalamus–pituitary–adrenal axis over time: a cross-lagged panel approach using hierarchical linear modeling. Dev. Psychopathol. 31, 1–15 (2018).
Google Scholar 
Blas, J., Bortolotti, G. R., Tella, J. L., Baos, R. & Marchant, T. A. Stress response during development predicts fitness in a wild, long lived vertebrate. Proc. Natl Acad. Sci. USA 104, 8880–8884 (2007).CAS 
PubMed 
PubMed Central 

Google Scholar 
Boonstra, R. Reality as the leading cause of stress: rethinking the impact of chronic stress in nature. Funct. Ecol. 27, 11–23 (2013).
Google Scholar 
Gunnar, M. R. & Vazquez, D. M. Low cortisol and a flattening of expected daytime rhythm: Potential indices of risk in human development. Dev. Psychopathol. 13, 515–538 (2001).CAS 
PubMed 

Google Scholar 
Perry, R. E. et al. Corticosterone administration targeting a hypo-reactive HPA axis rescues a socially-avoidant phenotype in scarcity-adversity reared rats. Dev. Cogn. Neurosci. 40, 100716 (2019).Fries, E., Hesse, J., Hellhammer, J. & Hellhammer, D. H. A new view on hypocortisolism. Psychoneuroendocrinology 30, 1010–1016 (2005).CAS 
PubMed 

Google Scholar 
Dorsey, C., Dennis, P., Guagnano, G., Wood, T. & Brown, J. L. Decreased baseline fecal glucocorticoid concentrations associated with skin and oral lesions in black rhinoceros (Diceros bicornis). J. Zoo. Wildl. Med. 41, 616–625 (2010).PubMed 

Google Scholar 
Pawluski, J. et al. Low plasma cortisol and fecal cortisol metabolite measures as indicators of compromised welfare in domestic horses (Equus caballus). PLoS ONE 12, 1–18 (2017).
Google Scholar 
Feng, X. et al. Maternal separation produces lasting changes in cortisol and behavior in rhesus monkeys. Proc. Natl Acad. Sci. USA 108, 14312–14317 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
González Ramírez, C. et al. The NR3C1 gene expression is a potential surrogate biomarker for risk and diagnosis of posttraumatic stress disorder. Psychiatry Res. 284, 112797 (2020).PubMed 

Google Scholar 
Cluver, L., Fincham, D. S. & Seedat, S. Posttraumatic stress in AIDS-orphaned children exposed to high levels of trauma: the protective role of perceived social support. J. Trauma. Stress 22, 106–112 (2009).PubMed 

Google Scholar 
Bastille-Rousseau, G. et al. Landscape-scale habitat response of African elephants shows strong selection for foraging opportunities in a human dominated ecosystem. Ecography 43, 149–160 (2020).
Google Scholar 
Foley, C. A. H., Papageorge, S. & Wasser, S. K. Noninvasive stress and reproductive measures of social and ecological pressures in free-ranging African elephants. Conserv. Biol. 15, 1134–1142 (2001).
Google Scholar 
Wittemyer, G., Getz, W. M., Vollrath, F. & Douglas-Hamilton, I. Social dominance, seasonal movements, and spatial segregation in African elephants: a contribution to conservation behavior. Behav. Ecol. Sociobiol. 61, 1919–1931 (2007).
Google Scholar 
Wittemyer, G., Daballen, D. & Douglas‐Hamilton, I. Differential influence of human impacts on age‐specific demography underpins trends in an African elephant population. Ecosphere 12, e03720 (2021).Brown, J. L. et al. Individual and environmental risk factors associated with fecal glucocorticoid metabolite concentrations in zoo-housed Asian and African elephants. PLoS ONE 14, 1–18 (2019).
Google Scholar 
Goldenberg, S. Z. et al. Increasing conservation translocation success by building social functionality in released populations. Glob. Ecol. Conserv. 18, e00604 (2019).Dantzer, B., Fletcher, Q. E., Boonstra, R. & Sheriff, M. J. Measures of physiological stress: a transparent or opaque window into the status, management and conservation of species? Conserv. Physiol. 2, 1–18 (2014).
Google Scholar 
Kaisin, O., Fuzessy, L., Poncin, P., Brotcorne, F. & Culot, L. A meta-analysis of anthropogenic impacts on physiological stress in wild primates. Conserv. Biol. 0, 1–14 (2020).CAS 

Google Scholar 
Ganswindt, A., Rasmussen, H. B., Heistermann, M. & Hodges, J. K. The sexually active states of free-ranging male African elephants (Loxodonta africana): defining musth and non-musth using endocrinology, physical signals, and behavior. Horm. Behav. 47, 83–91 (2005).CAS 
PubMed 

Google Scholar 
Santymire, R. M. et al. Using ACTH challenges to validate techniques for adrenocortical activity analysis in various African wildlife species. Int. J. Anim. Vet. Adv. 4, 99–108 (2012).CAS 

Google Scholar 
Watson, R. et al. Development of a versatile enzyme immunoassay for non-invasive assessment of glucocorticoid metabolites in a diversity of taxonomic species. Gen. Comp. Endocrinol. 186, 16–24 (2013).CAS 
PubMed 

Google Scholar 
Oduor, S. et al. Differing physiological and behavioral responses to anthropogenic factors between resident and non-resident African elephants at Mpala Ranch, Laikipia County, Kenya. PeerJ 8, e10010 (2020).Brown, J. L., Kersey, D. C., Freeman, E. W. & Wagener, T. Assessment of diurnal urinary cortisol excretion in Asian and African elephants using different endocrine methods. Zoo. Biol. 29, 274–283 (2010).PubMed 

Google Scholar 
Justice, C. O. et al. The moderate resolution imaging spectroradiometer (MODIS): land remote sensing for global change research. IEEE Trans. Geosci. Remote Sens. 36, 1228–1249 (1998).
Google Scholar 
Lafferty, D. J. R., Zimova, M., Clontz, L., Hackländer, K. & Mills, L. S. Noninvasive measures of physiological stress are confounded by exposure. Sci. Rep. 9, 1–6 (2019).
Google Scholar 
O’Dwyer, K., Dargent, F., Forbes, M. R. & Koprivnikar, J. Parasite infection leads to widespread glucocorticoid hormone increases in vertebrate hosts: a meta-analysis. J. Anim. Ecol. 89, 519–529 (2020).PubMed 

Google Scholar 
Parker, J. M., Goldenberg, S. Z., Letitiya, D. & Wittemyer, G. Strongylid infection varies with age, sex, movement and social factors in wild African elephants. Parasitology 147, 348–359 (2020).PubMed 

Google Scholar 
Gibbons, L., Jacobs, D. E., Fox, M. T. & Hansen, J. The RVC/FAO guide to veterinary diagnostic parasitology. McMaster egg-counting technique. http://www.rvc.ac.uk/review/Parasitology/EggCount/Purpose.htm (2004)R Core Team. A language and environment for statistical computing. https://www.r-project.org/. (2020).Rstudio Team. RStudio: integrated development for R. http://www.rstudio.com/ (2020).Plummer, M. rjags: Bayesian graphical models using MCMC. https://cran.r-project.org/package=rjags (2019).Brooks, S. P. & Gelman, A. General methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 7, 434–455 (1998).
Google Scholar 
Gelman, A. & Rubin, D. B. Inference from iterative simulation using multiple sequences. Stat. Sci. 7, 457–511 (1992).
Google Scholar 
Wickham, H. ggplot2: elegant graphics for data analysis. https://ggplot2.tidyverse.org (2016).Youngflesh, C. MCMCvis: tools to visualize, manipulate, and summarize MCMC output. J. Open Source Softw. 3, 640 (2018).
Google Scholar 
Parker, J. M. The Physiological Condition of Orphaned African Elephants (Loxodonta africana). Doctoral dissertation, Colorado State University. (2021). More

The management of introduced species, whether kudzu or zebra mussels, is costly and complex. Now, a paper reports a workable, effective solution that harnesses network analyses of ecological phenomena.Invasive species can pose severe economic and environmental problems, costing more than US$1 trillion worldwide since 1970 (ref. 1). Yet managing this human-driven issue is difficult in itself. The regions involved can be vast — entire continents or countries, for instance — while budgets are typically limited. As well, the sites potentially affected and management options can be numerous. Real systems (for example, all the lakes in the United States) can have thousands of locations that could potentially be infested. By contrast, considering just 40 locations means dealing theoretically with over 1 trillion unique combinations (240) where management could be applied (for instance, to reduce the number of invasive species leaving infested areas or entering uninfested ones). Given these constraints, a key problem is how and where to deploy control measures such as invasive-species removal. While sophisticated optimization approaches exist2, which use mathematical rules to exclude most suboptimal combinations and quickly zoom in to which locations should be managed to minimize new invasions, these algorithms are generally unfeasible for very large systems. Now, writing in Nature Sustainability, Ashander et al.3 demonstrate that simpler network metrics revealing linkages between patches can provide solutions that are often comparable to the more complex optimization algorithms. More

Webb, T. J., vanden Berghe, E. & O’Dor, R. Biodiversity’s big wet secret: The global distribution of marine biological records reveals chronic under-exploration of the deep pelagic ocean. PLoS ONE 5, https://doi.org/10.1371/journal.pone.0010223 (2010).Drazen, J. C. & Sutton, T. T. Dining in the Deep: The Feeding Ecology of Deep-Sea Fishes. Annual Review of Marine Science 9, 337–366, https://doi.org/10.1146/annurev-marine-010816-060543 (2017).ADS 
Article 
PubMed 

Google Scholar 
Brierley, A. S. Diel vertical migration. Current Biology 24, R1074–R1076, https://doi.org/10.1016/j.cub.2014.08.054 (2014).CAS 
Article 
PubMed 

Google Scholar 
Irigoien, X. et al. Large mesopelagic fishes biomass and trophic efficiency in the open ocean. Nature Communications 5, 10, https://doi.org/10.1038/ncomms4271 (2014).CAS 
Article 

Google Scholar 
Anderson, T. R. et al. Quantifying carbon fluxes from primary production to mesopelagic fish using a simple food web model. ICES Journal of Marine Science 76, 690–701, https://doi.org/10.1093/icesjms/fsx234 (2018).Article 

Google Scholar 
Saba, G. K. et al. Toward a better understanding of fish-based contribution to ocean carbon flux. Limnology and Oceanography 66, 1639–1664, https://doi.org/10.1002/lno.11709 (2021).ADS 
CAS 
Article 

Google Scholar 
Koslow, J. A., Kloser, R. J. & Williams, A. Pelagic biomass and community structure over the mid-continental slope off southeastern Australia based upon acoustic and midwater trawl sampling. Marine Ecology Progress Series 146, 21–35, https://doi.org/10.3354/meps146021 (1997).ADS 
Article 

Google Scholar 
Kaartvedt, S., Staby, A. & Aksnes, D. L. Efficient trawl avoidance by mesopelagic fishes causes large underestimation of their biomass. Marine Ecology Progress Series 456, 1–6, https://doi.org/10.3354/meps09785 (2012).ADS 
Article 

Google Scholar 
Lehodey, P., Murtugudde, R. & Senina, I. Bridging the gap from ocean models to population dynamics of large marine predators: A model of mid-trophic functional groups. Progress in Oceanography 84, 69–84, https://doi.org/10.1016/j.pocean.2009.09.008 (2010).ADS 
Article 

Google Scholar 
Van de Putte, A., Flores, H., Volckaert, F. & van Franeker, J. A. Energy content of Antarctic mesopelagic fishes: Implications for the marine food web. Polar Biology 29, 1045–1051, https://doi.org/10.1007/s00300-006-0148-z (2006).Article 

Google Scholar 
Stowasser, G. et al. Food web dynamics in the Scotia Sea in summer: A stable isotope study. Deep-Sea Research Part II-Topical Studies in Oceanography 59, 208–221, https://doi.org/10.1016/j.dsr2.2011.08.004 (2012).ADS 
CAS 
Article 

Google Scholar 
McCormack, S. A. et al. Decades of dietary data demonstrate regional food web structures in the Southern Ocean. Ecology and Evolution 11, 227–241, https://doi.org/10.1002/ece3.7017 (2021).Article 
PubMed 

Google Scholar 
Griffiths, S. P., Olson, R. J. & Watters, G. M. Complex wasp-waist regulation of pelagic ecosystems in the Pacific Ocean. Reviews in Fish Biology and Fisheries 23, 459–475, https://doi.org/10.1007/s11160-012-9301-7 (2013).Article 

Google Scholar 
Saunders, R. A., Hill, S. L., Tarling, G. A. & Murphy, E. J. Myctophid Fish (Family Myctophidae) Are Central Consumers in the Food Web of the Scotia Sea (Southern Ocean). Frontiers in Marine Science 6, https://doi.org/10.3389/fmars.2019.00530 (2019).Dornan, T., Fielding, S., Saunders, R. A. & Genner, M. J. Swimbladder morphology masks Southern Ocean mesopelagic fish biomass. Proceedings of the Royal Society B-Biological Sciences 286, 8, https://doi.org/10.1098/rspb.2019.0353 (2019).Article 

Google Scholar 
Freer, J. J., Tarling, G. A., Collins, M. A., Partridge, J. C. & Genner, M. J. Predicting future distributions of lanternfish, a significant ecological resource within the Southern Ocean. Diversity and Distributions 25, 1259–1272, https://doi.org/10.1111/ddi.12934 (2019).Article 

Google Scholar 
Hidalgo, M. & Browman, H. I. Developing the knowledge base needed to sustainably manage mesopelagic resources Introduction. ICES Journal of Marine Science 76, 609–615, https://doi.org/10.1093/icesjms/fsz067 (2019).Article 

Google Scholar 
Proud, R. et al. From siphonophores to deep scattering layers: Uncertainty ranges for the estimation of global mesopelagic fish biomass. ICES Journal of Marine Science 76, 718–733, https://doi.org/10.1093/icesjms/fsy037 (2019).Article 

Google Scholar 
Caccavo, J. A. et al. Productivity and Change in Fish and Squid in the Southern Ocean. Frontiers in Ecology and Evolution 9, https://doi.org/10.3389/fevo.2021.624918 (2021).Davison, P., Lara-Lopez, A. & Anthony Koslow, J. Mesopelagic fish biomass in the southern California current ecosystem. Deep-Sea Research Part II: Topical Studies in Oceanography 112, 129–142, https://doi.org/10.1016/j.dsr2.2014.10.007 (2015).ADS 
Article 

Google Scholar 
Pakhomov, E. & Yamamura, O. Report of the Advisory Panel on Micronekton Sampling Inter-calibration Experiment. Tech. Rep., PICES (2010).Cheung, W. W. L. et al. Projecting global marine biodiversity impacts under climate change scenarios. Fish and Fisheries 10, 235–251, https://doi.org/10.1111/j.1467-2979.2008.00315.x (2009).Article 

Google Scholar 
Saunders, R. A. & Tarling, G. A. Southern Ocean Mesopelagic Fish Comply with Bergmann’s Rule. American Naturalist 191, 343–351, https://doi.org/10.1086/695767 (2018).Article 

Google Scholar 
Proud, R., Cox, M. J. & Brierley, A. S. Biogeography of the Global Ocean’s Mesopelagic Zone. Current Biology 27, 113–119, https://doi.org/10.1016/j.cub.2016.11.003 (2017).CAS 
Article 
PubMed 

Google Scholar 
Robison, B. H. Conservation of Deep Pelagic Biodiversity. Conservation Biology 23, 847–858, https://doi.org/10.1111/j.1523-1739.2009.01219.x (2009).Article 
PubMed 

Google Scholar 
Constable, A. J. et al. Developing priority variables (“ecosystem Essential Ocean Variables” – eEOVs) for observing dynamics and change in Southern Ocean ecosystems. Journal of Marine Systems 161, 26–41, https://doi.org/10.1016/j.jmarsys.2016.05.003 (2016).ADS 
Article 

Google Scholar 
St John, M. A. et al. A Dark Hole in Our Understanding of Marine Ecosystems and Their Services: Perspectives from the Mesopelagic Community. Frontiers in Marine Science 3, 6, https://doi.org/10.3389/fmars.2016.00031 (2016).Article 

Google Scholar 
Newman, L. et al. Delivering Sustained, Coordinated, and Integrated Observations of the Southern Ocean for Global Impact. Frontiers in Marine Science 6, https://doi.org/10.3389/fmars.2019.00433 (2019).Costello, M. J. & Vanden Berghe, E. ‘Ocean biodiversity informatics’: a new era in marine biology research and management. Marine Ecology Progress Series 316, 203–214, https://doi.org/10.3354/meps316203 (2006).ADS 
Article 

Google Scholar 
Van de Putte, A. et al. From data to marine ecosystem assessments of the Southern Ocean, achievements, challenges, and lessons for the future. Frontiers in Marine Science 8, https://doi.org/10.3389/fmars.2021.637063 (2021).Duhamel, G. et al. Biogeographic Patterns of Fish. In Biogeographic Atlas of the Southern Ocean, 328–362 (Scientific Committee of Antarctic Research, Cambridge, UK, 2014).Piatkowski, U., Rodhouse, P. G., White, M. G., Bone, D. G. & Symon, C. Nekton community of the Scotia Sea as sampled by the RMT-25 during the austral summer. Marine Ecology Progress Series 112, 13–28, https://doi.org/10.3354/meps112013 (1994).ADS 
Article 

Google Scholar 
Collins, M. A. et al. Patterns in the distribution of myctophid fish in the northern Scotia Sea ecosystem. Polar Biology 31, 837–851, https://doi.org/10.1007/s00300-008-0423-2 (2008).Article 

Google Scholar 
Collins, M. A. et al. Latitudinal and bathymetric patterns in the distribution and abundance of mesopelagic fish in the Scotia Sea. Deep-Sea Research Part II-Topical Studies in Oceanography 59, 189–198, https://doi.org/10.1016/j.dsr2.2011.07.003 (2012).ADS 
Article 

Google Scholar 
Loeb, V. J., Hofmann, E. E., Klinck, J. M., Holm-Hansen, O. & White, W. B. ENSO and variability of the Antarctic Peninsula pelagic marine ecosystem. Antarctic Science 21, 135–148, https://doi.org/10.1017/s0954102008001636 (2009).ADS 
Article 

Google Scholar 
Reiss, C. S. et al. Overwinter habitat selection by Antarctic krill under varying sea-ice conditions: implications for top predators and fishery management. Marine Ecology Progress Series 568, 1–16, https://doi.org/10.3354/meps12099 (2017).ADS 
CAS 
Article 

Google Scholar 
Flores, H. et al. Distribution, abundance and ecological relevance of pelagic fishes in the Lazarev Sea, Southern Ocean. Marine Ecology Progress Series 367, 271–282, https://doi.org/10.3354/meps07530 (2008).ADS 
Article 

Google Scholar 
Flores, H. et al. Seasonal changes in the vertical distribution and community structure of Antarctic macrozooplankton and micronekton. Deep-Sea Research Part I-Oceanographic Research Papers 84, 127–141, https://doi.org/10.1016/j.dsr.2013.11.001 (2014).ADS 
Article 

Google Scholar 
Duhamel, G. The Pelagic Fish Community of the Polar Frontal Zone off the Kerguelen Islands. In Fishes of Antarctica, 63–74, https://doi.org/10.1007/978-88-470-2157-0_5 (Springer, Milano, 1998).Duhamel, G., Koubbi, P. & Ravier, C. Day and night mesopelagic fish assemblages off the Kerguelen Islands (Southern Ocean). Polar Biology 23, 106–112, https://doi.org/10.1007/s003000050015 (2000).Article 

Google Scholar 
Duhamel, G., Gasco, N. & Davaine, P. Poissons des îles Kerguelen et Crozet: Guide régional de l’océan Austral (Muséum national d’Histoire naturelle, Paris, 2005).Trebilco, R. et al. Mesopelagic community struture on the southern Kerguelen Axis. In The Kerguelen Plateau: Marine Ecosystem and Fisheries, 49–54 (Australian Antarctic Division, Kingston, Tasmania, 2019).Constable, A. J. & Swadling, K. M. Ecosystem drivers of food webs on the Kerguelen Axis of the Southern Ocean. Deep-Sea Research Part II-Topical Studies in Oceanography 174, 6, https://doi.org/10.1016/j.dsr2.2020.104790 (2020).Article 

Google Scholar 
Van de Putte, A. P., Jackson, G. D., Pakhomov, E., Flores, H. & Volckaert, F. A. M. Distribution of squid and fish in the pelagic zone of the Cosmonaut Sea and Prydz Bay region during the BROKE-West campaign. Deep-Sea Research Part II-Topical Studies in Oceanography 57, 956–967, https://doi.org/10.1016/j.dsr2.2008.02.015 (2010).ADS 
CAS 
Article 

Google Scholar 
Flynn, A. J. & Williams, A. Lanternfish (Pisces: Myctophidae) biomass distribution and oceanographic-topographic associations at Macquarie Island, Southern Ocean. Marine and Freshwater Research 63, 251–263, https://doi.org/10.1071/mf11163 (2012).Article 

Google Scholar 
Sutton, C. A., Kloser, R. J. & Gershwin, L. A. Micronekton in southeastern Australian and the Southern Ocean; A collation of the biomass, abundance, biodiversity and distribution data from CSIRO’s historical mesopelagic depth stratified new samples. CSIRO, Aust. http://hdl.handle.net/102.100.100/365479?index=1 (2018).Gon, O. & Heemstra, P. C. Fishes of the Southern Ocean (J.L.B. Smith Institute of Ichthyology, Grahamstown, South Africa, 1990).Darwin Core Maintenance Group. List of Darwin Core terms (2021).R Core Team. R: A language and environment for statistical computing (2021).Holstein, J. worms: Retrieving Aphia Information from World Register of Marine Species (2018).Bivand, R. et al. maptools: Tools for handling spatial objects. R package version 1.1-1 (2021).Orsi, A. H., Whitworth, T. & Nowlin, W. D. On the meridional extent and fronts of the Antarctic Circumpolar Current. Deep-Sea Research Part I-Oceanographic Research Papers 42, 641–673, https://doi.org/10.1016/0967-0637(95)00021-w (1995).ADS 
Article 

Google Scholar 
Constable, A. J. et al. Climate change and Southern Ocean ecosystems I: how changes in physical habitats directly affect marine biota. Global Change Biology 20, 3004–3025, https://doi.org/10.1111/gcb.12623 (2014).ADS 
Article 
PubMed 

Google Scholar 
Woods, B. et al. Myctobase. Zenodo https://doi.org/10.5281/zenodo.5590999 (2021).Saunders, R. A., Collins, M. A., Stowasser, G. & Tarling, G. A. Southern Ocean mesopelagic fish communities in the Scotia Sea are sustained by mass immigration. Marine Ecology Progress Series 569, 173–185, https://doi.org/10.3354/meps12093 (2017).ADS 
Article 

Google Scholar 
Provoost, P. & Bosch, S. obistools: Tools for data enhancement and quality control (2021).Murphy, E. J. et al. Understanding the structure and functioning of polar pelagic ecosystems to predict the impacts of change, https://doi.org/10.1098/rspb.2016.1646 (2016).McCormack, S. A., Melbourne-Thomas, J., Trebilco, R., Blanchard, J. L. & Constable, A. Alternative energy pathways in Southern Ocean food webs: Insights from a balanced model of Prydz Bay, Antarctica. Deep-Sea Research Part II-Topical Studies in Oceanography 174, https://doi.org/10.1016/j.dsr2.2019.07.001 (2020).Rodhouse, P. G. K. Role of squid in the Southern Ocean pelagic ecosystem and the possible consequences of climate change. Deep-Sea Research Part II-Topical Studies in Oceanography 95, 129–138, https://doi.org/10.1016/j.dsr2.2012.07.001 (2013).ADS 
CAS 
Article 

Google Scholar 
The MathWorks Inc., V.. MATLAB (2019).Potter, D. C., Lough, R. G., Perry, R. I. & Neilson, J. D. Comparison of the mocness and iygpt pelagic samplers for the capture of 0-group cod (gadus morhua) on georges bank. ICES Journal of Marine Science 46, https://doi.org/10.1093/icesjms/46.2.121 (1990).Elith, J., Leathwick, J. R. & Hastie, T. A working guide to boosted regression trees. Journal of Animal Ecology 77, 802–813, https://doi.org/10.1111/j.1365-2656.2008.01390.x (2008).CAS 
Article 
PubMed 

Google Scholar 
Oppel, S. et al. Comparison of five modelling techniques to predict the spatial distribution and abundance of seabirds. Biological Conservation 156, https://doi.org/10.1016/j.biocon.2011.11.013 (2012).McClatchie, S., Thorne, R. E., Grimes, P. & Hanchet, S. Ground truth and target identification for fisheries acoustics. Fisheries Research 47, 173–191, https://doi.org/10.1016/s0165-7836(00)00168-5 (2000).Article 

Google Scholar 
Collins, M., Piatkowski, U. & Saunders, R. A. Distribution of mesopelagic fish in the Scotia Sea from RMT25 and pelagic trawls deployed from RRS James Clark Ross and RRS John Biscoe, UK Polar Data Centre https://doi.org/10.5285/f4dfc0ee-4f61-47c5-a5a8-238e02ff2fdd (2021).Hoddell, R. J., Crossley, C., Hosie, G. & Williams, D. Fish and zooplankton from RMT-8 net hauls on the BROKE voyage. Australian Antarctic Data Centre https://doi.org/10.4225/15/57BA97EA8A22D (2016).Constable, A., Williams, D. & Lamb, T. Heard Island and McDonald Islands (HIMI) Marine Ecosystem. Australian Antarctic Data Centre https://doi.org/10.4225/15/5b31be45e8977 (2018).Van de Putte, A. Fish catches from Rectangular Midwater Trawl – data collected from the BROKE-West voyage of the Aurora Australis, 2006. Australian Antarctic Data Centre https://doi.org/10.4225/15/598d453109182 (2010).Flynn, A. J., Kloser, R. J. & Sutton, C. Micronekton assemblages and bioregional setting of the Great Australian Bight: A temperate northern boundary current system. Deep-Sea Research Part II: Topical Studies in Oceanography 157–158, https://doi.org/10.1016/j.dsr2.2018.08.006 (2018).Oozeki, Y., Hu, F., Tomatsu, C. & Kubota, H. Development of a new multiple sampling trawl with autonomous opening/closing net control system for sampling juvenile pelagic fish. Deep-Sea Research Part I-Oceanographic Research Papers 61, https://doi.org/10.1016/j.dsr.2011.12.001 (2012). More

Climate driversTo explore the impact of climate on forest resilience (see the following sections), we used monthly averaged total precipitation, 2-m air temperature, evapotranspiration deficit and surface solar radiation downwards acquired from the ERA5-Land reanalysis product at 0.1° spatial resolution for the 2000–2020 period (https://cds.climate.copernicus.eu/cdsapp#!/home). Evapotranspiration deficit was quantified as the total precipitation minus evapotranspiration. In this study, we referred to climate regions as defined by the Köppen–Geiger world map of climate classification51 (http://koeppen-geiger.vu-wien.ac.at/present.htm). The original 31 climatic zones were merged into major zones and only those characterized by vegetation cover were included in our study (tropical, arid, temperate and boreal; Extended Data Fig. 8).Vegetation dynamicsNDVI data acquired from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument aboard the Terra satellite was used to derive changes in global vegetation for the period 2000–2020. We used cloud-free spatial composites provided at 16-day temporal resolution and 0.05° spatial resolution (MOD13C1 Version 6; https://lpdaac.usgs.gov/products/mod13c1v006/) and retained only pixels with good and marginal overall quality. The MODIS-derived NDVI dataset represents a state-of-the-art product of vegetation state whose retrieval algorithm is constantly improved52, and being derived from a unique platform and sensor, it is temporally and spatially consistent. Vegetation dynamics were analysed in terms of kNDVI, a nonlinear generalization of the NDVI based on ref. 22 and derived as follows:$$text{kNDVI=}tanh left({text{NDVI}}^{2}right)$$
(1)
kNDVI has recently been proposed as a strong proxy for ecosystem productivity that shows high correlations with both plot level measurements of primary productivity and satellite retrievals of sun-induced fluorescence22. In addition, kNDVI has been documented to be more closely related to primary productivity, to be resistant to saturation, bias and complex phenological cycles, and to show enhanced robustness to noise and stability across spatial and temporal scales compared to alternative products (for example, NDVI and near-infrared reflectance of vegetation). For these reasons, it has been retained in this study as the preferred metric to describe the state of the forest ecosystem.To obtain an accurate estimate of resilience indicators, vegetation time series need to be stationary without seasonal periodic patterns or long-term trends53. To this aim, vegetation anomalies were obtained from kNDVI data by first subtracting the multi-year 16-day sample mean and then removing linear trends from the resulting time series. Missing data, due for instance to snow cover affecting the retrieval of reflectance properties, have been gap-filled by climatological kNDVI values. The time series of kNDVI-based vegetation anomalies was used to derive resilience indicators and assess their spatial and temporal variations (see next sections).Interannual changes in vegetation were assessed in terms of growing-season-averaged kNDVI. To this end, a climatological growing season that spanned months with at least 75% of days in the greenness phase was derived from the Vegetation Index and Phenology satellite-based product54 (https://vip.arizona.edu/) and acquired for the 2000–2016 period at 0.05° spatial resolution. In addition, forest cover (FC) fraction was derived from the annual land-cover maps of the European Space Agency’s Climate Change Initiative (https://www.esa-landcover-cci.org/)55 over the 2000–2018 period at 300-m spatial resolution. FC was retrieved by summing the fraction of broadleaved deciduous, broadleaved evergreen, needle leaf deciduous and needle leaf evergreen forest. FC was resampled to 0.05° to match the kNDVI spatial resolution.Spatial patterns of slowness and its dependence on environmental factorsIn this study, we quantified the resilience of forest ecosystems—their ability to recover from external perturbations—by the use of the 1-lag TAC (refs. 3,4,5). Such an indicator was initially computed on the whole time series of vegetation anomalies (2000–2020) for forest pixels with less than 50% missing data in the original NDVI and FC greater than 0.05 and referred to in the text as long-term TAC. This analysis was used to assess the spatial patterns of the forest slowness mediated by environmental factors that affect plant growth rates and capacity to recover from perturbations. The long-term TAC was explored both in the geographic and climate space (Extended Data Fig. 1). In the climate space, long-term TAC was binned in a 50 × 50 grid as a function of average annual precipitation and temperature, both computed over the 2000–2020 period, using the average as an aggregation metric weighted by the areal extents of each record. We retained only bins with at least 50 records.To explore the potential drivers of long-term TAC, we developed an RF regression model23 and predicted the observed long-term TAC (response variable) based on a set of environmental features (predictors). The use of machine learning in general and of RF in particular, being nonparametric and nonlinear data-driven methods, does not require a priori assumptions about the functional form relating the key drivers and the response functions. The environmental variables include vegetation properties (FC and growing-season-averaged kNDVI) and climate variables (total precipitation, 2-m air temperature, evapotranspiration deficit and surface solar radiation downwards). Each of the climate variables was expressed in terms of average, coefficient of variation and 1-lag autocorrelation and resampled to 0.05° spatial resolution to match the spatial resolution of kNDVI. All environmental variables were computed annually and then averaged over time, except the autocorrelation that was computed directly for the whole period, analogously to the long-term TAC. This resulted in a set of 14 predictors representing the forest density, the background climate, the climate variability and its TAC in the observational period (Extended Data Table 1). The RF model was developed by splitting the observed long-term TAC into two separate samples: 60% of records were used for model calibration, and the remaining 40% were used to validate model performances in terms of coefficient of determination (R2), mean squared error and percentage bias (PBIAS). Each record refers to a 0.05° pixel. The RF implemented here uses 100 regression trees, whose depth and number of predictors to sample at each node were identified using Bayesian optimization. The general model formulation is as follows:$$text{TAC},=,fleft(Xright)+{varepsilon }_{{rm{f}}}$$
(2)
in which f is the RF regression model, X are the environmental predictors and εf are the residuals. We found that the model explains 87% of the spatial variance (R2) of the observed long-term TAC with a mean squared error of 0.007 and an average overestimation of 0.058 (PBIAS; Extended Data Fig. 2a). By definition, machine learning methods are not based on the mechanistic representation of the phenomena and therefore cannot provide direct information on the underlying processes influencing the system response to drivers. However, some model-agnostic methods can be applied to gain insights into the outputs of RF models. Here we used variable importance metrics to quantify and rank how individual environmental factors influence TAC (Extended Data Fig. 2b). Furthermore, using partial dependence plots derived from the machine learning algorithm RF, we explored the ecosystem response function (TAC) across gradients of vegetation and climate features (Supplementary Discussion 1 and Extended Data Fig. 2c–f).CSD indicatorsTo explore the temporal variation in forest resilience, we used CSD indicators, here quantified in terms of temporal changes in TAC retrieved for two consecutive and independent periods ranging from 2000 to 2010 and from 2011 to 2020, and assessed the significance of the change in the sampled mean aggregated for different climate regions through a two-sided t-test (Fig. 1c). This analysis was complemented by the computation of TAC on the annual scale over a 2-year lagged temporal window (3-year window size) to track the temporal changes in CSD. This resulted in a time series of TAC with an annual time step.We point out that temporal dynamics of annual TAC are driven by two processes: the changes in the resilience of the system that affect the velocity of the recovery from external perturbations and the confounding effects of the changes in autocorrelation of the climate drivers (Xac) that directly affect the autocorrelation of NDVI. Given the specific goals of this study, we factored out the second process from the total TAC signal to avoid that an increasing autocorrelation in the drivers would affect our analysis and conclusions about the resilience and the potential increase in instability56. For this purpose, we disentangled the temporal changes in TAC due to variations in autocorrelation in the climate drivers (({rm{TAC}}| {X}_{{rm{ac}}})) by adopting the space-for-time analogy and applied the RF model (f) at an annual time step (t) in a set of factorial simulations as follows:$${text{TAC}}^{t},{rm{| }},{X}_{{rm{ac}}}=fleft({X}^{t}right)-fleft({X}_{-{rm{ac}}}^{t},{X}_{{rm{ac}}}^{2000}right)$$
(3)
The first term on the right side of equation (3) is the RF model simulation obtained by accounting for the dynamics of all predictors, and the second term is the RF model simulation generated by considering all predictors dynamic except the factors of autocorrelation in climate that are kept constant to their first-year value (year 2000). For such runs, we used predictors computed on an annual scale over a 2-year lagged temporal window, consistently to the TAC time series. We found that the direct effects of autocorrelation in climate have led to a positive trend of TAC in dry zones (due to the increasing autocorrelation of the drivers in these regions) and to an opposite effect in temperate humid forests (Supplementary Fig. 3). To remove these confounding effects, the estimated term ({{rm{TAC}}}^{t}| {X}_{{rm{ac}}}) is factored out from the TACt by subtraction to derive an enhanced estimate of annual resilience that is independent of autocorrelation in climate (Extended Data Fig. 3).Long-term linear trends computed on the resulting enhanced TAC time series (δTAC) represent our reference CSD indicator used in this study to explore the changes in forest resilience. δTAC was quantified for each grid cell (Fig. 1a) and represented in the climate space following the methodology previously described (Fig. 1b). We then assessed the significance of the trends at bin level by applying a two-sided t-test for the sampled trend distributions within each bin. This significance test is independent from the structural temporal dependencies originating from the use of a 2-year lagged temporal window to compute the TAC time series.Following an analogous approach described in equation (3), we disentangled the effect of the variation in forest density, background climate and climate variability on temporal changes in TAC (Fig. 1d,e). We recognize that other environmental factors not explicitly accounted for in our RF model could play a role in modulating the temporal variations in TAC. However, given the comprehensiveness of the suite of predictors used in equation (2) (Extended Data Table 1), it seems plausible that residuals mostly reflect the intrinsic forest resilience, the component intimately connected to the short-term responses of forests to perturbations, which is not directly related to climate variability. Forest ecosystem evolutionary processes could also play a role, but longer time series would be required to reliably capture these dynamics. Furthermore, abrupt declines (ADs) in the vegetation state and following recoveries, similarly to those potentially originating from forest disturbances (for example, wildfires and insect outbreaks), could influence the TAC changes. However, such occurrences, being distributed across the globe throughout the whole period, are expected to only marginally affect the resulting trend in TAC time series.Sensitivity analysisTo assess the robustness of our results with respect to the modelling choices described above, we performed a series of sensitivity analyses for the difference in TAC retrieved for the two independent periods (2000–2010 and 2011–2020). To this aim, we tested their dependence on: the quality flag of the NDVI data used for the analyses (good, good and marginal); the gap-filling procedure tested on different periods (year and growing season); the inclusion or exclusion of forest areas affected by ADs; the threshold on the maximum percentage of missing NDVI data allowed at the pixel level (20%, 50% and 80%); the threshold on the minimum percentage of FC allowed at the pixel level (5%, 50% and 90%); and the pixel spatial resolution used for the analyses (0.05°, 0.25° and 1°). In addition, we tested the sensitivity of the trend in total TAC signal on the moving temporal window length used to calculate autocorrelation at lag 1. Results obtained for the different configurations were compared in terms of frequency distributions, separately for climate regions (Extended Data Fig. 4), and further explored in the climate space (Extended Data Figs. 5 and 6). Outcomes of the sensitivity analysis are discussed in Supplementary Discussion 2.Interplay between GPP and CSDResilience and GPP interact with each other through mutual causal links. On one hand, a reduction in forest resilience makes the system more sensitive to perturbations with potential consequent losses in GPP (ref. 26). On the other hand, a reduction in GPP may lead to a decline in resilience according to the carbon starvation hypothesis, and may be associated with increasing hydraulic failure46. To explore the link between forest resilience and primary productivity, we quantified the correlation between TAC and GPP. Estimates of GPP were derived from the FluxCom Model Tree Ensemble for the 2001–2019 period at 8-daily temporal resolution and 0.0833° spatial resolution and generated using ecosystem GPP fluxes from the FLUXNET network and MODIS remote sensing data as predictor variables36 (http://www.fluxcom.org/). Annual maps of GPP were quantified and resampled to 0.05° to match the temporal and spatial resolution of TAC time series. The Spearman rank correlation (ρ) was then computed between annual GPP and TAC over a 1° spatial moving window to better sample the empirical distribution of the two variables (Fig. 2d). The significance of ρ(GPP,TAC) was assessed over the climate space separately for each bin (Fig. 2e), similarly to the approach used to test the significance of δTAC. Furthermore, we explored the relationships between the trend in GPP (δGPP) and the trend in TAC (δTAC) by clustering the globe according to the directions of the long-term trajectories of the above-mentioned variables (Fig. 2f).Disentangling the impact of forest managementTo characterize TAC on different forest types and disentangle the potential effects originating from forest management, results were separately analysed for intact forests and managed forests. Intact forests were considered those forest pixels constituting the Intact Forest Landscapes57 dataset (https://intactforests.org/). Intact Forest Landscapes identifies the forest extents with no sign of significant human activity over the period 2000–2016 based on Landsat time series. The remaining forests pixels—not labelled as intact—were considered as managed forests (Extended Data Fig. 8). The resulting forest type map is consistent with those used for United Nations Framework Convention on Climate Change reporting58, although with more conservative estimates of intact forests in the boreal zone due to the masking based on FC and percentage of missing data applied in this study.We analysed the differences in long-term TAC (computed for the whole 2000–2020 period) between managed and intact forests by masking out the potential effect of climate background. To this aim, we compared the climate spaces generated separately for managed and intact forests by extracting only those bins that are covered by both forest classes. The resulting distributions—one for each forest class—have the same sample size, and each pair of elements shares the same climate background. Potential confounding environmental effects on average recovery rates are, therefore, minimized. We then applied a two-sided t-test for analysing the significance of the difference in the sampled means (Fig. 2a). An analogous approach was used to test the differences in δTAC and ρ(GPP,TAC) between managed and intact forests (Fig. 2b,c).Early-warning signals of abrupt forest declinesWhen forest ecosystems are subject to an extended and progressive degradation, the loss of resilience can lead to an AD (refs. 3,4,5). Such abrupt changes can trigger a regime shift (tipping point) depending on the capacity of the system to recover from the perturbations (Supplementary Methods 1 and 2). We investigated the potential of changes in TAC as early-warning signals of ADs in intact forests over the 2010–2020 period. To this aim, we quantified at the pixel level ADs as the events occurring on a certain year when the corresponding growing-season average kNDVI was more than n-times local standard deviation below the local mean. Local mean and standard deviation (σ) were computed over the 10-year antecedent temporal window (undisturbed) period and n varies between 1 and 6 with higher values reflecting more severe changes in the state of the system. For each pixel and for each fixed n value, we recorded only the first AD occurrence, thus imposing a univocal record for each abrupt change in the state of the system.We then explored whether the retrieved ADs were statistically associated with antecedent high values of δTAC. To avoid confusion with the attribution of causality, for each AD that occurred at time t (over the 2010–2020 period), we derived the δTAC over the temporal window 2000 − (t − 1). The resulting trend in δTAC is therefore antecedent and independent of the changes in vegetation associated with the AD. Then, for each pixel with an AD at time t, we also extracted randomly one of the undisturbed (with no AD) adjacent pixels and retrieved δTAC over the same temporal window. This analysis produced two distributions of δTAC associated with pixels with and without ADs (AD and no AD, respectively). The two distributions have the same size and each pair of elements shares similar background climate. We calculated the probability of occurrence of AD conditional on the trend in δTAC (({rm{AD}}| delta {rm{TAC}})) as the frequency of ADs for which (delta {rm{TAC}}left(mathrm{AD}right)| > delta {rm{TAC}}left(mathrm{no; AD}right)), and the significance of the difference in the two sampled means (AD and no AD) was evaluated through a two-sided t-test. Probability and significance were assessed for different climate regions and severity of ADs (Fig. 3a). High statistically significant probabilities suggest that the AD is following the drifting towards a critical resilience threshold plausibly associated with changes in environmental drivers.We complemented the aforementioned analyses by retrieving the tolerance and proximity to AD, hereafter determined for a 3σ severity. We first quantified the TAC that proceeded the occurrence of an AD and followed a progressive loss of resilience as captured by positive δTAC. This value, hereafter referred to as abrupt decline temporal autocorrelation (TACAD), reflects the TAC threshold over which we observed an abrupt change in the forest state (Fig. 3b). The tolerance to AD was quantified as the difference between the local TACAD and the TAC value averaged over the 2000–2009 period to characterize the pre-disturbance conditions. The tolerance metric was explored across a gradient of aridity index59 (Fig. 3c).TACAD can be directly retrieved only on those forest pixels that have already experienced an AD. As a considerable fraction of undisturbed forests could potentially be close to their critical TAC threshold, or even have already passed it, it is important to determine their TACAD. To this aim, we developed an RF regression model that expresses the TACAD as a function of the set X of environmental variables used in model f (equation (2)) but excluding the autocorrelation in climate drivers (Xreduced) already disentangled in the TAC signal. The general formulation is as follows:$${{rm{TAC}}}_{{rm{AD}}}=gleft({X}_{text{reduced}}right)+{varepsilon }_{{rm{g}}}$$
(4)
in which g is the RF regression model, Xreduced are the environmental predictors and εg are the residuals. Implementation, calibration and validation of g follow the same rationale described before for the f model. We found that the RF model explains 50% of the variance (R2) of the observed TACAD, with a mean squared error of 0.019 and an average underestimation of 0.86 (PBIAS).The RF model was then used to predict the TACAD over the whole domain of intact forests and served as input to quantify the proximity to AD of undisturbed forest pixels at the end of the observational period (year 2020). Here we defined the proximity metric as the difference between the value of TAC in 2020 and TACAD. Proximity takes negative or zero values when TACAD has already been reached (({{{rm{TAC}}}^{2020}ge {rm{TAC}}}_{{rm{AD}}})) and positive values when there are still margins before reaching the critical threshold (({{{rm{TAC}}}^{2020} < {rm{TAC}}}_{{rm{AD}}})). Together (delta {rm{TAC}} > 0) and ({{{rm{TAC}}}^{2020}ge {rm{TAC}}}_{{rm{AD}}}) therefore represent the most critical conditions, as they indicate that the critical resilience threshold for AD has already been reached and the ecosystem is continuing to lose its capacity to respond to external perturbations. We finally quantified the amount of GPP potentially exposed to such critical conditions by linearly extrapolating the GPP for the year 2020 (available GPP data stop in 2019) and overlaying it on the map of critical conditions (proximity to ({rm{AD}} < 0) and (delta {rm{TAC}} > 0)).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this paper. More

Karlsson, O. et al. Pesticide-induced multigenerational effects on amphibian reproduction and metabolism. Sci. Total Environ. 775, 145771 (2021).ADS 
CAS 
PubMed 
Article 

Google Scholar 
IUCN. The IUCN Red List of Threatened Species. Version 2021-3. https://www.iucnredlist.org (2022).Wake, D. B. & Koo, M. S. Amphibians. Curr. Biol. 28, R1237–R1241 (2018).CAS 
PubMed 
Article 

Google Scholar 
Campbell Grant, E. H., Miller, D. A. & Muths, E. A synthesis of evidence of drivers of amphibian declines. Herpetologica 76, 101–107 (2020).Article 

Google Scholar 
Green, D. M., Lannoo, M. J., Lesbarrères, D. & Muths, E. Amphibian population declines: 30 years of progress in confronting a complex problem. Herpetologica 76, 97–100 (2020).Article 

Google Scholar 
Mason, R., Tennekes, H., Sánchez-Bayo, F. & Jepsen, P. U. Immune suppression by neonicotinoid insecticides at the root of global wildlife declines. J. Environ. Immunol. Toxicol. 1, 3–12 (2013).Article 

Google Scholar 
Adams, E., Leeb, C. & Brühl, C. A. Pesticide exposure affects reproductive capacity of common toads (Bufo bufo) in a viticultural landscape. Ecotoxicology 30, 213–223 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Frost, D. R. Amphibian species of the world 6,1, an online reference. Electron. Datab. https://doi.org/10.5531/db.vz.0001 (American Museum of Natural History, 2021).Article 

Google Scholar 
Eterovick, P. C., Souza, A. M. & Sazima, I. Anfíbios da Serra do Cipó [Amphibians from the Serra do Cipó]. http://herpeto.org/wp-content/uploads/2020/11/ANFIBIOS-DA-SERRA-DO-CIPO.pdf (PUCMINAS, 2020).Mijares, A., Rodrigues, M. T. & Baldo, D. Physalaemus cuvieri The IUCN Red List of Threatened Species, version 2014.3. http://www.iucnredlist.org (2010). Accessed 9 Jan 2015.de Sá, F. P., Zina, J. & Haddad, C. F. B. Reproductive dynamics of the Neotropical treefrog Hypsiboas albopunctatus (Anura, Hylidae). J. Herpetol. 48, 181–185 (2014).Article 

Google Scholar 
Herek, J. S. et al. Can environmental concentrations of glyphosate affect survival and cause malformation in amphibians? Effects from a glyphosate-based herbicide on Physalaemus cuvieri and P. gracilis (Anura: Leptodactylidae). Environ. Sci. Pollut. Res. 27, 22619–22630 (2020).CAS 
Article 

Google Scholar 
Silva, F. L. et al. Swimming ability in tadpoles of Physalaemus cf. cuvieri, Scinax x-signatus and Leptodactylus latrans (Amphibia: Anura) exposed to the insecticide chlorpyrifos. Ecotoxicol. Environ. Contam. 16, 13–18 (2021).
Google Scholar 
Pavan, F. A. et al. Morphological, behavioral and genotoxic effects of glyphosate and 2,4-D mixture in tadpoles of two native species of South American amphibians. Environ. Toxicol. Pharmacol. 85, 103637 (2021).CAS 
PubMed 
Article 

Google Scholar 
Simon-Delso, N. et al. Systemic insecticides (Neonicotinoids and fipronil): Trends, uses, mode of action and metabolites. Environ. Sci. Pollut. Res. 22, 5–34 (2015).CAS 
Article 

Google Scholar 
Pietrzak, D., Kania, J., Malina, G., Kmiecik, E. & Wątor, K. Pesticides from the EU first and second watch lists in the water environment. Clean 47, 1–10 (2019).
Google Scholar 
IBAMA: Instituto Brasileiro do Meio Ambiente e dos Recursos Naturais Renováveis. Relatório de comercialização de agrotóxicos 2019 [Brazilian Pesticide Marketing Report 2019] https://www.ibama.gov.br/agrotoxicos/relatorios-de-comercializacao-de-agrotoxicos#boletinsanuais (2021).IBAMA: Instituto Brasileiro do Meio Ambiente e dos Recursos Naturais Renováveis. Vendas de ingredientes ativos por UF [Active ingredient sales by UF in Brazil]. http://ibama.gov.br/phocadownload/qualidadeambiental/relatorios/2019/Vendas_ingredientes_ativos_UF_2019.x (2021).IBAMA – Instituto Brasileiro do Meio Ambiente e dos Recursos Naturais Renováveis. Boletins anuais de produção, importação, exportação e vendas de agrotóxicos no Brasil [Annual bulletins of production, import, export and sales of pesticides in Brazil]. http://ibama.gov.br/index.php?option=com_content&view=article&id=594&Itemid=54 (2021).Pietrzak, D., Kania, J., Kmiecik, E., Malina, G. & Wątor, K. Fate of selected neonicotinoid insecticides in soil–water systems: Current state of the art and knowledge gaps. Chemosphere 255, 126981 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
ANVISA: Agência Nacional de Vigilância Sanitária; Índice Monográfico I13. Imidacloprido. http://portal.anvisa.gov.br/documents/111215/117782/I13+%E2%80%93+Imidacloprido/9d08c7e5-8979-4ee9-b76c-1092899514d7 (2021).Kagabu, S. Discovery of imidacloprid and further developments from strategic molecular designs. J. Agric. Food Chem. 59, 2887–2896 (2011).CAS 
PubMed 
Article 

Google Scholar 
Tomizawa, M. & Casida, J. E. Neonicotinoid insecticide toxicology: Mechanisms of selective action. Annu. Rev. Pharmacol. Toxicol. 45, 247–268 (2005).CAS 
PubMed 
Article 

Google Scholar 
Hashimoto, F. et al. Occurrence of imidacloprid and its transformation product (imidacloprid-nitroguanidine) in rivers during an irrigating and soil puddling duration. Microchem. J. 153, 12 (2020).Article 
CAS 

Google Scholar 
Hladik, M. L. et al. Year-round presence of neonicotinoid insecticides in tributaries to the Great Lakes, USA. Environ. Pollut. 235, 1022–1029 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Jurado, A., Walther, M. & Díaz-Cruz, M. Occurrence, fate and environmental risk assessment of the organic microcontaminants included in the Watch Lists set by EU Decisions 2015/495 and 2018/840 in the groundwater of Spain. Sci. Total Environ. 663, 285–296 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Montagner, C. C. et al. Ten years-snapshot of the occurrence of emerging contaminants in drinking, surface and ground waters and wastewaters from São Paulo State, Brazil. J. Braz. Chem. Soc. 30, 614–632 (2019).CAS 

Google Scholar 
CCME. Council of Ministers of the Environment. Canadian water quality guidelines for the protection of aquatic life. Imidacloprid. In Canadian water quality guidelines, Council of Ministers of the Environment. Winnipeg. https://ccme.ca/en/res/imidacloprid-en-canadian-water-quality-guidelines-for-the-protection-of-aquatic-life.pdf (2007).RIVM. Water quality standards for imidacloprid: Proposal for an update according to the Water Framework Directive in National Institute for Public Health and the Environment. https://www.rivm.nl/bibliotheek/rapporten/270006001.pdf (2014).PAN. Pesticide Action Network. International Consolidated List of Banned Pesticides. https://pan-international.org/pan-international-consolidated-list-of-banned-pesticides/ (2021).Brazil. Secretaria Estadual da Saúde do Rio Grande do Sul. Portaria SES RS nº 320, de 28 de abril de 2014. https://www.cevs.rs.gov.br/upload/arquivos/201705/11110603-portaria-agrotoxicos-n-320-de-28-de-abril-de-2014.pdf. (2014).Kobashi, K. et al. Comparative ecotoxicity of imidacloprid and dinotefuran to aquatic insects in rice mesocosms. Ecotoxicol. Environ. Saf. 138, 122–129 (2017).CAS 
PubMed 
Article 

Google Scholar 
Islam, M. A., Hossen, M. S., Sumon, K. A. & Rahman, M. M. Acute toxicity of imidacloprid on the developmental stages of common carp Cyprinus carpio. Toxicol. Environ. Health Sci. 11, 244–251 (2019).Article 

Google Scholar 
Pérez-Iglesias, J. M. et al. The genotoxic effects of the imidacloprid-based insecticide formulation Glacoxan Imida on Montevideo tree frog Hypsiboas pulchellus tadpoles (Anura, Hylidae). Ecotoxicol. Environ. Saf. 104, 120–126 (2014).PubMed 
Article 
CAS 

Google Scholar 
Sievers, M., Hale, R., Swearer, S. E. & Parris, K. M. Contaminant mixtures interact to impair predator-avoidance behaviours and survival in a larval amphibian. Ecotoxicol. Environ. Saf. 161, 482–488 (2018).CAS 
PubMed 
Article 

Google Scholar 
USEPA. United States Environmental Protection Agency. Aquatic Life Benchmarks and Ecological Risk Assessments for Registered Pesticides. https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/aquatic-life-benchmarks-and-ecological-risk. (2021).Feng, S., Kong, Z., Wang, X., Zhao, L. & Peng, P. Acute toxicity and genotoxicity of two novel pesticides on amphibian, Rana N. Hallwell. Chemosphere 56, 457–463 (2004).ADS 
CAS 
PubMed 
Article 

Google Scholar 
De Arcaute, C. R. et al. Genotoxicity evaluation of the insecticide imidacloprid on circulating blood cells of Montevideo tree frog Hypsiboas pulchellus tadpoles (Anura, Hylidae) by comet and micronucleus bioassays. Ecol. Indic. 45, 632–639 (2014).Article 
CAS 

Google Scholar 
Nkontcheu, D. B. K., Tchamadeu, N. N., Ngealekeleoh, F. & Nchase, S. Ecotoxicological effects of imidacloprid and lambda-cyhalothrin (insecticide) on tadpoles of the African common toad, Amietophrynus regularis (Reuss, 1833) (Amphibia: Bufonidae). Emerg. Sci. J. 1, 49–53 (2017).
Google Scholar 
Bortoluzzi, E. C. et al. Contaminação de águas superficiais por agrotóxicos em função do uso do solo numa microbacia hidrográfica de Agudo, RS. Rev. Bras. Eng. Agric. Ambient. 10, 881–887 (2006).Article 

Google Scholar 
Bortoluzzi, E. C. et al. Investigation of the occurrence of pesticide residues in rural wells and surface water following application to tobacco. Quim. Nova 30, 1872–1876 (2007)CAS 
Article 

Google Scholar 
La, N., Lamers, M., Bannwarth, M., Nguyen, V. V. & Streck, T. Imidacloprid concentrations in paddy rice fields in northern Vietnam: measurement and probabilistic modeling. Paddy Water Environ. 13, 191–203 (2015).Article 

Google Scholar 
Sweeney, M. R., Thompson, C. M. & Popescu, V. D. Sublethal, behavioral, and developmental effects of the neonicotinoid pesticide imidacloprid on larval wood frogs (Rana sylvatica). Environ. Toxicol. Chem. 40, 1838–1847 (2021).Article 
CAS 

Google Scholar 
Gibbons, D., Morrissey, C. & Mineau, P. A review of the direct and indirect effects of neonicotinoids and fipronil on vertebrate wildlife. Environ. Sci. Pollut. Res. 22, 103–118 (2015).CAS 
Article 

Google Scholar 
Morrissey, C. A. et al. Neonicotinoid contamination of global surface waters and associated risk to aquatic invertebrates: A review. Environ. Int. 74, 150920 (2015).Article 
CAS 

Google Scholar 
Stinson, S. A. et al. Agricultural surface water, imidacloprid, and chlorantraniliprole result in altered gene expression and receptor activation in Pimephales promelas. Sci. Total Environ. 806, 150920. (2022).ADS 
CAS 
PubMed 
Article 

Google Scholar 
DiGiacopo, D. G. & Hua, J. Evaluating the fitness consequences of plasticity in tolerance to pesticides. Ecol. Evol. 10, 4448–4456 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Carlson, B. E. & Langkilde, T. Body size variation in aquatic consumers causes pervasive community effects, independent of mean body size. Ecol. Evol. 7, 9978–9990 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Phung, T. X., Nascimento, J. C. S., Novarro, A. J. & Wiens, J. J. Correlated and decoupled evolution of adult and larval body size in frogs. Proc. Royal Soc. B 287, 20201474 (2020).Article 

Google Scholar 
Beasley, V. R. Direct and indirect effects of environmental contaminants on amphibians. In Reference Module in Earth Systems and Environmental Sciences https://doi.org/10.1016/b978-0-12-409548-9.11274-6 (Elsevier, 2020).Toledo, L. F., Sazima, I. & Haddad, C. F. B. Behavioural defences of anurans: An overview. Ethol. Ecol. Evol. 23, 1–25 (2011).Article 

Google Scholar 
Hartmann, M. T., Hartmann, P. A. & Haddad, C. F. B. Reproductive modes and fecundity of an assemblage of anuran amphibians in the Atlantic rainforest, Brazil. Inheringia 100, 207–215 (2010).Article 

Google Scholar 
Pupin, N. C., Gasparini, J. L., Bastos, R. P., Haddad, C. F. B. & Prado, C. P. A. Reproductive biology of an endemic Physalaemus of the Brazilian Atlantic forest, and the trade-off between clutch and egg size in terrestrial breeders of the P. signifer group. Herpetol. J. 20, 147–156 (2010).
Google Scholar 
Pereira, G. & Maneyro, R. Size-fecundity relationships and reproductive investment in females of Physalaemus riograndensis Milstead, 1960 (Anura, Leiuperidae) in Uruguay. Herpetol. J. 22, 145–150 (2012).
Google Scholar 
Tolledo, J., Silva, E. T., Nunes-de-Almeida, C. H. L. & Toledo, L. F. Anomalous tadpoles in a Brazilian oceanic archipelago: implications of oral anomalies on foraging behaviour, food intake and metamorphosis. Herpetol. J. 24, 237–243 (2014).
Google Scholar 
Annibale, F. S. et al. Smooth, striated, or rough: how substrate textures affect the feeding performance of tadpoles with different oral morphologies. Zoomorphology 139, 97–110 (2020).Article 

Google Scholar 
Venesky, M. D., Wassersug, R. J. & Parris, M. J. The impact of variation in labial tooth number on the feeding kinematics of tadpoles of southern leopard frog (Lithobates sphenocephalus). Copeia 3, 481–486 (2010).Article 

Google Scholar 
Venesky, M. D. et al. Comparative feeding kinematics of tropical hylid tadpoles. J. Exp. Biol. 216, 1928–1937 (2013).PubMed 

Google Scholar 
Jones, S. K. C., Munn, A. J., Penman, T. D. & Byrne, P. G. Long-term changes in food availability mediate the effects of temperature on growth, development and survival in striped marsh frog larvae: implications for captive breeding programmes. Conserv. Physiol. 3, cov029 (2015).Article 
CAS 

Google Scholar 
Bach, N. C., Natale, G. S., Somoza, G. M. & Ronco, A. E. Effect on the growth and development and induction of abnormalities by a glyphosate commercial formulation and its active ingredient during two developmental stages of the South-American Creole frog, Leptodactylus latrans. Environ. Sci. Pollut. Res. 23, 23959–23971 (2016).CAS 
Article 

Google Scholar 
Capellán, E. & Nicieza, A. G. Non-equivalence of growth arrest induced by predation risk or food limitation: context-dependent compensatory growth in anuran tadpoles. J. Anim. Ecol. 76, 1026–1035 (2007).PubMed 
Article 

Google Scholar 
Chin, A. M., Hill, D. R., Aurora, M. & Spence, J. R. Morphogenesis and maturation of the embryonic and postnatal intestine. Semin. Cell Dev. Biol. 66, 81–93 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Sun, Y., Zhang, J., Song, W. & Shan, A. Vitamin E alleviates phoxim-induced toxic effects on intestinal oxidative stress, barrier function, and morphological changes in rats. Environ. Sci. Pollut. Res. 25, 26682–26692 (2018).
Google Scholar 
Ouellet, M. Amphibian deformities: current state of knowledge. In Ecotoxicology of Amphibians and Reptiles (eds Sparling, D. W. et al.) 617–661 (Society of Environmental Toxicology and Chemistry, 2000).Hussein, M. & Singh, V. Effect on chick embryos development after exposure to neonicotinoid insecticide imidacloprid. J. Anat. Soc. India 65, 83–89 (2016).Article 

Google Scholar 
Crosby, E. B., Bailey, J. M., Oliveri, A. N. & Levin, E. D. Neurobehavioral impairments caused by developmental imidacloprid exposure in zebrafish. Neurotoxicol. Teratol. 49, 81–90 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Lonare, M. et al. Evaluation of imidacloprid-induced neurotoxicity in male rats: A protective effect of curcumin. Neurochem. Int. 78, 122–129 (2014).CAS 
PubMed 
Article 

Google Scholar 
Žegura, B., Lah, T. T. & Filipič, M. The role of reactive oxygen species in microcystin-LR-induced DNA damage. Toxicology 200, 59–68 (2004).PubMed 
Article 
CAS 

Google Scholar 
Odetti, L. M., López González, E. C., Romito, M. L., Simoniello, M. F. & Poletta, G. L. Genotoxicity and oxidative stress in Caiman latirostris hatchlings exposed to pesticide formulations and their mixtures during incubation period. Ecotoxicol. Environ. Saf. 193, 110312 (2020).CAS 
PubMed 
Article 

Google Scholar 
Rutkoski, C. F. et al. Morphological and biochemical traits and mortality in Physalaemus gracilis (Anura: Leptodactylidae) tadpoles exposed to the insecticide chlorpyrifos. Chemosphere 250, 126162 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Herek, J. S. et al. Genotoxic effects of glyphosate on Physalaemus tadpoles. Environ. Toxicol. Pharmacol. 81, 103516 (2021).CAS 
PubMed 
Article 

Google Scholar 
Natale, G. S. et al. Lethal and sublethal effects of the pirimicarb-based formulation Aficida® on Boana pulchella (Duméril and Bibron, 1841) tadpoles (Anura, Hylidae). Ecotoxicol. Environ. Saf. 147, 471–479 (2018)
Google Scholar 
Gilbert, S. F. Developmental Biology, 8th edn. (Sinauer Associates, 2006).Soto, M., García-Santisteban, I., Krenning, L., Medema, R. H. & Raaijmakers, J. A. Chromosomes trapped in micronuclei are liable to segregation errors. J. Cell Sci. 131, 214742 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Crott, J. & Fenech, M. Preliminary study of the genotoxic potential of homocysteine in human lymphocytes in vitro. Mutagenesis 16, 213–217 (2001).CAS 
PubMed 
Article 

Google Scholar 
Benvindo-Souza, M. et al. Micronucleus test in tadpole erythrocytes: Trends in studies and new paths. Chemosphere 240, 124910 (2020).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Fenech, M. The in vitro micronucleus technique. Mutat. Res. 455, 81–95 (2000).CAS 
PubMed 
Article 

Google Scholar 
Podratz, J. L. et al. Drosophila melanogaster: A new model to study cisplatin-induced neurotoxicity. Neurobiol. Dis. 43, 330–337 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Iturburu, F. G. et al. Uptake, distribution in different tissues, and genotoxicity of imidacloprid in the freshwater fish Australoheros facetus. Environ. Toxicol. Chem. 36, 699–708 (2017).CAS 
PubMed 
Article 

Google Scholar 
Vieira, C. E. D., Pérez, M. R., Acayaba, R. D. A., Raimundo, C. C. M. & Martinez, C. B. R. DNA damage and oxidative stress induced by imidacloprid exposure in different tissues of the Neotropical fish Prochilodus lineatus. Chemosphere 195, 125–134 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Sanchéz-Bayo, F., Goka, K. & Hayasaka, D. Contamination of the aquatic environment with neonicotinoids and its implication for ecosystems. Front. Environ. Sci. 4, 71 (2016).Article 

Google Scholar 
Wood, T. & Goulson, D. The environmental risks of neonicotinoid pesticides: a review of the evidence post-2013. Environ. Sci. Pollut. Res. 24, 17285–17325 (2017).CAS 
Article 

Google Scholar 
Craddock, H. A., Huang, D., Turner, P.C., Quirós-Alcalá, L. & Payne-Sturges, D. C. Trends in neonicotinoid pesticide residues in food and water in the United States, 1999–2015. Environ. Health 18, 7 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Heyer, R. et al. Leptodactylus latrans. IUCN Red List https://doi.org/10.2305/IUCN.UK.2010-2.RLTS.T57151A11592655.en (2010).Ade, C. M., Boone, M. D. & Puglis, H. J. Effects of an insecticide and potential predators on green frogs and northern cricket frogs. J. Herpetol. 44, 591–600 (2010).Article 

Google Scholar 
Sarkar, M. A., Roy, S., Kole, R. K. & Chowdhury, A. Persistence and metabolism of imidacloprid in different soils of West Bengal. Pest Manag. Sci. 57, 598–602 (2001).CAS 
PubMed 
Article 

Google Scholar 
Goulson, D. Review: An overview of the environmental risks posed by neonicotinoid insecticides. J. Appl. Ecol. 50, 977–987 (2013).Article 

Google Scholar 
Mineau, P. Neonic insecticides and invertebrate species endangerment. In Reference Module in Earth Systems and Environmental Sciences https://doi.org/10.1016/B978-0-12-821139-7.00126-4 (2021).Yamamuro, M. et al. Neonicotinoids disrupt aquatic food webs and decrease fishery yields. Science 366, 620–623 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Gosner. K. L. A simplified table for staging anuran embryos and larvae with notes on identification. Herpetologica 16, 183–190 (1960).
Google Scholar 
Percie-du-Sert, N. et al. The ARRIVE guidelines 2.0: Updated guidelines for reporting animal research. PLoS Biol. 18, e3000410 (2020). CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Herkovits, J. & Pérez-Coll, C. S. AMPHITOX: A customized set of toxicity tests employing amphibian embryos. Symposium on multiple stressor effects in relation to declining amphibian populations. In Multiple Stressor Effects in Relation to Declining Amphibian Populations (eds Linder, G. et al.) 46–60 (ASTM International STP 1443, 2003).Merga, L. B. & Van den Brink, P. J. Ecological effects of imidacloprid on a tropical freshwater ecosystem and subsequent recovery dynamics. Sci. Total Environ. 784, 147167 (2021).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Bonmatin, J.-M. et al. Environmental fate and exposure; neonicotinoids and fipronil. Environ. Sci. Pollut. Res. 22, 35–67 (2015).CAS 
Article 

Google Scholar 
Sumon, K. A. et al. Effects of imidacloprid on the ecology of sub-tropical freshwater microcosms. Environ. Pollut. 236, 432–441 (2018).CONCEA – Conselho Nacional de Controle e Experimentação Animal. Resolução normativa Nº 25, 29 de setembro de 2015. Guia Brasileiro de Produção, Manutenção ou Utilização de Animais para Atividades de Ensino ou Pesquisa Científica do Conselho Nacional de Controle e Experimentação Animal. http://www.mctic.gov.br/mctic/export/sites/institucional/institucional/concea/arquivos/legislacao/resolucoes_normativas/Resolucao-Normativa-CONCEA-n-27-de-23.10.2015-D.O.U.-de-27.10.2015-Secao-I-Pag.-10.pdf. (2015).Rutkoski, C. F. et al. Lethal and sublethal effects of the herbicide atrazine in the early stages of development of Physalaemus gracilis (Anura: Leptodactylidae). Arch. Environ. Contam. Toxicol. 74, 587–593 (2018).CAS 
PubMed 
Article 

Google Scholar 
Pérez-Iglesias, J. M., Soloneski, S., Nikoloff, N., Natale, G. S. & Larramendy, M. L. Toxic and genotoxic effects of the imazethapyr-based herbicide formulation Pivot H® on montevideo tree frog Hypsiboas pulchellus tadpoles (Anura, Hylidae). Ecotoxicol. Environ. Saf. 119, 15–24 (2015).PubMed 
Article 
CAS 

Google Scholar 
Montalvão, M. F. et al. The genotoxicity and cytotoxicity of tannery effluent in bullfrog (Rana catesbeianus). Chemosphere 183, 491–502 (2017).ADS 
PubMed 
Article 
CAS 

Google Scholar  More

Brandes, N. & Linial, M. Giant viruses—big surprises. Viruses 11, 404 (2019).CAS 
Article 

Google Scholar 
Jover, L. F., Effler, T. C., Buchan, A., Wilhelm, S. W. & Weitz, J. S. The elemental composition of virus particles: implications for marine biogeochemical cycles. Nat. Rev. Microbiol. 12, 519–528 (2014).CAS 
Article 

Google Scholar 
Madsen, E. L. Microorganisms and their roles in fundamental biogeochemical cycles. Curr. opinion biotechnology 22, 456–464 (2011).CAS 
Article 

Google Scholar 
Gummow, B. Challenges posed by new and re-emerging infectious diseases in livestock production, wildlife and humans. Livest. Sci. 130, 41–46 (2010).CAS 
Article 

Google Scholar 
Becker, D. J., Streicker, D. G. & Altizer, S. Linking anthropogenic resources to wildlife–pathogen dynamics: a review and meta-analysis. Ecol. letters 18, 483–495 (2015).Article 

Google Scholar 
Woolhouse, M. E. & Gowtage-Sequeria, S. Host range and emerging and reemerging pathogens. Emerg. infectious diseases 11, 1842 (2005).Article 

Google Scholar 
Foster, R. et al. Pathogens co-transported with invasive non-native aquatic species: implications for risk analysis and legislation. NeoBiota 69, 79–102 (2021).Article 

Google Scholar 
Brasier, C. The biosecurity threat to the uk and global environment from international trade in plants. Plant Pathol. 57, 792–808 (2008).Article 

Google Scholar 
Ruiz, G. M. et al. Global spread of microorganisms by ships. Nature 408, 49–50 (2000).ADS 
CAS 
Article 

Google Scholar 
Essl, F. et al. Which taxa are alien? criteria, applications, and uncertainties. BioScience 68, 496–509 (2018).Article 

Google Scholar 
Blackburn, T. M., Bellard, C. & Ricciardi, A. Alien versus native species as drivers of recent extinctions. Front. Ecol. Environ. 17, 203–207 (2019).Article 

Google Scholar 
Hawkins, C. L. et al. Framework and guidelines for implementing the proposed iucn environmental impact classification for alien taxa (eicat). Divers. Distributions 21, 1360–1363 (2015).Article 

Google Scholar 
Corrales, X. et al. Advances and challenges in modelling the impacts of invasive alien species on aquatic ecosystems. Biol. Invasions 22, 907–934 (2020).Article 

Google Scholar 
Regulation, E. Regulation (eu) no 1143/2014 of the European Parliament and of the Council of 22 October 2014 on the prevention and management of the introduction and spread of invasive alien species. Off. J. Eur. Union 57, 35–55 (2014).
Google Scholar 
EU. Regulation (eu) 2016/2031 of the European Parliament of the Council of 26 October 2016 on protective measures against pests of plants, amending regulations (eu) 228/2013,(eu) 652/2014 and (eu) 1143/2014 and repealing council directives 69/464/eec, 74/647/eec, 93/85/eec, 98/57/ec, 2000/29/ec, 2006/91/ec and 2007/33/ec. Off. J. 317, 4–104 (2016).
Google Scholar 
Murtaugh, M. P. et al. The science behind one health: at the interface of humans, animals, and the environment. Tech. Rep. (2017).Ogden, N. H. et al. Emerging infectious diseases and biological invasions: a call for a one health collaboration in science and management. Royal Soc. open science 6, 181577 (2019).ADS 
Article 

Google Scholar 
Roy, H. E. et al. Alien pathogens on the horizon: Opportunities for predicting their threat to wildlife. Conserv. Lett. 10, 477–484 (2017).Article 

Google Scholar 
Ikner, L. A., Gerba, C. P. & Bright, K. R. Concentration and recovery of viruses from water: a comprehensive review. Food Environ. Virol. 4, 41–67 (2012).
Google Scholar 
Taylor, M. W. Introduction: A short history of virology. In Viruses and Man: A History of Interactions, 1–22 (Springer, 2014).Thakur, M. P., Van der Putten, W. H., Cobben, M. M., van Kleunen, M. & Geisen, S. Microbial invasions in terrestrial ecosystems. Nat. Rev. Microbiol. 17, 621–631 (2019).CAS 
Article 

Google Scholar 
Desprez-Loustau, M.-L. et al. The fungal dimension of biological invasions. Trends ecology & evolution 22, 472–480 (2007).Article 

Google Scholar 
Rivett, D. W. et al. Elevated success of multispecies bacterial invasions impacts community composition during ecological succession. Ecol. Lett. 21, 516–524 (2018).Article 

Google Scholar 
Dunn, A. M. & Hatcher, M. J. Parasites and biological invasions: parallels, interactions, and control. TRENDS Parasitol. 31, 189–199 (2015).Article 

Google Scholar 
Pyšek, P. et al. Macroecological framework for invasive aliens (mafia): disentangling large-scale context dependence in biological invasions. (2020).Hulme, P. E. et al. Blurring alien introduction pathways risks losing the focus on invasive species policy. Conserv. Lett. 10, 265–266 (2017).Article 

Google Scholar 
Gilroy, J. J., Avery, J. D. & Lockwood, J. L. Seeking international agreement on what it means to be “native”. Conserv. Lett. 10, 238–247 (2017).Article 

Google Scholar 
Webber, B. L. & Scott, J. K. Rapid global change: implications for defining natives and aliens. Glob. Ecol. Biogeogr. 21, 305–311 (2012).Article 

Google Scholar 
CBD Secretariat. Decision VI/23: Alien species that threaten ecosystems, habitats and species. Document UNEP/CBD/COP/6/23 (2002).World Health Organization. A brief guide to emerging infectious diseases and zoonoses. Tech. Rep. https://apps.who.int/iris/handle/10665/204722 (2014).Firrao, G. et al. Candidatus phytoplasma’, a taxon for the wall-less, non-helical prokaryotes that colonize plant phloem and insects. Int. J. Syst. Evol. Microbiol. 54, 1243–1255 (2004).CAS 
Article 

Google Scholar 
CBD. Pathways of introduction of invasive species, their prioritization and management (Secretariat of the Convention on Biological Diversity Montreal, 2014).OIE. Terrestrial Animal Health Code 2021 (OIE, 2021).Magliozzi, C. et al. bacteria and viruses traits and species-related factors. figshare https://doi.org/10.6084/m9.figshare.18550907.v2 (2022).Katsanevakis, S. et al. Implementing the European policies for alien species: networking, science, and partnership in a complex environment. Manag. Biol. Invasions 4, 3–6 (2013).Article 

Google Scholar 
Tsiamis, K. et al. The EASIN Editorial Board: quality assurance, exchange and sharing of alien species information in europe. Manag. Biol. invasions 7, 321–328 (2016).Article 

Google Scholar 
Wieczorek, J. et al. Darwin core: an evolving community-developed biodiversity data standard. PloS one 7, e29715 (2012).ADS 
CAS 
Article 

Google Scholar 
Darwin Core. Darwin Core quick reference guide. https://dwc.tdwg.org/terms/ (2018).R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, https://www.R-project.org/ (2020).Wickham, H. ggplot2: Elegant Graphics for Data Analysis, https://ggplot2.tidyverse.org (Springer-Verlag New York, 2016).Schwarzl, T. ggBubbles: Mini Bubble Plots for Comparison of Discrete Data with ‘ggplot2’ R package version 0.1.4 (2019).Moon, K. R statistics and graphs for medical papers (Hannarae Seoul, 2015).Current, C. Invasive species compendium. Wallingford, UK: CAB Int. Available online: www.cabi.org/isc (accessed on 19 August 2020) (2011).Adams, M. J. & Antoniw, J. F. Dpvweb: An open access internet resource on plant viruses and virus diseases. Outlooks on Pest Manag. 16, 268 (2005).Article 

Google Scholar 
Adams, M. J. & Antoniw, J. F. Dpvweb: a comprehensive database of plant and fungal virus genes and genomes. Nucleic acids research 34, D382–D385 (2006).CAS 
Article 

Google Scholar 
Benson, D. A. et al. Genbank. Nucleic acids research 41, D36–D42 (2012).Article 

Google Scholar  More

Darwin, C. On the various Contrivances by which British and Foreign Orchids are Fertilised by Insects, and on the good effects of Intercrossing. (John Murray, 1862).Darwin, C. The various Contrivances by which Orchids are Fertilised by Insects. Second edition, revised., (D. Appleton and Company, 1877).Sprengel, C. K. Das entdeckte Geheimnis der Natur im Bau und in der Befruchtung der Blumen. (Vieweg, 1793).Müller, H. Befruchtung der Blumen durch Insekten (Verlag Von Wilhelm Englemann, 1873).Book 

Google Scholar 
Riley, C. V. The yucca moth and yucca pollination. Rep. Missouri Botan. Garden 3, 99–159 (1892).Article 

Google Scholar 
Faegri, K. & Van Der Pijl, L. Principles of Pollination Ecology 3rd edn. (Pergamon, Berlin, 1979).
Google Scholar 
Fenster, C. B., Armbruster, W. S., Wilson, P., Dudash, M. R. & Thomson, J. D. Pollination syndromes and floral specialization. Annu. Rev. Ecol. Evol. Syst. 35, 375–403. https://doi.org/10.1146/annurev.ecolsys.34.011802.132347 (2004).Article 

Google Scholar 
Inouye, D. W. In The Biology of Nectaries (eds Elias, T. S. & Bentley, B. L.) 153–173 (Columbia University Press, 1983).
Google Scholar 
Irwin, R. E., Bronstein, J. L., Manson, J. S. & Richardson, L. Nectar robbing: ecological and evolutionary perspectives. Annu. Rev. Ecol. Evol. Syst. 41, 271–292. https://doi.org/10.1146/annurev.ecolsys.110308.120330 (2010).Article 

Google Scholar 
Irwin, R. E. & Maloof, J. E. Variation in nectar robbing over time, space, and species. Oecologia 133, 525–533. https://doi.org/10.1007/s00442-002-1060-z (2002).ADS 
Article 
PubMed 

Google Scholar 
Maloof, J. E. & Inouye, D. W. Are nectar robbers cheaters or mutualists?. Ecology 81, 2651–2661. https://doi.org/10.1890/0012-9658(2000)081[2651:ANRCOM]2.0.CO;2 (2000).Article 

Google Scholar 
Inouye, D. W. The terminology of floral larceny. Ecology 61, 1251–1253. https://doi.org/10.2307/1936841 (1980).Article 

Google Scholar 
Lyon, D. L. & Chadek, C. Exploitation of nectar resources by hummingbirds, bees (Bombus), and Diglossa baritula and Its role in the evolution of Penstemon kunthii. Condor 73, 246–248. https://doi.org/10.2307/1365847 (1971).Article 

Google Scholar 
Colwell, R. K., Betts, B. J., Bunnell, P., Carpenter, F. L. & Feinsinger, P. Competition for the nectar of Centropogon valerii by the hummingbird Colibri thalassinus and the flower-piercer Diglossa plumbea, and Its evolutionary implications. Condor 76, 447–452. https://doi.org/10.2307/1365817 (1974).Article 

Google Scholar 
Arizmendi, M. C., Dominguez, C. A. & Dirzo, R. The role of an avian nectar robber and of hummingbird pollinators in the reproduction of two plant species. Funct. Ecol. 10, 119–127. https://doi.org/10.2307/2390270 (1996).Article 

Google Scholar 
Arizmendi, M. C. Multiple ecological interactions: Nectar robbers and hummingbirds in a highland forest in Mexico. Can. J. Zool. 79, 997–1006. https://doi.org/10.1139/z01-066 (2001).Article 

Google Scholar 
Navarro, L. Pollination ecology and effect of nectar removal in Macleania bullata (Ericaceae)1. Biotropica 31, 618–625. https://doi.org/10.1111/j.1744-7429.1999.tb00410.x (1999).Article 

Google Scholar 
Traveset, A., Willson, M. F. & Sabag, C. Effect of nectar-robbing birds on fruit set of Fuchsia magellanica in Tierra Del Fuego: A disrupted mutualism. Funct. Ecol. 12, 459–464. https://doi.org/10.1046/j.1365-2435.1998.00212.x (1998).Article 

Google Scholar 
Skutch, A. F. Life histories of Central American birds. Families Fringillidae, Thraupidae Parulidae and Coerebidae. Pacific Coast Avifauna 31, 1–448 (1954).
Google Scholar 
Vuilleumier, F. Systematics and evolution in Diglossa (Aves, Coerebidae). Am. Mus. Novit. 2381, 1–44 (1969).
Google Scholar 
Graves, G. R. Pollination of a Tristerix mistletoe (Loranthaceae) by Diglossa (Aves: Thraupidae). Biotropica 14, 315–317. https://doi.org/10.2307/2388094 (1982).Article 

Google Scholar 
Hernández, H. M. & Toledo, V. M. The role of nectar robbers and pollinators in the reproduction of Erythrina leptorhiza. Ann. Mo. Bot. Gard. 66, 512–520. https://doi.org/10.2307/2398843 (1979).Article 

Google Scholar 
Neill, D. A. Trapliners in the trees: Hummingbird pollination of Erythrina Sect Erythrina (Leguminosae: Papilionoideae). Ann. Missouri Botan. Garden 74, 27–41. https://doi.org/10.2307/2399259 (1987).Article 

Google Scholar 
Hazlehurst, J. A. & Karubian, J. O. Nectar robbing impacts pollinator behavior but not plant reproduction. Oikos 125, 1668–1676. https://doi.org/10.1111/oik.03195 (2016).CAS 
Article 

Google Scholar 
Cuta-Pineda, J. A., Arias-Sosa, L. A. & Pelayo, R. C. The flowerpiercers interactions with a community of high Andean plants. Avian Res. 12, 22. https://doi.org/10.1186/s40657-021-00256-7 (2021).Article 

Google Scholar 
Askins, R. A., Karen, M. E. & Jeffrey, D. W. Flower destruction and nectar depletion by avian nectar robbers on a tropical tree, Cordia sebestena. J. Field Ornithol. 58, 345–349 (1987).
Google Scholar 
McDade, L. A. & Kinsman, S. The impact of floral parasitism in two Neotropical hummingbird-pollinated plant species. Evolution 34, 944–958. https://doi.org/10.2307/2408000 (1980).Article 
PubMed 

Google Scholar 
Ingels, J. Observations of the hummingbirds Orthorhynchus cristatus and Eulampis jugularis of Martinique (West Indies). Gerfaut 66, 129–132 (1976).
Google Scholar 
Feinsinger, P., Beach, J. H., Linhart, Y. B., Busby, W. H. & Murray, K. G. Disturbance, pollinator predictability, and pollination success among Costa Rican cloud forest plants. Ecology 68, 1294–1305. https://doi.org/10.2307/1939214 (1987).Article 

Google Scholar 
Kodric-Brown, A., Brown, J. H., Byers, G. S. & Gori, D. F. Organization of a tropical island community of hummingbirds and flowers. Ecology 65, 1358–1368. https://doi.org/10.2307/1939116 (1984).Article 

Google Scholar 
Lara, C. & Ornelas, J. F. Preferential nectar robbing of flowers with long corollas: Experimental studies of two hummingbird species visiting three plant species. Oecologia 128, 263–273. https://doi.org/10.1007/s004420100640 (2001).ADS 
Article 
PubMed 

Google Scholar 
Hazlehurst, J. A. & Karubian, J. O. Impacts of nectar robbing on the foraging ecology of a territorial hummingbird. Behav. Proc. 149, 27–34. https://doi.org/10.1016/j.beproc.2018.01.001 (2018).Article 

Google Scholar 
Boehm, M. A. Biting the hand that feeds you: Wedge-billed hummingbird is a nectar robber of a sicklebill-adapted Andean bellflower. Acta Amazon. 48, 146–150. https://doi.org/10.1590/1809-4392201703932 (2018).Article 

Google Scholar 
Igić, B., Nguyen, I. & Fenberg, P. B. Nectar robbing in the trainbearers (Lesbia, Trochilidae). PeerJ 8, e9561. https://doi.org/10.7717/peerj.9561 (2020).Article 

Google Scholar 
Lunardi, V. D. O., Silva, É. E., Silva, S. T. A. & Lunardi, D. G. Handroanthus impetiginosus (Bignoniaceae) as an important floral resource for synanthropic birds in the Brazilian semiarid. Oecol. Austr. https://doi.org/10.4257/oeco.2019.2301.12 (2019).Article 

Google Scholar 
Almeida, J. M., Missagia, C. C. C. & Alves, M. A. S. Effects of the availability of floral resources and neighboring plants on nectar robbery in a specialized pollination system. Curr. Zool. https://doi.org/10.1093/cz/zoab083 (2021).Article 

Google Scholar 
Rodríguez-Rodríguez, M. C. & Valido, A. Opportunistic nectar-feeding birds are effective pollinators of bird-flowers from Canary Islands: experimental evidence from Isoplexis canariensis (Scrophulariaceae). Am. J. Bot. 95, 1408–1415. https://doi.org/10.3732/ajb.0800055 (2008).Article 
PubMed 

Google Scholar 
Lohmann, L. G. Untangling the phylogeny of neotropical lianas (Bignonieae, Bignoniaceae). Am. J. Bot. 93, 304–318. https://doi.org/10.3732/ajb.93.2.304 (2006).CAS 
Article 
PubMed 

Google Scholar 
Olmstead, R. G., Zjhra, M. L., Lohmann, L. G., Grose, S. O. & Eckert, A. J. A molecular phylogeny and classification of Bignoniaceae. Am. J. Bot. 96, 1731–1743. https://doi.org/10.3732/ajb.0900004 (2009).CAS 
Article 
PubMed 

Google Scholar 
Lohmann, L. G. & Taylor, C. M. A new generic classification of tribe Bignonieae (Bignoniaceae). Ann. Mo. Bot. Gard. 99, 348–489. https://doi.org/10.3417/2003187 (2014).Article 

Google Scholar 
Gentry, A. H. Coevolutionary patterns in Central American bignoniaceae. Ann. Mo. Bot. Gard. 61, 728–759. https://doi.org/10.2307/2395026 (1974).Article 

Google Scholar 
Bertin, R. I. Floral biology, hummingbird pollination and fruit production of trumpet creeper (Campsis radicans, Bignoniaceae). Am. J. Bot. 69, 122–134. https://doi.org/10.2307/2442837 (1982).Article 

Google Scholar 
Bertin, R. I. Paternity and fruit production in trumpet creeper (Campsis radicans). Am. Nat. 119, 694–709. https://doi.org/10.1086/283943 (1982).Article 

Google Scholar 
Bertin, R. I. & Sullivan, M. Pollen interference and cryptic self-fertility in Campsis radicans. Am. J. Bot. 75, 1140–1147. https://doi.org/10.1002/j.1537-2197.1988.tb08827.x (1988).Article 

Google Scholar 
Bertin, R. I. Paternal success following mixed pollinations of Campsis radicans. Am. Midl. Nat. 124, 153–163. https://doi.org/10.2307/2426088 (1990).Article 

Google Scholar 
Bertin, R. I. Effects of pollination intensity in Campsis radicans. Am. J. Bot. 77, 178–187. https://doi.org/10.1002/j.1537-2197.1990.tb13544.x (1990).Article 
PubMed 

Google Scholar 
Bertin, R. I. & Peters, P. J. Paternal effects on offspring quality in Campsis radicans. Am. Nat. 140, 166–178. https://doi.org/10.1086/285408 (1992).Article 

Google Scholar 
Kartesz, J. T. Campsis radicans. Floristic Synthesis of North America, Version 1.0. Biota of North America Program (BONAP) http://bonap.net/MapGallery/County/Campsis%20radicans.png. (2015).Kolodziejska-Degorska, I. & Zych, M. Bees substitute birds in pollination of ornitogamous climber Campsis radicans [L.] Seem in Poland. Acta Soc. Botanicorum Poloniae 75, 79–85 (2006).Article 

Google Scholar 
Catesby, M. The Natural History of Carolina, Florida and the Bahama islands. Volume 1. (Published by the author, 1731).Audubon, J. J. Ornithological Biography Vol. 3, 638 (Adam and Charles Black, 1835).
Google Scholar 
Audubon, J. J. Ruby-throated Hummingbird, plate CCLIII, The Birds of America Vol. 3 (Havell, 1835).
Google Scholar 
Nuttall, T. Manual of the Ornithology of the United States and of Canada. The Land Birds (Hilliard and Brown, 1832).
Google Scholar 
Stiles, F. G. & Freeman, C. E. Patterns in floral nectar characteristics of some bird-visited plant species from Costa Rica. Biotropica 25, 191–205. https://doi.org/10.2307/2389183 (1993).Article 

Google Scholar 
Stiles, F. G. Ecology, flowering phenology, and hummingbird pollination of some Costa Rican Heliconia species. Ecology 56, 285–301. https://doi.org/10.2307/1934961 (1975).Article 

Google Scholar 
McDade, L. A. & Weeks, J. A. Nectar in hummingbird-pollinated Neotropical plants I: Patterns of production and variability in 12 species. Biotropica 36, 196–215. https://doi.org/10.1111/j.1744-7429.2004.tb00312.x (2004).Article 

Google Scholar 
Wunderle, J. M. Jr. Nectar robbing by Orchard Orioles. Chat 44, 107–108 (1980).
Google Scholar 
Tyler, W. M. in Life histories of North American blackbirds, orioles, tanagers, and allies. Order Passeriformes: Families Ploceidae, Icteridae, and Thraupidae. United States National Museum Bulletin 211 (ed Arthur Cleveland Bent) 247–270 (United States Government Printing Office, 1958).George, F. W. Baltimore Orioles destroying trumpet vine blossoms. Wilson Bull. 46, 64 (1934).
Google Scholar 
Ridgway, R. The birds of North and Middle America, Part 2. Bull. U.S. Natl. Mus. 50, 1–834 (1902).
Google Scholar 
Scharf, W. C. & Kren, J. In Birds of the World (ed. Poole, A. F.) (Cornell Lab of Ornithology, 2020).
Google Scholar 
Morton, E. S. Effective pollination of Erythrina fusca by the Orchard Oriole (Icterus spurius): Coevolved behavioral manipulation?. Ann. Mo. Bot. Gard. 66, 482–489. https://doi.org/10.2307/2398840 (1979).Article 

Google Scholar 
Dickey, D. R. & van Rossem, A. J. The birds of El Salvador. Field Mus. Publ. Zool. 23, 1–609 (1938).
Google Scholar 
Baumel, J. J., King, A. S., Breazile, J. E., Evans, H. E. & Vanden Berge, J. C. (eds). Handbook of Avian Anatomy: Nomina Anatomica Avium, Second Edition. Publications of the Nuttall Ornithological Club no. 23 (Nuttall Ornithological Club, 1993).Beecher, W. J. Adaptations for food-getting in the American blackbirds. Auk 68, 411–440. https://doi.org/10.2307/4080840 (1951).Article 

Google Scholar 
Zusi, R. The role of the depressor mandibulae muscle in kinesis of the avian skull. Proc. U.S. Natl. Mus. 123, 1–28 (1967).Article 

Google Scholar 
Remsen, J. V. Jr. & Robinson, S. K. A classification scheme for foraging behavior of birds in terrestrial habitats. Stud. Avian Biol. 13, 144–160 (1990).
Google Scholar 
Skutch, A. F. Orioles, Blackbirds, and Their Kin (University of Arizona Press, 1996).
Google Scholar 
Hansell, M. P. Bird nests and Construction Behaviour 294 (Cambridge University Press, 2000).Book 

Google Scholar 
Bent, A. C. Life histories of North American blackbirds, orioles, tanagers, and allies. Bull. U.S. Natl. Museum 211, 1–531 (1958).
Google Scholar 
Dennis, J. V. Observations on the orchard oriole in lower Mississippi Delta. Bird-Banding 19, 12–21. https://doi.org/10.2307/4509997 (1948).Article 

Google Scholar 
Wunderle, J. M. & Lodge, D. J. The effect of age and visual cues on floral patch use by bananaquits (Aves: Emberizidae). Anim. Behav. 36, 44–54. https://doi.org/10.1016/S0003-3472(88)80248-3 (1988).Article 

Google Scholar 
Edge, A. A. Characteristics of nectar production and standing crop in Campsis radicans (Bignoniaceae). MSc thesis. (East Tennessee State University, 2010).Galetto, L. Nectary structure and nectar characteristics in some Bignoniaceae. Plant Syst. Evol. 196, 99–121. https://doi.org/10.1007/BF00985338 (1995).Article 

Google Scholar 
Elias, T. S. & Gelband, H. Nectar: Its production and functions in trumpet creeper. Science 189, 289–291. https://doi.org/10.1126/science.189.4199.289 (1975).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Elias, T. S. & Gelband, H. Morphology and anatomy of floral and extrafloral nectaries in Campsis (Bignoniaceae). Am. J. Bot. 63, 1349–1353. https://doi.org/10.1002/j.1537-2197.1976.tb13220.x (1976).Article 

Google Scholar 
Hermans, M. & Rasson, J. P. A new Sobolev test for uniformity on the circle. Biometrika 72, 698–702. https://doi.org/10.2307/2336748 (1985).MathSciNet 
Article 
MATH 

Google Scholar 
Landler, L., Ruxton, G. D. & Malkemper, E. P. The Hermans-Rasson test as a powerful alternative to the Rayleigh test for circular statistics in biology. BMC Ecol. 19, 30. https://doi.org/10.1186/s12898-019-0246-8 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
RStudio Team. RStudio: Integrated Development for R. PBC, Boston, MA http://www.rstudio.com/. (RStudio 2020).Beecher, W. J. Convergent evolution in the American orioles. Wilson Bulletin 62, 50–86 (1950).
Google Scholar 
Wolf, L. L., Hainsworth, F. R. & Stiles, F. G. Energetics of foraging: Rate and efficiency of nectar extraction by hummingbirds. Science 176, 1351–1352. https://doi.org/10.1126/science.176.4041.1351 (1972).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Wolf, L. L., Hainsworth, F. R. & Gill, F. B. Foraging efficiencies and time budgets in nectar-feeding birds. Ecology 56, 117–128. https://doi.org/10.2307/1935304 (1975).Article 

Google Scholar 
Alcantara, S. & Lohmann, L. G. Evolution of floral morphology and pollination system in Bignonieae (Bignoniaceae). Am. J. Bot. 97, 782–796. https://doi.org/10.3732/ajb.0900182 (2010).Article 
PubMed 

Google Scholar 
Gentry, A. H. Bignoniaceae: Part II (Tribe Tecomeae). Flora Neotrop. 25, 1–370 (1992).
Google Scholar 
Grant, V. Historical development of ornithophily in the western North American flora. Proc. Natl. Acad. Sci. 91, 10407–10411. https://doi.org/10.1073/pnas.91.22.10407 (1994).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
James, R. L. Some hummingbird flowers east of the Mississippi. Castanea 13, 97–109 (1948).
Google Scholar 
Van Nest, B. N., Edge, A. A., Feathers, M. V., Worley, A. C. & Moore, D. Bees provide pollination service to Campsis radicans (Bignoniaceae), a primarily ornithophilous trumpet flowering vine. Ecol. Entomol. 46, 117–127. https://doi.org/10.1111/een.12947 (2021).Article 

Google Scholar 
Patuxent Wildlife Research Center. Orchard oriole Icterus spurius. BBS summer distribution map, 2011–2015 (relative abundance map). https://www.mbr-pwrc.usgs.gov/bbs/ra2015/ra2015_red_v3.shtml (accessed 7 March 2021) (2021). More

Palaeogeographical modelWe use palaeogeographical reconstructions describing Earth’s palaeotopography and palaeobathymetry for a series of time slices from 541 Ma to the present day. The reconstructions merge existing models from two published global reconstruction datasets—those of ref. 32 and ref. 33 (https://doi.org/10.5281/zenodo.5348492), which themselves are syntheses of a wealth of previous work.For continental regions, estimates of palaeoelevation and continental flooding rely on a diverse range of geological evidence, such as sedimentary depositional environments and the spatiotemporal distribution of volcanic activity. For a full description, see a recent review34. Together, these data can be used to define the past locations of mountain ranges and palaeoshorelines34. For this part of our reconstruction, we used the compilation of ref. 33 with updated palaeoshorelines based on depositional environment information in current fossil databases35. This compilation comprises 82 palaeotopography maps covering the entire Phanerozoic. Note that each palaeogeographical map is a time slice representing the concatenation of geological data over several million years36.We quantified the impact of using the original compilation of ref. 33 on our model results and found only small changes with respect to using the reconstructions with updated palaeoshorelines (Extended Data Fig. 3a–c). Similarly, eustatic sea level is thought to have varied by around 100 m at timescales much shorter than the duration of the time-slices throughout the Phanerozoic37, such that the extent of continental flooding could have varied within each time slice by an amount significant for our analysis. For this reason, and to assess the uncertainty of our results to continental palaeogeography in general, we computed additional maps of continental flooding in which the sea level is raised or lowered by 100 m compared with the original palaeo–digital elevation model grids of ref. 33 (Extended Data Fig. 3d–f).For deep-ocean regions, the primary control on seafloor depth is the age of the seafloor, so reconstructing palaeobathymetry relies on constructing maps of seafloor age back in time38. As a consequence, we rely on reconstruction models that incorporate a continuous network of plate boundaries. For this part, we used the reconstruction of ref. 32 and derived maps of seafloor age from the plate tectonic model using the method of ref. 39 for which source code is available at GitHub (https://github.com/siwill22/agegrid-0.1). Palaeobathymetry is derived from the seafloor age maps following the steps outlined in ref. 38. It is important to note that seafloor age maps for most of the Phanerozoic (that is, pre-Pangaea times) are not directly constrained by data due to recycling of oceanic crust at subduction zones. Rather, they are model predictions generated by constructing plate motions and plate boundary configurations from the geological and palaeomagnetic record of the continents. Nonetheless, the first-order trends in ocean-basin volume and mean seafloor age are consistent with independent estimates for at least the last 410 million years (Myr)39.The reconstructions of refs. 32,33 differ in the precise locations of the continents through time. To resolve this discrepancy, we reverse reconstructed the continental palaeoelevation model of ref. 33 to present-day coordinates using their rotation parameters, then reconstructed them back in time using the rotations of ref. 32. Owing to the differences in how the continents are divided into different tectonic units, this process leads to some gaps and overlaps in the results40, which we resolved primarily through a combination of data interpolation and averaging. Manual adjustments were made to ensure that the flooding history remained consistent with the original palaeotopography in areas in which interpolation gives a noticeably different history of seafloor ages. The resulting palaeotopography maps are therefore defined in palaeomagnetic reference frame32 appropriate for use in Earth system models.For the biodiversity modelling, we generate estimates of the age of the seafloor for discrete points within the oceans and flooded continents, and track these ages through the lifetime of each point (Supplementary Video 5). For the oceans, this is achieved using the method described in ref. 39 in which the seafloor is represented by points that are incrementally generated at the mid-ocean ridges for a series of time steps 1 Myr apart, with each point tracked through subsequent time steps based on Euler poles of rotation until either the present-day is reached, or they arrive at a subduction zone and are considered to be destroyed.For the continents, tracking the location of discrete points is generally simpler as most crust is conserved throughout the timespan of the reconstruction. In contrast to the deep oceans (where we assume that crust is at all times submerged), we model the ‘age’ of the seafloor from the history of continental flooding and emergence within the palaeogeographical interpretation33. The continents are seeded with uniformly distributed points at the oldest timeslice (541 Ma) at which they are assigned an age of zero. These points are tracked to subsequent time slices of which the palaeogeography is used to determine whether the point lies within a flooded or emergent region. Points within flooded regions of continents are considered to be seafloor, and the age of this seafloor is accumulated across consecutive time slices where a given point lies within a flooded region. When a point is within an emergent region, the seafloor age is reset to zero. Following this approach, individual points within stable continents may undergo several cycles of seafloor age increasing from zero before being reset. At the continental margins formed during the Pangaea breakup, the age of the seafloor continuously grows from the onset of rifting. Intraoceanic island arcs represent an additional case, which can appear as new tectonic units with the reconstructions at various times. In these cases, we assume that the seafloor has a zero-age at the time at which the intraoceanic arc first develops, then remains predominantly underwater for the rest of its lifetime.Thus, for each of the 82 palaeogeographical reconstructions, we annotate 0.5° by 0.5° grids as continental, flooded continental shelf or oceanic for later use in model coupling and production of regional diversity maps.Palaeoenvironmental conditions under the cGENIE Earth system modelWe use cGENIE41, an Earth system model of intermediate complexity, to simulate palaeoenvironmental conditions of seawater temperature and organic carbon export production (as a surrogate for food supply) throughout the Phanerozoic (from 541 Ma to the present day).cGENIE is based on a three-dimensional (3D) ocean circulation model coupled to a 2D energy–moisture balance atmospheric component and a sea-ice module. We configured the model on a 36 × 36 (latitude, longitude) equal area grid with 17 unevenly spaced vertical levels in depth, down to a maximum ocean depth of 5,900 m. The cycling of carbon and associated tracers in the ocean is based on a size-structured plankton ecosystem model with a single (phosphate) nutrient42,43, and adopts an Arrhenius-type temperature-dependent scheme for the remineralization of organic matter exported to the ocean interior44.cGENIE provides a spatially resolved representation of ocean physics and biogeochemistry, which is a prerequisite for the present study to be able to reconstruct the spatial patterns of biodiversity in deep time. However, owing to the computational impracticality of generating a single transient simulation of physics (that is, temperature) and biogeochemistry (that is, export production) over the entire Phanerozoic, we therefore generate 30 model equilibria at regular time intervals throughout the Phanerozoic that are subsequently used as inputs for the regional diversification model (see the ‘Model coupling’ section).We used 30 Phanerozoic palaeogeographical reconstructions through time (~20 Myr evenly spaced time intervals) to represent key time periods. For each continental configuration corresponding to a given age in Earth history, we generate idealized 2D (but zonally averaged) wind speed and wind stress, and 1D zonally averaged albedo forcing fields45 required by the cGENIE model using the ‘muffingen’ open-source software (see the ‘Code availability’ section). For each palaeogeographical reconstruction, the climatic forcing (that is, solar irradiance and carbon dioxide concentration) is adapted to match the corresponding geological time interval. The partial pressure of CO2 is taken from the recent update of the GEOCARB model46. Solar luminosity is calculated using the model of stellar physics of ref. 47. We impose modern-day orbital parameters (obliquity, eccentricity and precession). The simulations are initialized with a sea-ice-free ocean, homogeneous oceanic temperature (5 °C) and salinity (34.9‰). As variations in the oceanic concentration of bio-available phosphate remain challenging to reconstruct in the geological past48,49, we impose a present-day mean ocean phosphate concentration (2.159 μmol kg−1) in our baseline simulations. We quantify the impact of this uncertainty on our model results by conducting additional simulations using half and twice the present-day ocean phosphate concentration (Extended Data Fig. 3g–i). For each ocean phosphate scenario (that is, 0.5×, 1× and 2× the present-day value), each of the 30 model simulations is then run for 20,000 years, a duration ensuring that deep-ocean temperature and geochemistry reach equilibrium. For each model simulation, the results of the mean annual values of the last simulated year are used for the analysis. Note that, although cGENIE makes projections of the distribution of dissolved oxygen ([O2]) in the ocean, our diversification model does not currently consider oxygenation to be a limit on diversity. Thus, we assumed a modern atmospheric partial pressure of O2 in all 30 palaeo simulations and did not use the resulting projected [O2] fields.Regional diversification modelWe tested two models of diversification—the logistic model and the exponential model—describing the dynamics of regional diversity over time. In both models, the net diversification rate (ρ), with units of inverse time (Myr−1), varies within a pre-fixed range of values as a function of seawater temperature and food availability. The net diversification rate is then calculated for a given location and time according to the following equation:$$rho ={rho }_{max }-({rho }_{max }-{rho }_{min })(1-{Q}_{{rm{temp}}}{Q}_{{rm{food}}})$$
(1)
where ρmin and ρmax set the lower and upper net diversification rate limits within which ρ is allowed to vary, and Qtemp and Qfood are non-dimensional limitation terms with values between 0 and 1 that define the dependence of ρ on temperature and food, respectively (Extended Data Table 1).The model considers a direct relationship between seawater temperature, food supply and the rate of net diversification on the basis of the theoretical control that temperature and food supply exert on the rates of origination and extinction (Supplementary Fig. 1). Temperature rise is expected to accelerate the biochemical kinetics of metabolism50 and shorten the development times of individuals51, leading to higher rates of mutation and origination. Greater food availability increases population sizes, which increases the rates of mutation and reduces the probability of extinction52. Furthermore, a large body of observations shows the existence of a positive relationship between resource availability (that is, food supply) and the standing stock of species in marine and terrestrial communities53,54. A larger food supply would support a greater number of individuals. A greater diversity of food resources could also lead to a finer partitioning of available resources55.The temperature dependence of ρ is calculated using the following equation:$${Q}_{{rm{temp}}}=frac{{Q}_{10}^{frac{T-Tmin }{10}}}{{Q}_{10}^{frac{Tmax -Tmin }{10}}}$$
(2)
where the Q10 coefficient measures the temperature sensitivity of the origination rate. In equation (2) above, T is the seawater temperature (in °C) at a given location and time, and Tmin and Tmax are the 0.01 percentile and the 0.99 percentile, respectively, of the temperature frequency distribution in each time interval. In the model, the values of Tmin and Tmax used to calculate Qtemp are therefore recomputed for every time interval (~5 Myr) according to the temperature frequency distribution of the corresponding time interval. This enables us to use updated Tmin and Tmax values in each Phanerozoic time interval and to account for the thermal adaptation of organisms to ever changing climate conditions.The food limitation term is parameterized using a Michaelis–Menten formulation as follows:$${Q}_{{rm{food}}}=frac{text{POC flux}}{left({K}_{{rm{food}}},+,text{POC flux}right)}$$
(3)
where POC flux (mol m−2 yr−1) is the particulate organic carbon export flux, which is used as a surrogate for food availability, at a given location and time of the simulated seafloor. The parameter Kfood (mol m−2 yr−1) in equation (3) is the half-saturation constant, that is, the POC flux at which the diversification rate is half its maximum value, provided that other factors were not limiting. These temperature and food supply limitation terms vary in space and time as a result of changes in seawater temperature and particulate organic carbon export rate, respectively, thereby controlling the spatial and temporal variability of ρ (Supplementary Video 6).The net diversification rate becomes negative (1) in the event of mass extinctions or (2) in response to regional-scale processes, such as sea-level fall and/or seafloor deformation along convergent plate boundaries. Mass extinction events are imposed as external perturbations to the diversification model by imputing negative net diversification rates to all active seafloor points (ocean points and flooded continental points) and assuming non-selective extinction. The percentage of diversity loss as well as the starting time and duration of mass extinctions are extracted from three fossil diversity curves of reference20,21,22 (Source Data for Fig. 1). Each of these fossil diversity curves provides different insights into the Phanerozoic history of marine animal diversity based on uncorrected range-through genus richness estimates20,22 and sampling standardized estimates21. Regional-scale processes—such as sea level fall during marine regressions and/or seafloor destruction at plate boundaries, either by subduction or uplift—are simulated by the combined plate tectonic–palaeoelevation model, and constrain the time that seafloor habitats have to accumulate diversity.The model assumes non-selective extinction during mass extinction events (that is, the field of bullets model of extinction; everything is equally likely to die, no matter the age of the clade and regardless of adaptation)56. However, there is much fossil evidence supporting extinction selectivity57,58. It could be argued that higher extinction rates at diversity hotspots would have delayed their subsequent recovery, flattening global diversity trends. This argument is difficult to reconcile with Sepkoski’s genus-level global diversity curve but could be consistent with the standardized diversity curve of ref. 21. Similarly, the model is also not suitable for reproducing the explosive radiations of certain taxonomic groups after mass extinctions, which could explain the offset between the model and fossil observations in the early Mesozoic (Fig. 1).Letting D represent regional diversity (number of genera within a given seafloor point) and t represent time, the logistic model is formalized by the following differential equation:$$frac{partial Dleft(tright)}{partial t}=rho Dleft[1-frac{D}{{K}_{{rm{eff}}}}right]$$
(4)
where D(t) is the number of genera at time t and Keff is the effective carrying capacity or maximum number of genera that a given seafloor point (that is, grid cell area after gridding) can carry at that time, t. In our logistic model, Keff is allowed to vary within a fixed range of values (from Kmin to Kmax) as a positive linear function of the POC flux at a given location and time as follows:$${K}_{{rm{eff}}}={K}_{max }-left({K}_{max }-{K}_{min }right)frac{{text{POC flux}}_{max }-text{POC flux}}{{text{POC flux}}_{max }-{text{POC flux}}_{min }}$$
(5)
where POC fluxmin and POC fluxmax correspond to the 0.01 and 0.99 quantiles of the POC flux range in the whole Phanerozoic dataset.In the logistic model, the net diversification rate decreases as regional diversity approaches its Keff. The exponential model is a particular case of the logistic model when Keff approaches infinity and, therefore, neither the origination rate nor the extinction rate depend on the standing diversities. In this scenario, diversity grows in an unlimited manner over time only truncated by the effect of mass extinctions and/or by the dynamics of the seafloor (creation versus destruction). Thus, the exponential model is as follows:$$frac{partial Dleft(tright)}{partial t}=rho D$$
(6)
where the rate of change of diversity (the time derivative) is proportional to the standing diversity D such that the regional diversity will follow an exponential increase in time at a speed controlled by the temperature- and food-dependent net diversification rate. Even if analytical solutions exist for the steady-state equilibrium of the logistic and exponential functions, we solved the ordinary differential equations (4) and (6) using numerical methods with a time lag of 1 Myr to account for the spatially and temporally varying environmental constraints, seafloor dynamics and mass extinction events.As the analysis of global fossil diversity curves is unable to discern the causes of diversity loss during mass extinctions, our imputation of negative diversification rates could have overestimated the loss of diversity in those cases in which sea level fall, a factor already accounted for by our model, contributed to mass extinction. This effect was particularly recognizable across the Permian–Triassic mass extinction (Extended Data Fig. 6d–f), and supports previous claims that the decline in the global area of the shallow water shelf exacerbated the severity of the end-Permian mass extinction34.Model couplingAs stated above, the coupled plate tectonic–palaeoelevation (palaeogeographical) model corresponds to a tracer-based model (a Lagrangian-based approach) that simulates and tracks the spatiotemporal dynamics of ocean and flooded continental points. The diversification models start at time 541 Ma with all active points having a D0 = 1 (one single genus everywhere) and we let points accumulate diversity heterogeneously with time according to seafloor age distributions (for ocean points) and the time that continents have been underwater (for flooded continental points). The ocean points are created at mid-ocean ridges and disappear primarily at subduction zones. Between their origin and demise, the points move following plate tectonic motions and we trace their positions while accumulating diversity. The flooded continental points begin to accumulate diversity from the moment that they are submerged, starting with a D value equal to the nearest neighbour flooded continental point with D  > 1, thereby simulating a process of coastal recolonization (or immigration). The diversification process remains active while the seafloor points remain underwater, but it is interrupted, and D set to 0, in those continental points that emerge above sea level. Similarly, seafloor points corresponding to ocean domains disappear in subduction zones, and their diversity is lost. We track the geographical position of the ocean and flooded continental points approximately every 5 Myr, from 541 Ma to the present. Each and every one of the tracked points accumulates diversity over time at a different rate, which is modulated by the environmental history (seawater temperature and food availability) of each point, as described in equations (1)–(3). When a point arrives in an environment with a carrying capacity lower than the diversity it has accumulated through time, we reset the diversity of the point to the value of the carrying capacity, thereby simulating local extinction.Seawater temperature (T) and food availability (POC flux) are provided by the cGENIE model, which has a spatial and temporal resolution coarser than the palaeogeographical model. The cGENIE model provides average seawater T and POC flux values in a 36 × 36 equal area grid (grid cell area equivalent to 2° latitude by 10° longitude at the equator) and 30 time slices or snapshots (from 541 Ma to the present: each ~20 Myr time intervals). To have environmental inputs for the 82 time slices of the plate tectonic–palaeoelevation model, we first interpolate the cGENIE original model output data on a 0.5° by 0.5° grid to match the annotated grids provided by the plate tectonic–palaeoelevation model. As the relatively coarse spatial resolution of the cGENIE model prevents rendering the coast–ocean gradients, we assign surface T and POC flux at the base of the euphotic zone to the flooded continental shelf grid cells, and deep ocean T and POC flux at the bottom of the ocean to the ocean grid cells. As there are time slices without input data of seawater T and POC flux, we interpolate/extrapolate seawater T and POC flux values into the 0.5° by 0.5° flooded continental shelf and ocean grids independently. Finally, we interpolate values from these 0.5° by 0.5° flooded continental shelf and ocean grids into the exact point locations in each time frame. Thus, each active point is tracked with its associated time-varying T and POC flux values throughout its lifetime. On average, 6,000 flooded continental points and 44,000 oceanic points were actively accumulating diversity in each time frame. The model cannot simulate the singularities of relatively small enclosed seas for which the spatial resolution of the palaeogeographical and Earth system models is insufficient to capture relevant features (such as palaeobathymetry, seawater temperature) in detail. The method is also likely to underestimate the diversity of epeiric (inland) seas due to the difficulty of simulating immigration, a process that is strongly influenced by the effect of surface ocean currents and is not considered here. However, as stated above, the model considers recolonization of recently submerged areas by marine biota from nearby coastal environments, which partially explains coastal immigration.Estimation of global diversity from regional diversityOur regional diversity maps are generated by separately interpolating ocean point diversity and flooded continental point diversity into the 0.5° by 0.5° annotated grids provided by the palaeogeographical model. We calculate global diversity at each time step from each of the regional diversity maps following a series of steps to integrate diversity along line transects from diversity peaks (maxima) to diversity troughs (minima) (Extended Data Fig. 1). To select the transects, first, we identify on each of the regional diversity maps the geographical position of the diversity peaks. We identify local maxima (that is, grid cells with diversity greater than their neighbour cells), and define the peaks as those local maxima with diversity greater than the 0.75 quantile of diversity values in all local maxima in the map. In the case of grid cells with equal neighbour diversity, the peak is assigned to the grid cell in the middle. We subsequently identify the geographical position of the diversity troughs, which are defined as newly formed ocean grid cells (age = 0 Myr) and, therefore, with diversities equal to one. The troughs are mostly located at mid-ocean ridges.On each of the 82 spatial diversity maps, we trace a line transect from each diversity peak to its closest trough, provided that the transect does not cross land in more than 20% of the grid cells along the linear path (Supplementary Video 7). On average, for each spatial diversity map, we trace 400 (σ = ±75) linear transects. This sampling design gives rise to transects of different lengths, which may bias the estimates of global diversity. To minimize this bias, we cut the tail of the transects to have a length of 555 km equivalent to 5° at the equator. We tested an alternative cut-off threshold, 1,110 km, and the results do not alter the study’s conclusions.We apply Bresenham’s line algorithm59 to detect the grid cells crossed by the transects and annotate their diversity. To integrate regional diversity along the transects, we developed a method to simplify the scenario of peaks and troughs heterogeneously distributed on the 2D diversity maps. The method requires (1) a vector (the transect) of genus richness (αn) at n different locations (grids) arranged in a line (1D) of L grids, and (2) a coefficient of similarity (Vn,n + 1) between each two neighbouring locations, n and n + 1. Vn,n + 1, the coefficient of similarity, follows a decreasing exponential function with distance between locations. The number of shared genera between n and n + 1 is Vn,n + 1 × min(αn; αn + 1). We integrate diversity from peaks to troughs and assume that, along the transect, αn + 1 is lower than αn. We further assume that the genera present in n and n + 2 cannot be absent from n + 1. Using this method, we integrate the transect’s diversity (γi) using the following equation:$${gamma }_{i}={ {mathbf{upalpha}}}_{1}+{sum }_{n=1}^{L-1}left(1-{V}_{n,n+1}right){ {mathbf{upalpha}}}_{n+1}$$
(7)
To integrate the diversity of all transects (γi) on each 2D diversity map (or time slice), we apply the same procedure as described above (Extended Data Fig. 1). We first sort the transects in descending order from the highest to the lowest diversity. We then assume that the number of shared genera between transect i and the rest of the transects with greater diversity {1, 2, …, i − 1} is given by the distance of its peak to the nearest neighbour peak (NN(i)) of those already integrated {1, 2, …, i − 1}. Thus, we perform a zigzag integration of transects’ diversities down gradient, from the greatest to the poorest, weighted by the nearest neighbour distance among the peaks already integrated. As a result, the contribution of each transect to global diversity will depend on its diversity and its distance to the closest transect out of all those transects already integrated. Using this method, we linearize the problem to simplify the cumbersome procedure of passing from a 2D regional diversity map to a global diversity estimate without knowing the identity (taxonomic affiliation) of the genera. If γtotal is the global diversity at time t:$${gamma }_{{rm{total}}}={gamma }_{1}+{sum }_{i=2}^{j}left(1-{V}_{{rm{NN}}left(iright),i}right){gamma }_{i}$$
(8)
Finally, the resulting global estimates are plotted against the midpoint value of the corresponding time interval to generate a synthetic global diversity curve. To compare the global diversity curves produced by the diversification models with those composed from the fossil record, Lin’s CCC60 is applied to the data normalized to the min–max values of each time series (that is, rescaled within the range 0–1). Lin’s CCC combines measures of both precision and accuracy to determine how far the observed data deviate from the line of perfect concordance (that is, the 1:1 line). Lin’s CCC increases in value as a function of the nearness of the data’s reduced major axis to the line of perfect concordance (the accuracy of the data) and of the tightness of the data around its reduced major axis (the precision of the data).The time series of global diversity generated from the fossil record and from the diversification model exhibit serial correlation and the resulting CCCs are therefore inflated. The use of methods for analysing non-zero autocorrelation time series data, such as first differencing or generalised least squares regression, enables high-frequency variations along the time series to be taken into account. However, the relative simplicity of our model, which was designed to reproduce the main Phanerozoic trends in global diversity, coupled with the fact that biases in the fossil data would introduce uncertainty into the analysis, leads us to focus our analysis on the long-term trends, obviating the effect of autocorrelation.Model parameterization and calibrationThe diversification models are parameterized assuming a range of values that constrain the lower and upper limits of the genus-level net diversification rate (ρmin and ρmax, respectively) (Extended Data Table 1) according to previously reported estimates from fossil records (figures 8 and 11 of ref. 5). A range of realistic values is assigned for the parameters Q10 and Kfood, determining, respectively, the thermal sensitivity and food dependence of the net diversification rate. We test a total of 40 different parameter combinations (Extended Data Table 2). The resulting estimates of diversity are then compared against the fossil diversity curves of ref. 20, ref. 21 or ref. 22, and the 15 parameter combinations providing the highest CCCs are selected.The results of the logistic diversification model rely on the values of the minimum and maximum carrying capacities (Kmin and Kmax, respectively) within which the spatially resolved effective carrying capacities (Keff) are allowed to vary. The values of Kmin and Kmax are therefore calibrated by running 28 simulations of pair-wise Kmin and Kmax combinations increasing in a geometric sequence of base 2, from 2 to 256 genera (Extended Data Fig. 4). We perform these simulations independently for each of the 15 parameter settings selected previously (Extended Data Fig. 4 and Extended Data Table 2). Each combination of Kmin and Kmax produces a global diversity curve, which is evaluated as described above using Lin’s CCC.Calculating estimates of global diversity from regional diversity maps in the absence of information on genus-level taxonomic identities requires that we assume a spatial turnover of taxa with geographical distance (Extended Data Fig. 1). Distance-decay curves are routinely fitted by calculating the ecological similarity (for example, the Jaccard similarity index) between each pair of sampling sites, and fitting an exponential decay function to the points on a scatter plot of similarity (y axis) versus distance (x axis). Following this method, we fit an exponential decay function to the distance–decay curves reported in ref. 61, depicting the decrease in the Jaccard similarity index (J) of fossil genera with geographical distance (great circle distance) at different Phanerozoic time intervals:$$J={J}_{{rm{o}}{rm{f}}{rm{f}}}+(,{J}_{max}text{-}{J}_{{rm{o}}{rm{f}}{rm{f}}}){{rm{e}}}^{-lambda times {rm{d}}{rm{i}}{rm{s}}{rm{t}}{rm{a}}{rm{n}}{rm{c}}{rm{e}}}$$
(9)
where Joff = 0.06 (n.d.) is a small offset, Jmax = 1.0 (n.d.) is the maximum value of the genus-based Jaccard similarity index and λ = 0.0024 (km−1) is the distance-decay rate.The Jaccard similarity index (J) between consecutive points n and n + 1 is bounded between 0 and min(αn; αn + 1)/max(αn; αn + 1). A larger value for J would mean that there are more shared genera between the two communities than there are genera within the least diverse community, which is ecologically absurd. However, using a single similarity decay function can lead the computed value of J to be locally larger than min(αn; αn + 1)/max(αn; αn + 1). To prevent this artefact, we use the Simpson similarity index or ‘overlap coefficient’ (V) instead of J. V corresponds to the percentage of shared genera with respect to the least diverse community (min(αn ; αn + 1)). V is bounded between 0 and 1, whatever the ratio of diversities. As the pre-existing estimates of similarity are expressed using J (ref. 61), we perform the conversion from J to V using the algebraic expression V = (1 + R) × J/(1 + J) where R = max(αn; αn + 1)/min(αn; αn + 1) (Supplementary Note 1). In the cases in which J exceeds the min(αn; αn + 1)/max(αn; αn + 1), V becomes >1 and, in those cases, we force V to be More