Ecology

Shiklomanov, I. A. & Rodda, J. C. World water resources at the beginning of the twenty-first century. (Cambridge University Press, 2003).Biggs, J., von Fumetti, S. & Kelly-Quinn, M. The importance of small waterbodies for biodiversity and ecosystem services: implications for policy makers. Hydrobiologia 793, 3–39 (2017).Article 

Google Scholar 
Heino, J. et al. Lakes in the era of global change: moving beyond single-lake thinking in maintaining biodiversity and ecosystem services. Biol. Rev. 96, 89–106 (2021).PubMed 
Article 

Google Scholar 
Janssen, A. B. G. et al. Shifting states, shifting services: linking regime shifts to changes in ecosystem services of shallow lakes. Freshw. Biol. 66, 1–12 (2021).Article 

Google Scholar 
Knoll, L. B. et al. Consequences of lake and river ice loss on cultural ecosystem services. Limnol. Oceanogr. Lett. 4, 119–131 (2019).Article 

Google Scholar 
Sterner, R. W. et al. Ecosystem services of Earth’s largest freshwater lakes. Ecosyst. Serv. 41, 101046 (2020).Article 

Google Scholar 
Reynaud, A. & Lanzanova, D. A global meta-analysis of the value of ecosystem services provided by lakes. Ecol. Econ. 137, 184–194 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Cooley, S. W., Ryan, J. C. & Smith, L. C. Human alteration of global surface water storage variability. Nature 591, 78–81 (2021).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Downing, J. A. Global limnology: up-scaling aquatic services and processes to planet Earth. SIL Proceedings, 1922–2010 30, 1149–1166 (2009).Article 

Google Scholar 
Tranvik, L. J., Cole, J. J. & Prairie, Y. T. The study of carbon in inland waters—from isolated ecosystems to players in the global carbon cycle. Limnol. Oceanogr. Lett. 3, 41–48 (2018).Article 

Google Scholar 
Balsamo, G. et al. On the contribution of lakes in predicting near-surface temperature in a global weather forecasting model. Tellus A Dyn. Meteorol. Oceanogr. 64, 15829 (2012).Article 

Google Scholar 
DelSontro, T., Beaulieu, J. J. & Downing, J. A. Greenhouse gas emissions from lakes and impoundments: upscaling in the face of global change. Limnol. Oceanogr. Lett. 3, 64–75 (2018).CAS 
Article 

Google Scholar 
Beaulieu, J. J. et al. Methane and carbon dioxide emissions from reservoirs: controls and upscaling. J. Geophys. Res. Biogeosciences 125, e2019JG005474 (2020).ADS 
CAS 
Article 

Google Scholar 
Slater, J. A. et al. The SRTM data “finishing” process and products. Photogramm. Eng. Remote Sens. 72, 237–247 (2006).Article 

Google Scholar 
Pekel, J.-F., Cottam, A., Gorelick, N. & Belward, A. S. High-resolution mapping of global surface water and its long-term changes. Nature 540, 418–422 (2016).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Verpoorter, C., Kutser, T., Seekell, D. A. & Tranvik, L. J. A global inventory of lakes based on high-resolution satellite imagery. Geophys. Res. Lett. 41, 6396–6402 (2014).ADS 
Article 

Google Scholar 
Pickens, A. H. et al. Mapping and sampling to characterize global inland water dynamics from 1999 to 2018 with full Landsat time-series. Remote Sens. Environ. 243, 111792 (2020).ADS 
Article 

Google Scholar 
Messager, M. L., Lehner, B., Grill, G., Nedeva, I. & Schmitt, O. Estimating the volume and age of water stored in global lakes using a geo-statistical approach. Nat. Commun. 7, 13603 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Tickner, D. et al. Bending the curve of global freshwater biodiversity loss: an emergency recovery plan. Bioscience 70, 330–342 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Downing, J. A., Polasky, S., Olmstead, S. M. & Newbold, S. C. Protecting local water quality has global benefits. Nat. Commun. 12, 1–6 (2021).Article 
CAS 

Google Scholar 
Hill, R. A., Weber, M. H., Debbout, R. M., Leibowitz, S. G. & Olsen, A. R. The Lake-Catchment (LakeCat) Dataset: characterizing landscape features for lake basins within the conterminous USA. Freshw. Sci. 37, 208–221 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Soranno, P. A. et al. LAGOS-NE: a multi-scaled geospatial and temporal database of lake ecological context and water quality for thousands of US lakes. Gigascience 6, 1–22 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Toptunova, O., Choulga, M. & Kurzeneva, E. Status and progress in global lake database developments. Adv. Sci. Res. 16, 57–61 (2019).Article 

Google Scholar 
Meyer, M. F., Labou, S. G., Cramer, A. N., Brousil, M. R. & Luff, B. T. The global lake area, climate, and population dataset. Sci. Data 7, 174 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Kling, G. W., Kipphut, G. W., Miller, M. M. & O’Brien, W. J. Integration of lakes and streams in a landscape perspective: the importance of material processing on spatial patterns and temporal coherence. Freshw. Biol. 43, 477–497 (2000).Article 

Google Scholar 
Fergus, C. E. et al. The freshwater landscape: lake, wetland, and stream abundance and connectivity at macroscales. Ecosphere 8, e01911 (2017).Article 

Google Scholar 
Lehner, B., Messager, ML., Korver, MC. & Linke, S. LakeATLAS Version 1.0, figshare, https://doi.org/10.6084/m9.figshare.19312001 (2022).Linke, S. et al. Global hydro-environmental sub-basin and river reach characteristics at high spatial resolution. Sci. data 6, 283 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Fergus, C. E. et al. National framework for ranking lakes by potential for anthropogenic hydro-alteration. Ecol. Indic. 122, 107241 (2021).Article 

Google Scholar 
Bracht-Flyr, B., Istanbulluoglu, E. & Fritz, S. A hydro-climatological lake classification model and its evaluation using global data. J. Hydrol. 486, 376–383 (2013).ADS 
Article 

Google Scholar 
Soranno, P. A. et al. Using landscape limnology to classify freshwater ecosystems for multi-ecosystem management and conservation. Bioscience 60, 440–454 (2010).Article 

Google Scholar 
McCullough, I. M., Skaff, N. K., Soranno, P. A. & Cheruvelil, K. S. No lake left behind: how well do U.S. protected areas meet lake conservation targets? Limnol. Oceanogr. Lett. 4, 183–192 (2019).Article 

Google Scholar 
Stanley, E. H. et al. Biases in lake water quality sampling and implications for macroscale research. Limnol. Oceanogr. 64, 1572–1585 (2019).ADS 
CAS 
Article 

Google Scholar 
Hanson, P. C., Weathers, K. C. & Kratz, T. K. Networked lake science: how the Global Lake Ecological Observatory Network (GLEON) works to understand, predict, and communicate lake ecosystem response to global change. Inl. Waters 6, 543–554 (2016).Article 

Google Scholar 
Lottig, N. R. & Carpenter, S. R. Interpolating and forecasting lake characteristics using long-term monitoring data. Limnol. Oceanogr. 57, 1113–1125 (2012).ADS 
Article 

Google Scholar 
Filazzola, A. et al. A database of chlorophyll and water chemistry in freshwater lakes. Sci. Data 2020 71 7, 1–10 (2020).
Google Scholar 
Lehner, B. & Messager, M. L. HydroLAKES – Technical Documentation Version 1.0. https://data.hydrosheds.org/file/technical-documentation/HydroLAKES_TechDoc_v10.pdf (2016).Natural Resources Canada. CanVec Hydrography: Waterbody Features. Version 12.0. https://ftp.maps.canada.ca/pub/nrcan_rncan/vector/canvec (2013).Lehner, B., Verdin, K. & Jarvis, A. New global hydrography derived from spaceborne elevation data. Eos, Trans. AGU 89, 93–94 (2008).ADS 
Article 

Google Scholar 
Farr, T. G. & Kobrick, M. Shuttle radar topography mission produces a wealth of data. Eos, Trans. AGU 81, 583–585 (2000).ADS 
Article 

Google Scholar 
Müller Schmied, H. et al. The global water resources and use model WaterGAP v2.2d: model description and evaluation. Geosci. Model Dev. 14, 1037–1079 (2021).ADS 
Article 

Google Scholar 
Beck, H. E. et al. Global evaluation of runoff from 10 state-of-the-art hydrological models. Hydrol. Earth Syst. Sci. 21, 2881–2903 (2017).ADS 
Article 

Google Scholar 
Alcamo, J. et al. Development and testing of the WaterGAP 2 global model of water use and availability. Hydrol. Sci. J. 48, 317–338 (2003).Article 

Google Scholar 
Döll, P., Kaspar, F. & Lehner, B. A global hydrological model for deriving water availability indicators: model tuning and validation. J. Hydrol. 270, 105–134 (2003).ADS 
Article 

Google Scholar 
Lehner, B. & Grill, G. Global river hydrography and network routing: baseline data and new approaches to study the world’s large river systems. Hydrol. Process. 27, 2171–2186 (2013).ADS 
Article 

Google Scholar 
Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
Hengl, T. et al. SoilGrids250m: Global gridded soil information based on machine learning. PLoS One 12, e0169748 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Zhang, X. et al. GLC_FCS30: Global land-cover product with fine classification system at 30 m using time-series Landsat imagery. Earth Syst. Sci. Data 13, 2753–2776 (2021).ADS 
Article 

Google Scholar 
Buchhorn, M. et al. Copernicus Global Land Service: Land Cover 100m: Collection 3: epoch 2019: Globe, Zenodo, https://doi.org/10.5281/zenodo.3939050 (2020).ESRI. ArcGIS Desktop: Release 10.4.1 (Environmental Systems Research Institute, Redlands, CA, USA, 2016).Soranno, P. A., Cheruvelil, K. S., Wagner, T., Webster, K. E. & Bremigan, M. T. Effects of land use on lake nutrients: the importance of scale, hydrologic connectivity, and region. PLoS One 10, e0135454 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Su, Z. H., Lin, C., Ma, R. H., Luo, J. H. & Liang, Q. O. Effect of land use change on lake water quality in different buffer zones. Appl. Ecol. Environ. Res. 13, 639–653 (2015).
Google Scholar 
Brakebill, J. W., Schwarz, G. E. & Wieczorek, M. E. An enhanced hydrologic stream network based on the NHDPlus medium resolution dataset. Scientific Investigations Report https://doi.org/10.3133/sir20195127 (2020).Carroll, M., Townshend, J., DiMiceli, C., Noojipady, P. & Sohlberg, R. Global raster water mask at 250 meter spatial resolution, Collection 5: MOD44W MODIS Water Mask. College Park, Maryland: University of Maryland (2009).Carroll, M. L., Townshend, J. R., DiMiceli, C. M., Noojipady, P. & Sohlberg, R. A. A new global raster water mask at 250 m resolution. Int. J. Digit. Earth 2, 291–308 (2009).ADS 
Article 

Google Scholar 
European Environment Agency (EEA). European Catchments and Rivers Network System (ECRINS), https://www.eea.europa.eu/data-and-maps/data/european-catchments-and-rivers-network (2012).Ouellet Dallaire, C., Lehner, B., Sayre, R. & Thieme, M. A multidisciplinary framework to derive global river reach classifications at high spatial resolution. Environ. Res. Lett. 14, 024003 (2019).ADS 
Article 

Google Scholar 
Global Runoff Data Centre (GRDC). River discharge data. Federal Institute of Hydrology, 56068 Koblenz, Germany, https://www.bafg.de/GRDC (2014).Openshaw, S. The modifiable areal unit problem. In Quantitative Geography: A British View (eds. Wrigley, N. & Bennett, R.) 60–69 (Routledge and Kegan Paul, Andover, 1981).United States Census Bureau. 2010 Census. ftp://ftp2.census.gov/geo/tiger (2010).Center for International Earth Science Information Network (CIESIN) & NASA Socioeconomic Data and Applications Center (SEDAC). Gridded Population of the World, Version 4 (GPWv4): Population Count and Density. https://doi.org/10.7927/H4JW8BX5 (2016).Grill, G. et al. Mapping the world’s free-flowing rivers. Nature 569, 215–221 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Allen, D. J. et al. The Diversity of Life in African Freshwaters: Under Water, Under Threat: an Analysis of the Status and Distribution of Freshwater Species Throughout Mainland Africa. (IUCN, 2011).Markovic, D. et al. Europe’s freshwater biodiversity under climate change: distribution shifts and conservation needs. Divers. Distrib. 20, 1097–1107 (2014).Article 

Google Scholar 
Fluet-Chouinard, E., Lehner, B., Rebelo, L.-M., Papa, F. & Hamilton, S. K. Development of a global inundation map at high spatial resolution from topographic downscaling of coarse-scale remote sensing data. Remote Sens. Environ. 158, 348–361 (2015).ADS 
Article 

Google Scholar 
Lehner, B. et al. High‐resolution mapping of the world’s reservoirs and dams for sustainable river‐flow management. Front. Ecol. Environ. 9, 494–502 (2011).Article 

Google Scholar 
Fan, Y., Li, H. & Miguez-Macho, G. Global patterns of groundwater table depth. Science 339, 940–943 (2013).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Robinson, N., Regetz, J. & Guralnick, R. P. EarthEnv-DEM90: A nearly-global, void-free, multi-scale smoothed, 90m digital elevation model from fused ASTER and SRTM data. ISPRS J. Photogramm. Remote Sens. 87, 57–67 (2014).ADS 
Article 

Google Scholar 
Metzger, M. J. et al. A high-resolution bioclimate map of the world: a unifying framework for global biodiversity research and monitoring. Glob. Ecol. Biogeogr. 22, 630–638 (2013).Article 

Google Scholar 
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).Article 

Google Scholar 
Zomer, R. J., Trabucco, A., Bossio, D. A. & Verchot, L. V. Climate change mitigation: a spatial analysis of global land suitability for clean development mechanism afforestation and reforestation. Agric. Ecosyst. Environ. 126, 67–80 (2008).Article 

Google Scholar 
Trabucco, A., Zomer, R. J., Bossio, D. A., van Straaten, O. & Verchot, L. V. Climate change mitigation through afforestation/reforestation: a global analysis of hydrologic impacts with four case studies. Agric. Ecosyst. Environ. 126, 81–97 (2008).Article 

Google Scholar 
Trabucco, A. & Zomer, R. J. Global soil water balance geospatial database. CGIAR Consortium for Spatial Information, https://cgiarcsi.community/data/global-high-resolution-soil-water-balance (2010).Hall, D. K., Riggs, G. A. & Salomonson, V. MODIS/Terra snow cover daily L3 global 500m grid, version 5, 2002–2015, https://doi.org/10.5067/MODIS/MOD10A1.006 (2016).Bartholomé, E. & Belward, A. S. GLC2000: a new approach to global land cover mapping from Earth observation data. Int. J. Remote Sens. 26, 1959–1977 (2005).Article 

Google Scholar 
Ramankutty, N. & Foley, J. A. Estimating historical changes in global land cover: Croplands from 1700 to 1992. Global Biogeochem. Cycles 13, 997–1027 (1999).ADS 
CAS 
Article 

Google Scholar 
Lehner, B. & Döll, P. Development and validation of a global database of lakes, reservoirs and wetlands. J. Hydrol. 296, 1–22 (2004).ADS 
Article 

Google Scholar 
Ramankutty, N., Evan, A. T., Monfreda, C. & Foley, J. A. Farming the planet: 1. Geographic distribution of global agricultural lands in the year 2000. Global Biogeochem. Cycles 22, (2008).Siebert, S. et al. A global data set of the extent of irrigated land from 1900 to 2005. Hydrol. Earth Syst. Sci. 19, 1521–1545 (2015).ADS 
Article 

Google Scholar 
GLIMS & NSIDC. Global land ice measurements from space (GLIMS) glacier database, v1. National Snow and Ice Data Center (NSIDC), https://doi.org/10.7265/N5V98602 (2012).Gruber, S. Derivation and analysis of a high-resolution estimate of global permafrost zonation. Cryosphere 6, 221–233 (2012).ADS 
Article 

Google Scholar 
UNEP-WCMC & IUCN. The World Database on Protected Areas, http://www.protectedplanet.net (2014).Dinerstein, E. et al. An ecoregion-based approach to protecting half the terrestrial realm. Bioscience 67, 534–545 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Abell, R. et al. Freshwater ecoregions of the world: a new map of biogeographic units for freshwater biodiversity conservation. Bioscience 58, 403–414 (2008).Article 

Google Scholar 
Hengl, T. et al. SoilGrids1km—global soil information based on automated mapping. PLoS One 9, e105992 (2014).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hartmann, J. & Moosdorf, N. The new global lithological map database GLiM: a representation of rock properties at the Earth surface. Geochem. Geophys. Geosyst. 13, Q12004 (2012).ADS 
Article 

Google Scholar 
Williams, P. W. & Ford, D. C. Global distribution of carbonate rocks. Zeitschrift für Geomorphologie Suppl. 147, 1–2 (2006).
Google Scholar 
Borrelli, P. et al. An assessment of the global impact of 21st century land use change on soil erosion. Nat. Commun. 8, 1–13 (2017).ADS 
CAS 
Article 

Google Scholar 
Pesaresi, M. & Freire, S. GHS Settlement grid following the REGIO model 2014 in application to GHSL Landsat and CIESIN GPW v4-multitemporal (1975-1990-2000-2015). European Commission, Joint Research Centre (JRC), https://data.europa.eu/data/datasets/jrc-ghsl-ghs_smod_pop_globe_r2016a (2016).Doll, C. N. H. CIESIN thematic guide to night-time light remote sensing and its applications. CIESIN http://sedac.ciesin.columbia.edu/binaries/web/sedac/thematic-guides/ciesin_nl_tg.pdf (2008).Meijer, J. R., Huijbregts, M. A. J., Schotten, K. C. G. J. & Schipper, A. M. Global patterns of current and future road infrastructure. Environ. Res. Lett. 13, 64006 (2018).Article 

Google Scholar 
Venter, O. et al. Global terrestrial Human Footprint maps for 1993 and 2009. Sci. data 3, 160067 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
University of Berkeley. Database of global administrative areas (GADM). University of Berkeley, Museum of Vertebrate Zoology and the International Rice Research Institute, http://www.gadm.org (2012).Kummu, M., Taka, M. & Guillaume, J. H. A. Gridded global datasets for gross domestic product and Human Development Index over 1990–2015. Sci. data 5, 180004 (2018).PubMed 
PubMed Central 
Article 

Google Scholar  More

Study speciesCane toads (Rhinella marina) are large (to  > 1 kg) bufonids (Fig. 1a). Although native to north-eastern South America, these toads have been translocated to many countries worldwide to control insect pests12. Adult cane toads forage at night for insect prey and retreat to moist shelter-sites per day13. Small body size (and thus, high desiccation rate) restricts young toads to the margins of natal ponds14, but adult toads can survive even in highly arid habitats if they have access to water13,15. Cane toads prefer open habitats for foraging12, and thus can thrive in post-fire landscapes16,17. Cane toads in post-fire landscapes tend to have lower parasite burdens, probably because free-living larvae of their lungworm parasites cannot survive either the fire or the more sun-exposed post-fire landscape18.Figure 1taken from study sites between Casino, Grafton, and surrounds, NSW, by S.W. Kaiser.The cane toad Rhinella marina (a), and unburned, (b) and burned (c) habitats in which toads were collected and radio-tracked. Photographs were Full size imageStudy areaEast of the Great Dividing Range, near-coastal Clarence Dry Sclerophyll Forests of north-eastern New South Wales (NSW) are dominated by Spotted gum (Corymbia variegata) and Pink bloodwood (Corymbia intermedia)19. Fires are common, but typically cover relatively small areas before they are extinguished. In the summer of 2019–2020, however, prolonged drought followed by an unusually hot summer resulted in massive fires across this region, burning almost 100,000 km2 of vegetation9. In the current study, the toads we measured and dissected came from several sites within 75 km of the city of Casino (for site locations, see Fig. 2, Table 1, and18). The impacts of fire on faunal abundance and attributes shift with time since fire; for example, the abundance of a particular species may be reduced by fire (due to mortality from flames) but then increase as individuals from surrounding areas migrate to the recently-burned site to exploit new ecological opportunities provided by that landscape8. We chose to study this system 1-year post-fire, to allow time for such longer-term effects to be manifested.Figure 2Sampling sites relative to fire history. Sample sites are burned (red circles), and unburned (green squares). See Table 1 for key to sites. The legend shows the extent of burn a year prior to our study. Map created in QGIS 3.22.3. Fire history available from https://datasets.seed.nsw.gov.au/dataset/fire-extent-and-severity-mapping-fesm CC BY 4.0.Full size imageTable 1 Sampling sites and sample sizes for dissected and radio-tracked cane toads (Rhinella marina) in New South Wales, Australia.Full size tableSurveys of toad abundanceTo quantify toad abundance in burned and unburned sites, one observer (MJG) walked 100-m transects along roads at night (N = 23 and 8 respectively), recording all toads and native frogs (both adult and juvenile). The smaller number of unburned sites reflects the massive spatial scale of the wildfires, which made it difficult to find unburned areas. The transect sites were not the same as those sampled by “toad-busters” (below). We sampled both burned and unburned sites on each night, to de-confound effects of weather conditions with fire treatment. We scored frogs as well as toads to provide an estimate of overall anuran abundance and activity, and so that we could examine toad abundance relative to frog abundance as well as absolute toad numbers.“Toad-buster” sampleBecause of their ecological impact on native fauna, cane toads are culled by community groups as well as by government authorities12,20. We asked “toad-buster” groups to record whether the sites at which they collected toads had been burned during the 2019–2020 fires, or had remained unburned (Table 1). The toads were humanely euthanized (cooled-then-pithed: see21). The euthanasia method is brief (a few hours in the refrigerator, followed by pithing) and thus should not have affected any of the traits that we measured. For all of these toads, we measured body length (snout-urostyle length = SUL) and mass, and determined sex based on external morphology (skin colour and rugosity, nuptial pads: see22). A subset of toads (chosen to provide relatively equal numbers of males and females, and with equal numbers from burned and unburned sites) was dissected to provide data on mass of internal organs (fat bodies, liver, ovaries), reproductive condition (state of ovarian follicle development) and diet (mass and identity of prey items). To select the subsample of toads for dissection, we took relatively equal numbers of male and female toads from each bag of toads that was provided to us by the “toad-busters”. For logistical reasons, we were unable to dissect all of the toads that had been collected. Overall, we obtained data on morphology, diets and other traits from 481 fully dissected and 1443 partially dissected cane toads.Radio-trackingTo explore habitat use and movement patterns, we radio-tracked 57 toads over the course of two fieldtrips (0900–1800 h from 20 Nov 2021 to 6 Dec 2021 and 25 Jan 2022 to 10 Feb 2022). We selected seven sites (4 burned, 3 unburned) within 28 km of Tabbimoble, NSW (see Table 1 for locations and sample sizes of tracked toads). We hand-captured toads found active at night. These were measured, and their sex determined by external morphology (see above) and behaviour (release calls, given only by males: see23). We then fitted the toads with radio-transmitters (PD-2; Holohil Systems, Ontario, Canada; weighing ≤ 3.8 g) on cotton waist-belts, and released them at the site of capture. Tracked toads were 88.2–160.9 mm SUL (mass 70.1–546.3 g); thus, transmitters weighed  20 mm thick) within the quadrat, and estimated exposure of the toad within its refuge (the percentage of the animal’s body exposed to the naked eye). We then selected a compass bearing at random and walked 20 m in that direction where we rescored all of the above habitat attributes, to quantify habitat features in the broader environment (i.e., not just in microhabitats used by toads). We used those “random” sites to quantify overall habitat attributes of burned and unburned sites. Temperature was recorded by directing a temperature gun (Digitech QM7221) on (or otherwise close-to) toads and at a random point on the ground for random replicates. In total, we gathered radio-tracking data on movements and habitat variables from 57 cane toads, each of which was tracked for 5 days. Recaptured toads were euthanized by cooling-then-pithing.Morphological traitsTo obtain an index of body condition of toads, we regressed ln mass against ln SUL, and used the residual scores from that general linear regression as our estimate of body condition. Negative residual scores show an individual that weighs less-than-expected based on its body length. Likewise, we regressed mass of the fat bodies, liver and stomach against body mass to obtain indices of energy stores and stomach-content volumes relative to body mass. We scored male secondary sexual characteristics using the system of Bowcock et al.22. In their system, three sexually dimorphic traits (nuptial pad size, skin roughness and skin colouration) are scored from 0 to 2, and the scores from those three traits are summed to create a final value (on a 6-point scale) for the degree of elaboration of male-specific secondary sexual characteristics. We scored reproductive condition in adult female toads based on whether or not egg masses were visible during dissection, based on dissected toads from both “toad-buster” and telemetry samples.Statistical methodsData were analysed in R version 4.2.025. We used Linear Mixed Models (LMMs), Generalised Linear Mixed Models (GLMMs) and logistic regressions for our analyses. The R packages ‘tidyverse’26, ‘lmerTest’27, and ‘performance’28 were used.Habitat dataWe compared habitat variables between burned and unburned sites, and attributes of toads in burned versus unburned sites, using GLMMs (with negative binomial distribution) for count data (models were checked for overdispersion29) and LMMs on distance data, using ln-transformations where required to achieve normality. LMMs were used on non-normal percentage data, which were ln- and then logit-transformed (using log[(P + e)/(1 − P + e)], where e is the lowest non-zero number, halved)30. We used toad id, site (sampling location) and sampling trip (2019 versus 2020) as random factors.Anuran transect dataCounts of toads in burned versus unburned areas were compared both directly via GLMMs with a negative binomial distribution and relative to the numbers of frogs sighted along the same transects (binding the columns in R as ‘number of toads, number of amphibians – number of toads’ and using a GLMM with a binomial distribution). We used site as a random factor.Telemetry dataFor telemetry data, we analysed response variables via LMMs, and ln-transformed data where relevant to achieve normality.Dissection dataWe used LMMs for SUL, body mass, body condition and organ mass residuals (e.g., fat body mass relative to body mass). For prey item data, we used a poisson distribution with row number as a random factor, as the negative binomial and beta distribution GLMMs were overdispersed (see31). We used LMM for number of prey items and number of prey groups, with site as a random factor. Where models failed to converge, we reduced or removed the error term(s). Analyses were restricted to toads ≥ 70 mm SUL, because animals below this size were difficult to sex. We also performed nominal logistic regression to explore variation in sex ratio and male secondary sexual traits.Reproductive conditionWe used LMM for male secondary sexual characteristic display, using site as a random factor. For ovary presence, we used a binomial GLMM with a logit link, using site as a random factor. We used a LMM of the residual values from ovary mass relative to body mass (ln-transformed), using site as a random factor.Ethics declarationsAll procedures were performed in accordance with the relevant guidelines and regulations approved by Macquarie University Animal Ethics Committee (ARA Number: 2019/040-2) and in accordance with ARRIVE guidelines. More

Thiel, M. et al. The Humboldt Current System of northern and central Chile—Oceanographic processes, ecological interactions and socioeconomic feedback. Oceanogr. Mar. Biol. Annu. Rev. 45, 195–344 (2007).
Google Scholar 
FAO. The State of World Fisheries and Aquaculture 2020. Sustainability in action 2020.Castilla, J. C. & Camus, P. A. The Humboldt-El Niño scenario: Coastal benthic resources and anthropogenic influences, with particular reference to the 1982/83 ENSO. S. Afr. J. Mar. Sci. 12, 703–712. https://doi.org/10.2989/02577619209504735 (1992).Article 

Google Scholar 
Alheit, J. & Niquen, M. Regime shifts in the Humboldt Current ecosystem. Prog. Oceanogr. 60, 201–222. https://doi.org/10.1016/j.pocean.2004.02.006 (2004).Article 

Google Scholar 
González, H. E. et al. Carbon fluxes within the epipelagic zone of the Humboldt Current System off Chile: The significance of euphausiids and diatoms as key functional groups for the biological pump. Prog. Oceanogr. 83, 217–227. https://doi.org/10.1016/j.pocean.2009.07.036 (2009).Article 

Google Scholar 
Quiñones, R. A., Levipan, H. A. & Urrutia, H. Spatial and temporal variability of planktonic archaeal abundance in the Humboldt Current System off Chile. Deep Sea Res. Part II 56, 1073–1082. https://doi.org/10.1016/j.dsr2.2008.09.012 (2009).Article 

Google Scholar 
Antezana, T. Euphausia mucronata: A keystone herbivore and prey of the Humboldt Current System. Deep Sea Res. Part II 57, 652–662. https://doi.org/10.1016/j.dsr2.2009.10.014 (2010).Article 

Google Scholar 
Anguita, C., Gelcich, S., Aldana, M. & Pulgar, J. Exploring the influence of upwelling on the total allowed catch and harvests of a benthic gastropod managed under a territorial user rights for fisheries regime along the Chilean coast. Ocean Coast. Manag. 195, 105256. https://doi.org/10.1016/j.ocecoaman.2020.105256 (2020).Article 

Google Scholar 
González, J. E., Yannicelli, B. & Stotz, W. The interplay of natural variability, productivity and management of the benthic ecosystem in the Humboldt Current System: Twenty years of assessment of Concholepas concholepas fishery under a TURF management system. Ocean Coast. Manag. 208, 105628. https://doi.org/10.1016/j.ocecoaman.2021.105628 (2021).Article 

Google Scholar 
Canales, T. M. et al. Endogenous, climate, and fishing influences on the population dynamics of Small Pelagic Fish in the Southern Humboldt Current Ecosystem. Front. Mar. Sci. https://doi.org/10.3389/fmars.2020.00082 (2020).Article 

Google Scholar 
González, J. E., Ortiz, M. Exploring harvest strategies in a benthic habitat in the Humboldt Current System (Chile): A study case. In Marine Coastal Ecosystems Modelling and Conservation: Latin American Experiences 127–141 (Springer International Publishing, 2021). https://doi.org/10.1007/978-3-030-58211-1_6.Ortiz, M. Pre-image population indices for anchovy and sardine species in the Humboldt Current System off Peru and Chile: Years decaying productivity. Ecol. Ind. 119, 106844. https://doi.org/10.1016/j.ecolind.2020.106844 (2020).Article 

Google Scholar 
Tognelli, M. F., Silva-Garcia, C., Labra, F. A. & Marquet, P. A. Priority areas for the conservation of coastal marine vertebrates in Chile. Biol. Conserv. 126, 420–428. https://doi.org/10.1016/j.biocon.2005.06.021 (2005).Article 

Google Scholar 
Bustamante, C., Vargas-Caro, C. & Bennett, M. B. Not all fish are equal: Functional biodiversity of cartilaginous fishes (Elasmobranchii and Holocephali) in Chile. J. Fish Biol. 85, 1617–1633. https://doi.org/10.1111/jfb.12517 (2014).CAS 
Article 
PubMed 

Google Scholar 
Sarmiento-Devia, R. A., Harrod, C. & Pacheco, A. S. Ecology and Conservation of Sea Turtles in Chile. Chelonian Conserv. Biol. 14, 21–33. https://doi.org/10.2744/ccab-14-01-21-33.1 (2015).Article 

Google Scholar 
Pérez-Álvarez, M. J., Alvarez, E., Aguayo-Lobo, A. & Olavarría, C. Occurrence and distribution of Chilean dolphin (Cephalorhynchus eutropia) in coastal waters of central Chile. N.Z. J. Mar. Freshw. Res. 41, 405–409. https://doi.org/10.1080/00288330709509931 (2007).Article 

Google Scholar 
Pacheco, A. S. et al. Cetacean diversity revealed from whale-watching observations in Northern Peru. Aquat. Mamm. 45, 116–122. https://doi.org/10.1578/AM.45.1.2019.116 (2019).Article 

Google Scholar 
Buchan, S. J., Vásquez, P., Olavarría, C. & Castro, L. R. Prey items of baleen whale species off the coast of Chile from fecal plume analysis. Mar. Mamm. Sci. 37, 1116–1127 (2021).Article 

Google Scholar 
Hucke-Gaete, R. et al. From Chilean Patagonia to Galapagos, Ecuador: Novel insights on blue whale migratory pathways along the Eastern South Pacific. PeerJ 6, e4695. https://doi.org/10.7717/peerj.4695 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Llapapasca, M. A. et al. Modeling the potential habitats of dusky, commons and bottlenose dolphins in the Humboldt Current System off Peru: The influence of non-El Niño vs. El Niño 1997–98 conditions and potential prey availability. Prog. Oceanogr. 168, 169–181. https://doi.org/10.1016/j.pocean.2018.09.003 (2018).Article 

Google Scholar 
Sepúlveda, M. et al. From whaling to whale watching: Identifying fin whale critical foraging habitats off the Chilean coast. Aquat. Conserv. Mar. Freshw. Ecosyst. 28, 821–829. https://doi.org/10.1002/aqc.2899 (2018).Article 

Google Scholar 
Williams, R. et al. Chilean blue whales as a case study to illustrate methods to estimate abundance and evaluate conservation status of rare species. Conserv. Biol. 25, 526–535. https://doi.org/10.1111/j.1523-1739.2011.01656.x (2011).Article 
PubMed 

Google Scholar 
Moore, J. E. & Barlow, J. Bayesian state-space model of fin whale abundance trends from a 1991–2008 time series of line-transect surveys in the California Current. J. Appl. Ecol. 48, 1195–1205. https://doi.org/10.1111/j.1365-2664.2011.02018.x (2011).Article 

Google Scholar 
Campbell, G. S. et al. Inter-annual and seasonal trends in cetacean distribution, density and abundance off southern California. Deep Sea Res. Part II 112, 143–157. https://doi.org/10.1016/j.dsr2.2014.10.008 (2015).Article 

Google Scholar 
Nichol, L. M., Wright, B. M., O’Hara, P. & Ford, J. K. B. Risk of lethal vessel strikes to humpback and fin whales off the west coast of Vancouver Island, Canada. Endanger. Species Res. 32, 373–390. https://doi.org/10.3354/esr00813 (2017).Article 

Google Scholar 
Pennino, M. G. et al. A spatially explicit risk assessment approach: Cetaceans and marine traffic in the Pelagos Sanctuary (Mediterranean Sea). PLoS One 12, e0179686. https://doi.org/10.1371/journal.pone.0179686 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Van Waerebeek, K. & Reyes, J. C. Catch of small cetaceans at Pucusana Port, central Peru, during 1987. Biol. Conserv. 51, 15–22. https://doi.org/10.1016/0006-3207(90)90028-N (1990).Article 

Google Scholar 
Mangel, J. C. et al. Small cetacean captures in Peruvian artisanal fisheries: High despite protective legislation. Biol. Conserv. 143, 136–143. https://doi.org/10.1016/j.biocon.2009.09.017 (2010).Article 

Google Scholar 
Campbell, E., Pasara-Polack, A., Mangel, J. C. & Alfaro-Shigueto, J. Use of small cetaceans as bait in small-scale fisheries in Peru. Front. Mar. Sci. https://doi.org/10.3389/fmars.2020.534507 (2020).Article 

Google Scholar 
Reyes, J. C. & Oporto, J. A. Gillnet fisheries and cetaceans in the southeast Pacific. Report of the International Whaling Commission 467–474 (1994).Aguayo-Lobo, A. Los cetáceos y sus perspectivas de conservación. Estudios Oceanológicos 18, 35–43 (1999).
Google Scholar 
Félix, F., Muñoz, M., Falconí, J., Botero, N., Haase, B., et al. Entanglement of humpback whales in artisanal fishing gear in Ecuador. J. Cetacean. Res. Manag. 283–290 (2020).Félix, F. et al. Challenges and opportunities for the conservation of marine mammals in the Southeast Pacific with the entry into force of the U.S. Marine Mammal Protection Act. Reg. Stud. Mar. Sci. 48, 102036. https://doi.org/10.1016/j.rsma.2021.102036 (2021).Article 

Google Scholar 
García-Cegarra, A. M. & Pacheco, A. S. Collision risk areas between fin and humpback whales with large cargo vessels in Mejillones Bay (23°S), northern Chile. Mar. Policy 103, 182–186. https://doi.org/10.1016/j.marpol.2018.12.022 (2019).Article 

Google Scholar 
Santos-Carvallo, M. et al. Impacts of whale-watching on the short-term behavior of Fin Whales (Balaenoptera physalus) in a marine protected area in the southeastern pacific. Front. Mar. Sci. https://doi.org/10.3389/fmars.2021.623954 (2021).Article 

Google Scholar 
Villagra, D., García-Cegarra, A., Gallardo, D. I. & Pacheco, A. S. Energetic effects of whale-watching boats on humpback whales on a breeding ground. Front. Mar. Sci. https://doi.org/10.3389/fmars.2020.600508 (2021).Article 

Google Scholar 
Buckland, S., Anderson, D., Burnham, K., Laake, J., Borchers, D., Thomas, L. Introduction to Distance Sampling Estimating Abundance of Biological Populations. (Oxford University Press, 2001).Hedley, S. L. & Buckland, S. T. Spatial models for line transect sampling. JABES 9, 181–199. https://doi.org/10.1198/1085711043578 (2004).Article 

Google Scholar 
Williams, R., Hedley, S. L., Hammond, P. S. Modeling distribution and abundance of Antarctic baleen whales using ships of opportunity (2006).DoniolValcroze, T., Berteaux, D., Larouche, P. & Sears, R. Influence of thermal fronts on habitat selection by four rorqual whale species in the Gulf of St. Lawrence. Mar. Ecol. Prog. Ser. 335, 207–216. https://doi.org/10.3354/meps335207 (2007).Article 

Google Scholar 
Scales, K. L. et al. Should I stay or should I go? Modelling year-round habitat suitability and drivers of residency for fin whales in the California Current. Divers. Distrib. 23, 1204–1215. https://doi.org/10.1111/ddi.12611 (2017).Article 

Google Scholar 
Bedriñana-Romano, L. et al. Integrating multiple data sources for assessing blue whale abundance and distribution in Chilean Northern Patagonia. Divers. Distrib. https://doi.org/10.1111/ddi.12739 (2018).Article 

Google Scholar 
Bedriñana-Romano, L. et al. Defining priority areas for blue whale conservation and investigating overlap with vessel traffic in Chilean Patagonia, using a fast-fitting movement model. Sci. Rep. 11, 2709. https://doi.org/10.1038/s41598-021-82220-5 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Pirotta, E., Matthiopoulos, J., MacKenzie, M., Scott-Hayward, L. & Rendell, L. Modelling sperm whale habitat preference: A novel approach combining transect and follow data. Mar. Ecol. Prog. Ser. 436, 257–272. https://doi.org/10.3354/meps09236 (2011).Article 

Google Scholar 
Mendelssohn, R. rerddapXtracto: Extracts Environmental Data from “ERDDAP” Web Services. (2020).Lau-Medrano, W. grec: Gradient-Based Recognition of Spatial Patterns in Environmental Data. (2020).Belkin, I. M. & O’Reilly, J. E. An algorithm for oceanic front detection in chlorophyll and SST satellite imagery. J. Mar. Syst. 78, 319–326. https://doi.org/10.1016/j.jmarsys.2008.11.018 (2009).Article 

Google Scholar 
Hijmans, R. J., van Etten, J., Cheng, J., Sumner, M., Mattiuzzi, M., Greenberg, J. A., et al. raster: Geographic Data Analysis and Modeling. (2018).Royle, J. A. N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115. https://doi.org/10.1111/j.0006-341X.2004.00142.x (2004).MathSciNet 
Article 
PubMed 
MATH 

Google Scholar 
Chelgren, N. D., Samora, B., Adams, M. J. & McCreary, B. Using spatiotemporal models and distance sampling to map the space use and abundance of newly metamorphosed Western Toads (Anaxyrus boreas). Herpetol. Conserv. Biol. 6, 16 (2011).
Google Scholar 
Hartig, F., Lohse, L. DHARMa: Residual Diagnostics for Hierarchical (Multi-Level/Mixed) Regression Models. (2022).Gelman, A., Meng, X.-L. & Stern, H. Posterior predictive assessment of model fitness via realized discrepancies. Stat. Sin. 6, 733–760. https://doi.org/10.2307/24306036 (1996).MathSciNet 
Article 
MATH 

Google Scholar 
Kery, M. & Royle, J. A. Applied Hierarchical Modeling in Ecology: Analysis of Distribution, Abundance and Species Richness in R and BUGS: Volume 1: Prelude and Static Models. (Academic Press, 2015).R DCT. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2015).Plummer, M. JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. (2003).Fonnesbeck, C. J., Garrison, L. P., Ward-Geiger, L. I. & Baumstark, R. D. Bayesian hierarchichal model for evaluating the risk of vessel strikes on North Atlantic right whales in the SE United States. Endanger. Species Res. 6, 87–94. https://doi.org/10.3354/esr00134 (2008).Article 

Google Scholar 
Vanderlaan, A. S. M., Taggart, C. T., Serdynska, A. R., Kenney, R. D. & Brown, M. W. Reducing the risk of lethal encounters: Vessels and right whales in the Bay of Fundy and on the Scotian Shelf. Endanger. Species Res. 4, 283–297. https://doi.org/10.3354/esr00083 (2008).Article 

Google Scholar 
Warren, D. L., Glor, R. E. & Turelli, M. Environmental niche equivalency versus conservatism: Quantitative approaches to niche evolution. Evolution 62, 2868–2883. https://doi.org/10.1111/j.1558-5646.2008.00482.x (2008).Article 
PubMed 

Google Scholar 
Hijmans, R. J., Phillips, S., Leathwick, J., Elith, J. & Hijmans, M. R. J. Package ‘dismo’. Circles 9, 1–68 (2017).
Google Scholar 
Daneri, G. et al. Primary production and community respiration in the Humboldt Current System off Chile and associated oceanic areas. Mar. Ecol. Prog. Ser. 197, 41–49. https://doi.org/10.3354/meps197041 (2000).Article 

Google Scholar 
Montecino, V. & Lange, C. B. The Humboldt Current System: Ecosystem components and processes, fisheries, and sediment studies. Prog. Oceanogr. 83, 65–79. https://doi.org/10.1016/j.pocean.2009.07.041 (2009).Article 

Google Scholar 
Escribano, R., Hidalgo, P. & Krautz, C. Zooplankton associated with the oxygen minimum zone system in the northern upwelling region of Chile during March 2000. Deep Sea Res. Part II 56, 1083–1094. https://doi.org/10.1016/j.dsr2.2008.09.009 (2009).Article 

Google Scholar 
Perez-Alvarez, M. et al. Fin whales (Balaenoptera physalus) feeding on Euphausia mucronata in nearshore waters off North-Central Chile. Aquat. Mamm. 32, 109–113. https://doi.org/10.1578/AM.32.1.2006.109 (2006).Article 

Google Scholar 
Riquelme-Bugueño, R. et al. Fatty acid composition in the endemic Humboldt Current krill, Euphausia mucronata (Crustacea, Euphausiacea) in relation to the phytoplankton community and oceanographic variability off Dichato coast in central Chile. Prog. Oceanogr. 188, 102425. https://doi.org/10.1016/j.pocean.2020.102425 (2020).Article 

Google Scholar 
Escribano, R., Marin, V. & Irribarren, C. Distribution of Euphausia mucronata at the upwelling area of Peninsula Mejillones, northern Chile: The influence of the oxygen minimum layer. Sci. Mar. 64, 69–77. https://doi.org/10.3989/scimar.2000.64n169 (2000).Article 

Google Scholar 
Riquelme-Bugueno, R., Escribano, R. & Gomez-Gutierrez, J. Somatic and molt production in Euphausia mucronata off central-southern Chile: The influence of coastal upwelling variability. Mar. Ecol. Prog. Ser. 476, 39–57 (2013).Article 

Google Scholar 
Savoca, M. S. et al. Baleen whale prey consumption based on high-resolution foraging measurements. Nature 599, 85–90. https://doi.org/10.1038/s41586-021-03991-5 (2021).CAS 
Article 
PubMed 

Google Scholar 
Roman, J. & McCarthy, J. J. The whale pump: Marine mammals enhance primary productivity in a coastal basin. PLoS One 5, e13255. https://doi.org/10.1371/journal.pone.0013255 (2010).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Hucke-Gaete, R. Whales might also be an important component in patagonian fjord ecosystems: Comment to Iriarte et al. Ambio 40, 104–105. https://doi.org/10.1007/s13280-010-0110-8 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Lavery, T. J. et al. Whales sustain fisheries: Blue whales stimulate primary production in the Southern Ocean. Mar. Mamm. Sci. https://doi.org/10.1111/mms.12108 (2014).Article 

Google Scholar 
Roman, J. et al. Whales as marine ecosystem engineers. Front. Ecol. Environ. 12, 377–385. https://doi.org/10.1890/130220 (2014).Article 

Google Scholar 
Hucke-Gaete, R., Osman, L. P., Moreno, C. A., Findlay, K. P. & Ljungblad, D. K. Discovery of a blue whale feeding and nursing ground in southern Chile. Proc. R. Soc. Lond. B 271, S170–S173. https://doi.org/10.1098/rsbl.2003.0132 (2004).Article 

Google Scholar 
Buchan, S. J. & Quiones, R. A. First insights into the oceanographic characteristics of a blue whale feeding ground in northern Patagonia, Chile. Mar. Ecol. Prog. Ser. 554, 183–199. https://doi.org/10.3354/meps11762 (2016).CAS 
Article 

Google Scholar 
Findlay, K., Pitman, R., Tsurui, T., Sakai, K., Ensor, P., Iwakami, H., et al. IWC-southern whale and ecosystem research (IWC/SOWER) blue whale Cruise, Chile. Documento Técnico, IWC 1998 (1998).Branch, T. A. et al. Past and present distribution, densities and movements of blue whales Balaenoptera musculus in the Southern Hemisphere and northern Indian Ocean. Mamm. Rev. 37, 116–175. https://doi.org/10.1111/j.1365-2907.2007.00106.x (2007).Article 

Google Scholar 
Barlow, D. R., Klinck, H., Ponirakis, D., Garvey, C. & Torres, L. G. Temporal and spatial lags between wind, coastal upwelling, and blue whale occurrence. Sci. Rep. 11, 6915. https://doi.org/10.1038/s41598-021-86403-y (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Galletti-Vernazzani, B., Jackson, J. A., Cabrera, E., Carlson, C. A. Jr. & RLB.,. Estimates of abundance and trend of chilean blue whales off Isla de Chiloé, Chile. PLoS One 12, e0168646. https://doi.org/10.1371/journal.pone.0168646 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Friedlaender, A. S., Goldbogen, J. A., Hazen, E. L., Calambokidis, J. & Southall, B. L. Feeding performance by sympatric blue and fin whales exploiting a common prey resource. Mar. Mamm. Sci. 31, 345–354. https://doi.org/10.1111/mms.12134 (2015).Article 

Google Scholar 
Abrahms, B. et al. Memory and resource tracking drive blue whale migrations. PNAS 116, 5582–5587 (2019).CAS 
Article 

Google Scholar 
Clarke, R., Aguayo, A. & Basulto, S. Whale observation and whale marking off the coast of Chile in 1964. Sci. Rep. Whales Res. Inst. Tokyo 30, 117–178 (1978).
Google Scholar 
Allison, C. IWC individual and summary catch databases Version 5.5 (12 February 2013). Available from the International Whaling Commission 135 (2013).Pastene, L. A., Acevedo, J. & Branch, T. A. Morphometric analysis of Chilean blue whales and implications for their taxonomy. Mar. Mamm. Sci. 36, 116–135. https://doi.org/10.1111/mms.12625 (2020).Article 

Google Scholar 
Rendell, L., Whitehead, H. & Escribano, R. Sperm whale habitat use and foraging success off northern Chile: Evidence of ecological links between coastal and pelagic systems. Mar. Ecol. Prog. Ser. 275, 289–295. https://doi.org/10.3354/meps275289 (2004).Article 

Google Scholar 
Jaquet, N. & Whitehead, H. Scale-dependent correlation of sperm whale distribution with environmental features and productivity in the South Pacific. Mar. Ecol. Prog. Ser. 135, 1–9. https://doi.org/10.3354/meps135001 (1996).Article 

Google Scholar 
O’Hern, J. E., Biggs, D. C. Sperm whale (Physeter macrocephalus) habitat in the Gulf of Mexico: Satellite observed ocean color and altimetry applied to small-scale variability in distribution. Aquat. Mamm. 35 (2009).Koen Alonso, M., Crespo, E. A., García, N. A., Pedraza, S. N. & Coscarella, M. A. Diet of dusky dolphins, Lagenorhynchus obscurus, in waters off Patagonia, Argentina. Fish. Bull. 96, 366–374 (1998).
Google Scholar 
García-Godos, I., Waerebeek, K. V., Reyes, J. C., Alfaro-Shigueto, J. & Arias-Schreiber, M. Prey occurrence in the stomach contents of four small cetacean species in Peru. Latin Am. J. Aquat. Mamm. 6, 171–183. https://doi.org/10.5597/lajam00122 (2007).Article 

Google Scholar 
Dans, S. L., Crespo, E. A., Koen-Alonso, M., Markowitz, T. M., Berón Vera, B., Dahood, A. D. Chapter 3—Dusky dolphin trophic ecology: Their role in the food web. In The Dusky Dolphin (eds. Würsig, B., Würsig, M.) 49–74 (Academic Press, 2010). https://doi.org/10.1016/B978-0-12-373723-6.00003-5.Romero, M. A. et al. Feeding habits of two sympatric dolphin species off North Patagonia, Argentina. Mar. Mamm. Sci. 28, 364–377 (2012).Article 

Google Scholar 
Loizaga de Castro, R. et al. Feeding ecology of dusky dolphins Lagenorhynchus obscurus: Evidence from stable isotopes. J. Mammal. 97, 310–320. https://doi.org/10.1093/jmammal/gyv180 (2016).Article 

Google Scholar 
Cipriano, F. W. Behavior and occurrence patterns, feeding ecology, and life history of dusky dolphins (Lagenorhynchus obscurus) off Kaikoura, New Zealand. (1992).Benoit-Bird, K. J., Würsig, B. & Mfadden, C. J. Dusky dolphin (lagenorhynchus obscurus) foraging in two different habitats: Active acoustic detection of dolphins and their prey. Mar. Mamm. Sci. 20, 215–231. https://doi.org/10.1111/j.1748-7692.2004.tb01152.x (2004).Article 

Google Scholar 
Van Waerebeek, K. Records of dusky dolphins Lagenorhynchus obscurus (Gray, 1828) in the eastern South Pacific. Beaufortia (1992).Selzer, L. A. & Payne, P. M. The distribution of white-sided (Lagenorhynchus acutus) and common dolphins (Delphinus delphis) vs. Environmental features of the continental shelf of the Northeastern United States. Mar. Mamm. Sci. 4, 141–153. https://doi.org/10.1111/j.1748-7692.1988.tb00194.x (1988).Article 

Google Scholar 
Neumann, D. R. Seasonal movements of short-beaked common dolphins (Delphinus delphis) in the north-western Bay of Plenty, New Zealand: Influence of sea surface temperature and El Niño/La Niña. N.Z. J. Mar. Freshw. Res. 35, 371–374. https://doi.org/10.1080/00288330.2001.9517007 (2001).Article 

Google Scholar 
Peters, K. J. et al. Foraging ecology of the common dolphin Delphinus delphis revealed by stable isotope analysis. Mar. Ecol. Prog. Ser. 652, 173–186. https://doi.org/10.3354/meps13482 (2020).CAS 
Article 

Google Scholar 
Brand, D. et al. Common dolphins, common in neritic waters off southern Israel, demonstrate uncommon dietary habits. Aquat. Conserv. Mar. Freshw. Ecosyst. 31, 15–21. https://doi.org/10.1002/aqc.3165 (2021).Article 

Google Scholar 
Barlow, J. & Taylor, B. L. Estimates of sperm whale abundance in the Northeastern temperate pacific from a combined acoustic and visual survey. Mar. Mamm. Sci. 21, 429–445. https://doi.org/10.1111/j.1748-7692.2005.tb01242.x (2005).Article 

Google Scholar 
Cañadas, A., Desportes, G. & Borchers, D. Estimation of g (0) and abundance of common dolphins (Delphinus delphis) from the NASS-95 Faroese survey. J. Cetac. Res. Manag. 6, 191–198 (2004).
Google Scholar 
Miller, D. L., Burt, M. L., Rexstad, E. A. & Thomas, L. Spatial models for distance sampling data: Recent developments and future directions. Methods Ecol. Evol. 4, 1001–1010. https://doi.org/10.1111/2041-210X.12105 (2013).Article 

Google Scholar 
Sigourney, D. B. et al. Developing and assessing a density surface model in a Bayesian hierarchical framework with a focus on uncertainty: Insights from simulations and an application to fin whales (Balaenoptera physalus). PeerJ 8, e8226. https://doi.org/10.7717/peerj.8226 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Panigada, S. et al. Mediterranean fin whales at risk from fatal ship strikes. Mar. Pollut. Bull. 52, 1287–1298. https://doi.org/10.1016/j.marpolbul.2006.03.014 (2006).CAS 
Article 
PubMed 

Google Scholar 
Ribeiro, S., Viddi, F. A. & Freitas, T. R. Behavioural responses of Chilean dolphins (Cephalorhynchus eutropia) to boats in Yaldad Bay, southern Chile. Aquat. Mamm. 31, 234 (2005).Article 

Google Scholar 
Bearzi, G. et al. Overfishing and the disappearance of short-beaked common dolphins from western Greece. Endanger. Species Res. 5, 1–12. https://doi.org/10.3354/esr00103 (2008).Article 

Google Scholar 
Reeves, R. R., McClellan, K. & Werner, T. B. Marine mammal bycatch in gillnet and other entangling net fisheries, 1990 to 2011. Endanger. Species Res. 20, 71–97. https://doi.org/10.3354/esr00481 (2013).Article 

Google Scholar 
van der Hoop, J. M. et al. Vessel strikes to large whales before and after the 2008 Ship Strike Rule. Conserv. Lett. 8, 24–32. https://doi.org/10.1111/conl.12105 (2015).Article 

Google Scholar 
Erbe, C., Reichmuth, C., Cunningham, K., Lucke, K. & Dooling, R. Communication masking in marine mammals: A review and research strategy. Mar. Pollut. Bull. 103, 15–38. https://doi.org/10.1016/j.marpolbul.2015.12.007 (2016).CAS 
Article 
PubMed 

Google Scholar 
González-But, J. C. & Sepúlveda, M. Captura incidental del delfín común (Delphinus delphis) en la pesquería industrial de cerco, norte de Chile. Rev. Biol. Mar. Oceanogr. 51, 429–433. https://doi.org/10.4067/S0718-19572016000200019 (2016).Article 

Google Scholar 
Alvarado-Rybak, M. et al. Pathological findings in cetaceans sporadically stranded along the Chilean Coast. Front. Mar. Sci. 7, 684. https://doi.org/10.3389/fmars.2020.00684 (2020).Article 

Google Scholar 
Dans, S. L., Koen, A. M., Pedraza, S. & Crespo, E. A. Incidental catch of dolphins in trawling fisheries off Patagonia, Argentina: Can populations persist?. Ecol. Appl. 13, 754–762. https://doi.org/10.1890/1051-0761(2003)013[0754:ICODIT]2.0.CO;2 (2003).Article 

Google Scholar 
Childerhouse S, Baxter A. Human interactions with dusky dolphins: A management perspective, Chapter 12. In The Dusky Dolphin (eds. Würsig, B. & Würsig, M.) 245–275 (Academic Press, 2010). https://doi.org/10.1016/B978-0-12-373723-6.00012-6.Mannocci, L. et al. Assessing the impact of bycatch on dolphin populations: The case of the common dolphin in the Eastern North Atlantic. PLoS One 7, e32615. https://doi.org/10.1371/journal.pone.0032615 (2012).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Thompson, F. N., Abraham, E. R. & Berkenbusch, K. Common dolphin (Delphinus delphis) Bycatch in New Zealand commercial trawl fisheries. PLoS One 8, e64438. https://doi.org/10.1371/journal.pone.0064438 (2013).Article 
PubMed 
PubMed Central 

Google Scholar  More

Proteobacteria are major contributors to marine microbial phosphonate cyclingDatabases for all putative sequences of genes for phosphonate production (pepM, aepY, phpC, mpnS, hepD), substrate-specific catabolism (phnAWXYZ, palA), and broad-specificity catabolism (phnIJM) were created using available public genomes from JGI IMG/MER and GORG-Tropics. Gene identity was verified by the presence of catalytically essential residues (Supplementary Table S2). Phosphonate genes were identified in 10,337 genomes of bacteria and archaea spanning over 100 unique classes, suggesting a wide variety of microorganisms mediate phosphonate production and catabolism (Fig. 2, Supplementary Dataset S1). A high proportion of all collected sequences affiliated with Proteobacteria (Gamma, Alpha, and Beta classes), averaging 52% of the production genes, 78% of substrate-specific catabolism genes, and 88% of broad-specificity catabolism genes before dereplication (Fig. 2).Fig. 2: Phosphonate gene and genome count with taxonomic distribution.Number of sequences and genomes collected for study (A, D, G) with distribution of class-level taxa for all redundant sequences (B, E, H) and marine redundant (C, F, I) sequences. Results are shown for selected genes representing phosphonate (A–C) production, (D–F) substrate-specific catabolism, and (G–I) broad-specificity catabolism. The taxa shown are the 15 classes with the highest representation across all databases.Full size imageOf the 10,337 genomes, 1556 (15%) were confirmed to be marine organisms from 35 different classes (Fig. 2, Supplementary Dataset S1). Proteobacteria had even greater representation in the subset of marine genomes, averaging 65% of marine production genes, 88% of marine substrate-specific catabolism genes, and 96% of marine broad-specificity catabolism genes from the redundant databases (Fig. 2). The dominance of Alphaproteobacteria in the marine subset may be attributed to the wide variety of Pelagibacterales bacterium captured in the database, making up 426 (27%) of the 1556 genomes involved in all three categories of phosphonate cycling. Rhodobacterales (Ascidiaceihabitans sp., Roseovarius sp., Sulfitobacter sp., Labrenzia sp., and Phaeobacter sp.) alongside Rhodospirillales (Thalassobaculum sp., Thalassospira sp., Roseospira sp., Varunaivibrio sp., and Oceanibaculum sp.) were also highly represented among the marine subset with 214 (14%) and 251 (16%) genomes, respectively (Supplemental Dataset S1), though these taxa primarily show potential for phosphonate catabolism rather than production. Vibrionales were well represented in the JGI IMG/MER marine genome subset with 107 (7%) genomes spanning 59 different species including Vibrio lentus, Vibrio breoganii, and Vibrio splendidus.Diverse taxa encode the capacity to produce phosphonate derivativesPhosphonate production is widespread and distributed throughout many different bacteria and archaea. Genes responsible for the first two steps in phosphonate production, pepM and aepY, had the broadest taxonomic distribution within the redundant databases (Shannon indices of 2.66 and 2.76) for all genes in this study, distributed with 0.59 and 0.61 evenness from 70 and 72 unique, verified classes, respectively. Their broad distribution further highlights the ubiquity and necessity of phosphonate compounds to microbial life and function across all environments. Within the marine setting, both pepM and aepY have reference sequences from 22 unique, verified classes which is the second highest class representation in the marine genome subset (Fig. 2). The marine subset of pepM and aepY also have the highest Shannon indices (1.76 and 1.92) distributed with 0.53 and 0.58 evenness, respectively. A majority (87%) of the Alphaproteobacteria phosphonate producers are Pelagibacterales bacterium with other notable taxa including Bacteria: Candidatus Actinomarinaceae, Prochlorococcus sp., Synechococcus sp., Nitrosococcus sp., and MG-I Archaea: Candidatus Nitrosomarinus catalina, Nitrosopumilus maritimus, alongside other unidentified Crenarchaeota and Thaumarchaeota genomes.The gene phpC was found in less than half the number of genomes than pepM and aepY, and encoded by fewer classes in both the general database (47) and marine subset (10). In the full databases, the distribution of retrieved phpC sequences are similar to pepM and aepY with respect to taxonomic ranking, Shannon index (2.49), and evenness (0.60) (Fig. 2A–C, Supplementary Table 5). Within the marine subset, phpC has less Shannon index (1.61) but greater evenness (0.61) than the marine subset of pepM and aepY. All three upstream phosphonate production genes (pepM, aepY, phpC) are found together within Pelagibacterales bacterium, Prochlorococcus sp., Thaumarchaeota, and Crenarchaoeta alongside other taxa such as Oceanospirillales sp., Arenimonas donghaensis, Desulfuromusa kysingii, and Cellulosilyticum lentocellum.We further investigated the relationship between pepM, aepY, and phpC by examining co-occurrence in genomes and synteny with the general, redundant databases. The first two steps in phosphonate biosynthesis are intimately linked (Fig. 3). Out of all genomes with pepM, 86% have aepY, and out of all genomes with aepY, 90% have pepM. By contrast, phpC is not as closely tied to pepM and phosphonate production. We found phpC in just over 20% of genomes with the capability of phosphonate production (Fig. 3), implying that a majority of bacterial and archaeal phosphonate production stops at the production of phosphonoacetaldehyde or 2-AEP (Fig. 1A). Furthermore, half of the phpC genes were not associated with phosphonate production, given 53% of genomes with phpC did not have pepM and 54% did not have aepY (Fig. 3). In these instances, microbes may use phpC within a 2-AEP substrate-specific catabolism operon (Fig. 3) that allows phosphonate compounds to be synthesized by transforming 2-AEP with phnW and phpC into 2-HEP (Figs. 1A and 3). By repurposing 2-AEP, individuals can still create the specific compound needed while bypassing the energetically unfavourable first step of phosphonate production.Fig. 3: Co-occurrence of phosphonate cycling genes within the same genome and examples of genetic organization of phosphonate cycling genes.The heatmap displays co-occurrence of phosphonate cycling genes. Each column represents the subset of all genomes which contain the source gene and the heatmap value represents the fraction of the source genomes which also contain the co-occurring gene. Heatmap values are not symmetrical due to differing number of genomes represented in each column, database size listed above each column. Examples for phosphonate cycling genomic neighbourhoods were chosen to maximize diversity in synteny with examples from both Bacteria and Archaea where applicable. Several phosphonate-specific ABC transport system clusters are labelled as follows: phnC = phosphonate transport system ATPase; phnD = phosphonate transport system substrate-binding; phnE = phosphonate transport system permease; phnS = 2-AEP transport system substrate-binding; phnT = 2-AEP transport system ATP-binding; phnV = 2-AEP transport system permease; palC = transport system permease; palD = transport system ATP-binding; palE = transport system permease. Genes are colour coded by: red = lyase; orange = transcriptional regulator; yellow = hydrolase; green = transferase; light blue = oxidoreductase; dark blue = transaminase; purple = kinase; pink = isomerase; brown = transport; white = synthase; black = uncharacterized protein; grey = unknown.Full size imageA narrow but diverse selection of taxa encoded MpnS, the marker gene for Mpn production and a key determinant in marine methane production. We observed distinct clades of this enzyme in autotrophic archaea and heterotrophic bacteria (Fig. 2B, C). Within the marine ecosystems, Pelagibacterales, Rhodospirillales, Rickettsiales, Oceanospirillales, Flavobacteriales, and Synechococcales are bacterial candidates for MPn production alongside Thaumarchaeota and Crenarchaeota archaeon (Fig. 2B, C). While six of the bacterial genomes with MpnS also encoded genes for phosphonate catabolism, none of the archaeal MPn producers showed capacity for catabolism (Supplementary Dataset S1). The genomic neighbourhoods for general phosphonate production (pepM, aepY, phpC) and MPn production (mpnS) in both bacteria and archaea include genes such as glycosyltransferase, lipopolysaccharide choline phosphotransferase, choline kinase, adenylyltransferase, and arylsulfatase A (Fig. 3) suggesting the potential for synthesis of (methyl)phosphonate esters [93]. This is consistent with previous analysis [29] of the Nitrosopumilus maritimus SCM1 MPn production genomic neighbourhood and biophysical evidence that MPn producing archaea synthesize an exopolysaccharide modified with MPn similar to 2-AEP modified polymers.Contrary to the diversity of the other phosphonate production databases, the hepD database has low Shannon index (0.62) and evenness (0.45) with 79% of sequences mapping to Actinomycetia including Streptomycetales and Corynebacteriales (Fig. 2B, C). The marine subset has lower Shannon index (0.28) and evenness (0.41) where all sequences derive from Pelagibacterales except one from Prochlorococcus sp. The genomic neighbourhood of HMP production may contain genes for cell surface modification such as acetyltransferase, peptidoglycan biosynthesis, and adenylylsulfate transferase, suggesting that some organisms may use HMP as a conjugate for membrane-associated or exported macromolecules similar to theories on MPn utilization. Other examples of hepD synteny contain more specific genes such as the HMP dehydrogenase or other enzymes for downstream modification (Fig. 3).Marine proteobacteria encode genes for substrate-specific and broad-specificity phosphonate catabolismGenes for marine substrate-specific phosphonate catabolism were widespread among Proteobacterial classes, and to a lesser extent amongst other classes including Bacilli, Planctomycetes, and Synechococcus (Fig. 2E,F). Marine substrate-specific catabolism has lower average Shannon index (1.00) and evenness (0.43) than the three general production genes (pepM, aepY, phpC). The most widespread of these genes was phnW, likely due to its pivotal role in 2-AEP transformations as a precursor reaction to phnAY or phnX (Fig. 1, Supplementary Table 5). Marine hydrolases for 2-AEP catabolism, phnA, phnX, and phnZ, have similar Shannon indices (mean: 1.11 ± 0.05) and evenness (mean: 0.41 ± 0.03) (Fig. 2E, F, Supplementary Table 5).While not exclusive, sequenced references demonstrate a strong taxonomic partition between Proteobacterial classes for 2-AEP catabolism pathways phnAWY and phnWX. Over 74% of marine genomes with phnAWY are Alphaproteobacteria, in particular Rhodobacterales species such as Roseovarius nubinhibens, Marivita geojedonensis, and Pelagicola litoralis. On the contrary, ~80% of marine genomes with phnWX are Gammaproteobacteria, specifically of Vibrionales, Oceanospirillales, and Alteromonadales including a wide range of species from Vibrio, Photobacterium, Marinobacterium, Halomonas, and Pseudoalteromonas.Taxonomic distribution for marine phnZ was 72% Alphaproteobacteria with Pelagibacterales making up 45% of marine phnZ sequences. Note that phnZ has the most (17) reference sequences from marine Cyanobacteriia, specifically Prochlorococcus sp., than any other phosphonate catabolizing gene. Lack of marine sequence representatives for catabolism of phosphonopyruvate by palA suggests that either the substrate is uncommon, therefore the function unnecessary, or marine microbes have other methods of catabolizing phosphonopyruvate, perhaps by the C-P lyase. Overall taxonomic distribution of phosphonate substrate-specific catabolism, specifically targeting 2-AEP, suggests said function is essential to many marine heterotrophs within Alphaproteobacteria and Gammaproteobacteria. However 2-AEP catabolism appears to be less universally important than phosphonate production to marine microbial life since the required genes are found in a less diverse selection of taxa.Genetic organization for substrate-specific catabolism genes, particularly those targeting 2-AEP, varied widely in line with the numerous options for 2-AEP catabolism (Fig. 3). Though some bacteria specialize in a single 2-AEP degradation pathway such as only containing phnWAY, others contained multiple hydrolases for 2-AEP catabolism with some incorporating phpC into a 2-AEP specific catabolism operon (Fig. 3). When a genome has two hydrolases for phosphonate catabolism, often phnZ was paired with either phnA or phnX. Co-occurrence between phnZ and either phnA or phnX ranged between 30-50%, whereas co-occurrence between phnA and phnX was between 6-12% (Fig. 3). This discrepancy in co-occurrence may be due to the metabolic similarity between phnA and phnX, where having both may be redundant. Both of these enzymes rely on phnW for 2-AEP catabolism and produce carbon metabolites, whereas phnZ does not need phnW and produces the amino acid glycine (Fig. 1B).C-P lyase genes representing substrate non-specific catabolism were overwhelmingly attributed to Alphaproteobacteria which consisted over 75% of all collected marine sequences for phnIJM (Fig. 2H, I). A wide variety of Rhodobacterales, spanning 55 different genus are the most numerous representatives, followed by Pelagibacterales and Rhodospirillales. The genes in all three databases have very high genome co-occurrence, 89–99%, as expected given all three operate within the same enzyme complex (Fig. 3). Gene co-occurrence, Shannon index, and evenness is lower for phnM than the other two C-P lyase components, phnI and phnJ, likely due to instances of organisms containing two copies of phnM where one copy lies outside the C-P lyase operon [94]. C-P lyase gene databases have lower Shannon index (mean 0.81 ± 0.11) than phosphonate production and 2-AEP substrate-specific catabolism genes (phnAWXZ) (Fig. 3D, G), suggesting broad-specificity phosphonate catabolism by the C-P lyase is a narrowly distributed function (Supplementary Table 5). Organization of C-P lyase operons held the most consistency between example genomes, likely due to the high number of genes simultaneously utilized for lyase construction. These operons encoded a consecutive string of lyase subunits, including a generic phosphonate transporter (phnCDE) and GntR transcriptional regulator (Fig. 3). C-P lyase genes had low genomic co-occurrence with all other phosphonate cycling genes with notable co-occurrence between phnW at 26%, phnZ at 21%, and phnX and 18% (Fig. 3). The low rate of co-occurrence may be due to redundancy in function for P harvesting between the C-P lyase and substrate-specific catabolism. In some cases there are instances of a substrate-specific hydrolase gene located within the C-P lyase operon (Fig. 3).Phosphonate biosynthesis genes are globally prevalent in oceans and increase in mesopelagic watersFollowing curation of phn-gene databases, we analysed 121 metagenomes and 91 metatranscriptomes from the publicly available TARA Oceans expedition (spanning samples from the Atlantic, Indian, Pacific, and Southern Ocean and Red Sea) to investigate the global potential for marine phosphonate cycling. Measuring the proportion of the community capable of performing specific tasks through metagenomics indicates the long-term selective pressures that shape P-cycling and microbial communities.Potential for phosphonate production (pepM) was globally ubiquitous across all depths, with 14–17% of the community encoding in the surface waters and deep chlorophyll maximum (DCM), increasing to 45% in mesopalagic waters (Fig. 4A, Supplementary Tables S6 and S7), highlighting the importance of phosphonate compounds to marine microbial communities. Relative abundance of phpC was 64–76% that of pepM and aepY across all depths (Fig. 4A). We observed significant increase in relative abundance between the surface and mesopelagic for pepM (ANOVA: F = 1262, p  More

DEGs under low-temperature stressThe results from the field experiments indicated that the daily mean values of A, Fvʹ/Fmʹ, ETR and Fv/Fm decreased in the LT group, the PSII function was impaired, and the photosynthetic capacity was weakened. Through the specific analysis of the “Photosynthesis” pathway (pathway ID ko00195) in the LT group, it was found that PSII, the cytochrome b6f. complex (Cyt b6f.), PSI and ATPase exhibited differential gene expressions. Figure 9 shows the structural pattern diagram for photosynthesis. The parts marked by white boxes indicate that the structure has DEGs. The gene expressions of CP43, CP47, D1 protein and Cytb559 of PSII changed. The inner peripheral antenna pigment proteins, CP43 and CP47, of PSII bind to chlorophyll. They accept the excitation energy transferred from the surrounding antenna complex and transfer this energy to the reaction centre complex. Changes in CP43 and CP47 affect the absorption and transmission of light energy. In the PSII reaction centre, light energy is converted into chemical energy. P680 absorbs light and is excited to become P680*, and then transfers electrons to pheophytin (Pheo). At the same time, the PSII oxygen-evolving complex obtains electrons from water molecules, the water molecules are split and releases oxygen and protons. As one of the two core proteins that compose the reaction centre complex, the D1 protein combines with various cofactors that are related to the original charge separation and electron transfer. The D1 protein plays an important role in the process of photosynthetic electron transfer. Studies have found that low temperatures can induce allosteric inactivation of the D1 protein, which results in changes in the structure of thylakoid membranes and hinders electron transfer8. As part of the reaction centre, Cytb559 can adjust the photoinhibition sensitivity of PSII through redox changes so that the PSII reaction centre is protected from damage9. The light energy absorption, energy conversion and electron transfer functions of PSII are impaired, which result in significant decreases in Fv/Fm to levels far below the normal value. The results of Xiangchun Song are similar to those presented in this paper: the PS II reaction centre of A. mongolicus seedlings is irreversibly inactivated or the thylakoid membrane is damaged under subzero low temperature stress, which may produce serious photoinhibition. However, Song believes that the peripheral antenna component of the optical system is more affected than the core complex at low temperatures, which was not observed in the corresponding results in this study10.Figure 9Photosynthesis of A. mongolicus under low-temperature stress. The areas outlined by white boxes indicate the differentially expressed genes in these structures.Full size imageThe gene expressions of Cyt b6, PrtD and Cyt f in Cyt b6f. changed. Cyt b6f. changes not only affect the electron transport function of photosynthesis but also affect ATP synthesis. Pheo transfers the received electrons to plastid quinone (PQ). PQ receives electrons and protons to form plastid hydroquinone (PQH2). Then, the electrons of PQH2 are transferred to plastid cyanin (PC) on PSI through Cyt b6f., and hydrogen protons are released into the cavity of the thylakoid to form a transmembrane proton gradient. The transmembrane proton gradient is the driving force for ATP synthesis.The function of PSI is to transfer electrons from PC to ferredoxin for the reduction of NADP+. Recent studies have found that PSI is more sensitive to light and more prone to selective photoinhibition than PS II under low temperature and weak light conditions11,12. The KEGG analysis results indicated that the LHCI complex, PsaF and PsaE subunits of PSI showed differential gene expressions. The main function of the LHCI light-harvesting pigment protein complex is to capture light energy. PsaF is a low-molecular-weight protein that is distributed in the membrane. Some studies have suggested that the N-terminal amino acid sequence of eukaryotic PsaF is involved in the binding of PSI and PC13. PsaE, PsaD and PsaC together form the docking site of ferredoxin on the PSI receptor side14,15. Ferredoxin and ferredoxin-NADP+ reductase in the photosynthetic electron transport chain are also affected, which results in hindrance of NADPH synthesis. The F-type H+/Na+ transport ATPase subunits also show differential gene expressions, which lead to impaired ATP synthesis. Low temperatures affect the ability to absorb light energy, transfer electrons, convert light energy into electric energy, and synthesize NADPH as well as ATP, which ultimately lead to declines in Fv’/Fm’ and ETR and impair the photosynthesis capacity of A. mongolicus.Compared with the light reaction, low temperatures have a greater impact on the dark reaction. Because the dark reaction process is composed of many complex enzymatic reactions, the enzyme activity is very susceptible to temperature. The KEGG results show that 13 related enzymes were differentially expressed in the “carbon sequestration of photosynthesis” (ko00710). The Rubisco enzyme is a key enzyme that determines the direction and efficiency of photosynthetic carbon metabolism in C3 plants and is sensitive to temperature16. The results also show that the expression levels of 10 differentially expressed genes of Rubisco enzymes all declined. In the Calvin cycle, the gene expressions of only transketolase and glyceraldehyde-3-phosphate dehydrogenase are not sensitive to temperature. In addition, the reduction phase of the dark reaction requires the use of NADPH and ATP that are produced by the light reaction. The inhibition of NADPH and ATP synthesis will inevitably affect the normal progression of the Calvin cycle.Chloroplast respiration is an O2-dependent electron transport pathway in chloroplasts. Chloroplast respiration includes the nonphotochemical reduction of PQ by NAD(P) H and the reoxidation of PQ by terminal oxidase, which can consume excess electrons to protect plants from damage due to photooxidation.Figure 10 shows the partial KEGG enrichment metabolic pathway in the LT group. There were three significant enrichment pathways related to carbohydrate metabolism: fructose and mannose metabolism (ko00051), butanoate metabolism (ko00650) and C5-branched dibasic acid metabolism (ko00660). The metabolism of fructose and mannose includes the ascorbic acid biosynthetic pathway. Ascorbic acid (ASA), also known as vitamin C, can be used as a cofactor of violaxanthin de-epoxidase to participate in the lutein cycle and consume excess light energy and protect plants from harm.Figure 10The regulatory mechanism of A. mongolicus under low-temperature stress. The white ovals represent the enriched metabolic pathways. The blue rectangles represent significantly enriched KEGG metabolic pathways. The pathways are followed by the physiological structures and substances or physiological processes in which the expressions of related genes change.Full size imageLow temperatures damage cell membranes first. Increasing the mass fraction of unsaturated fatty acids in the membrane is beneficial to improve the stability and fluidity of the membrane. Some studies have shown that the degree of unsaturation of fatty acids in adult leaves of A. mongolicus that grow naturally in the field is lower in summer and higher in autumn and winter17. The significantly enriched pathways related to unsaturated fatty acid metabolism were alpha-linolenic acid metabolism (ko00592), linoleic acid metabolism (ko00591) and arachidonic acid metabolism (ko00590). Various proteins, such as linoleate 13S-lipoxygenase and cytochrome P450 family 2 subfamily J (CYP2J), which are involved in the metabolism of linoleic acid, showed differences in their gene expressions. Linoleate 13S-lipoxygenase is a common lipoxygenase in plants that can catalyse the production of precursors of several important compounds, including jasmonic acid. CYP2J is a group of P450 haem thiolate proteins, which are mainly distributed on the endoplasmic reticulum and inner mitochondrial membrane and are involved in the synthesis of sterol hormones, including brassinosteroids. Because light systems are distributed on the thylakoid membrane, damage to this membrane will affect the progress of plant photosynthesis.Plant hormone signal transduction (ko04075) plays an important role in plant resistance to stress. Studies have shown that JAs have physiological functions, such as inducing stomatal closure, inhibiting photosynthesis, promoting respiration and promoting leaf senescence18,19. Treating plants with exogenous methyl jasmonate can induce the transcription of the heat shock protein family, increase the synthesis of antioxidants, reduce lipoxygenase activity and enhance the ability of plants to resist cold damage20.Figure 11 shows the regulatory mechanism of A. mongolicus in the HL group. The MapMan analysis results show that the DEGs of the LHCII complex and those for the assembly and maintenance of PSII are significantly changed. LHCII contains chlorophyll and carotenoids, which can capture and transmit light energy. Chlorophyll is an important photosynthetic pigment that captures light energy and drives electrons to the reaction centre. The chlorophyll molecule in the reaction centre is related to photochemical quenching. The entire chlorophyll biosynthesis process (e.g., L-glutamyl-tRNA → chlorophyll a → chlorophyll b) involves 15 enzymes. The analysis found that 4/5 of the enzymes’ expression genes were changed. Carotenoids include carotene and lutein, and their synthesis is affected by high temperatures. Lutein participates in the lutein cycle, which can dissipate excess light energy and prevent membrane lipids from being peroxidized and thus maintain the stability of the thylakoid membrane structure and protect A. mongolicus. from high temperature stress and strong light stress.Figure 11The regulatory mechanism of A. mongolicus. under high-temperature stress. The white ovals represent enriched metabolic pathways. The red rectangles represent significantly enriched KEGG metabolic pathways. The pathways are followed by the physiological structures and substances or physiological processes in which the expressions of related genes change.Full size imageThe D1 protein in the PSII reaction centre is rapidly degraded under strong light conditions. To maintain the normal physiological needs of plants, the degraded D1 protein will be replaced by the new D1 protein that is produced by the repair mechanism. The reversible inactivation of the PSII reaction centre can protect the photosynthetic system and avoid destruction. This may be the reason for the significant changes in the DEGs that are involved in the assembly and maintenance of PSII.Rubisco is the main site for high-temperature inhibition of the Calvin cycle16. The KEGG analysis found that there were 7 (4↑, 3↓) DEGs of Rubisco. SBPase catalyses the conversion of sedum heptulose-1,7-diphosphate (SBP) into sedum heptulose-7-phosphate (S7P) in the renewal phase. Under low-temperature stress, only transketolase and glyceraldehyde-3-phosphate dehydrogenase remained unchanged in the Calvin cycle. In addition, NDH-mediated cyclic electron transfer may decreased the photooxidation damage that is caused by high-temperature stress by shunting the excess electrons that were generated by the inhibition of CO2 assimilation to the chloroplast respiratory pathway21.In the HT group, the net photosynthetic rates of the leaves showed two peaks on the diurnal change curves, and there was an obvious phenomenon of midday photosynthesis depression. The daily average A values were greater than those of the CK group. These results show that A. mongolicus has a complete photosynthetic structure protection mechanism and can adapt to high-temperature environments. The pathway of significant enrichment related to carbohydrate metabolism in the HT group was the same as that in the LT group. The enrichment degrees of the fructose and mannose metabolic pathways were higher only in the HT group, and C5-branched dibasic acid metabolism and butanoate metabolism were higher in the LT group.Under high temperature and strong light conditions, the balance between production and removal of reactive oxygen species (ROS) in plant cells was broken, and large amounts of reactive oxygen species accumulated in the cells. Active oxygen can cause lipid peroxidation of the biomembrane, enlarge membrane pores, increase the permeability, and affect the spatial structures of enzymes on the membrane, which thus leads to chloroplast destruction. In severe cases, ROS will cause serious injury or even death to plants22. The gene expressions of FabH and acetyl-CoA carboxylase (ACCase) changed during the synthesis of unsaturated fatty acids in the HT group.There are two types of active oxygen scavenging mechanisms in plants. (1) The enzymatic detoxification system: superoxide dismutase (SOD), ascorbate peroxidase (APX), and catalase (CAT). (2) Nonenzymatic antioxidants: ASA, carotenoids, glutathione, mannitol, and flavonoids23.Secondary metabolites result from long-term adaptation of plants to their environments. They can improve the ability of plants to protect themselves, compete for survival, and coordinate the relationship between plants and the environment. The significant enrichment pathways related to the biosynthesis of secondary metabolites in the HT group consisted of phenylpropane biosynthesis (ko00940), flavonoid biosynthesis (ko00941) and isoflavone biosynthesis (ko00943). The phenylpropanoid biosynthesis pathway is one of the three main secondary metabolic pathways in plants. It starts from phenylalanine and generates different phenylpropane metabolites through multistep reactions, such as flavonoids, isoflavones, anthocyanins and lignin24,25. Anthocyanins can protect plants from light damage by quenching free oxygen radicals and reducing the absorption of light energy. Hughes studied 10 species of evergreen broad-leaved trees and found that red leaves containing anthocyanins always maintained higher Fv/Fm levels than green leaves. Fv’/Fm’ is related to nonphotochemical quenching. This means that trees with red leaves rely more on the light-damage defence function of anthocyanins than on the light-damage defence mediated by lutein26.Riboflavin metabolism (ko00740) and biotin metabolism (ko00780) are two significantly enriched cofactors and vitamin metabolic pathways. Riboflavin is the precursor of flavin mononucleotide (FMN) and flavin adenine dinucleotide (FAD). As a prosthetic group of flavinases, FAD participates in multiple biochemical processes, such as mitochondrial electron transport, photosynthesis, fatty acid oxidation and folate metabolism, in plants27. Riboflavin can induce antioxidant accumulations in plant cells and can also promote plant growth by affecting the ethylene signalling pathway28. Biotin (e.g., VH or VB7), as an essential cofactor for biotin-dependent carboxylase, plays an important role in the life activities of plants. Common biotin-dependent carboxylase enzymes are pyruvate carboxylase (PC) and ACCase. PC is present in the mitochondria and participates in the replenishment mechanism of the tricarboxylic acid cycle. ACCase plays a pivotal role in the feedback regulation of fatty acid synthesis and is the site of action for the feedback regulation of fatty acid synthesis29.The four pathways related to amino acid metabolism showed differences in the HT group. The enrichment degrees of each pathway were as follows: valine, leucine and isoleucine biosynthesis (ko00290)  > biosynthesis of amino acids (ko01230)  > lysine biosynthesis (ko00300)  > glycine, serine and threonine metabolism (ko00260). The branched chain amino acids, valine, leucine and isoleucine and their derivatives, are beneficial to plant growth and plant responses to stress30. As an essential amino acid, lysine metabolism affects many physiological reactions, such as the tricarboxylic acid cycle, abiotic and biotic stress responses, and starch metabolism31. The glycine, serine and threonine metabolic pathways combined with the GO enrichment results showed that the genes related to glycine catabolism and glycine dehydrogenation/decarboxylase activity changed greatly. It is known that when the activity of mitochondrial glycine decarboxylase increases, both photorespiration and photosynthesis will increase32.In terms of hormones, salicylic acid, cytokinin, and abscisic acid (ABA) can improve plant active oxygen scavenging ability. Salicylic acid can decrease the damage to seedlings due to high temperatures by improving the ability of plants to resist oxidative stress and increasing the contents of osmotic adjustment substances in cells33. Salicylic acid also has the function of delaying the degradation of D1 protein and speeding up the recovery of D1 protein when high temperatures are no longer present34. ABA can improve the heat tolerance of plants by regulating the expressions of heat stress-induced genes at the transcriptional level35.In conclusion, A. mongolicus has weak resistance to low temperatures and good adaptation to high temperatures. At the physiological level, under low-temperature stress, the proportion of Y (NO) increased, the function of PSII was damaged, and photosynthesis was inhibited. A. mongolica maintains normal physiological activities by regulating the circadian rhythm, increasing the synthesis of unsaturated fatty acids and changing the effects of plant hormones. Under high-temperature stress, A. mongolicus maintains normal photosynthesis by adjusting gsw as well as water utilization and by increasing the proportion of Y (NPQ). At the same time, A. mongolicus uses LHCII to consume excess energy, continuously assembles and maintains the normal function of PSII, and changes the types of antioxidants, such as by synthesizing anthocyanins, flavonoids, and isoflavones, to protect itself from injury. In addition, the porphyrin and chlorophyll metabolisms, carotenoid metabolism, plant hormones, amino acid metabolism, unsaturated fatty acid synthesis and other metabolic pathways that are related to the differentially expressed genes changed greatly. More

The avian plumage color gamut is much more diverse than previously estimated2. We demonstrate that hummingbird barbule structural colors contribute substantially to the total color diversity of living birds, occurring in areas of the avian color space that were sparsely occupied in Stoddard and Prum2, which most notably included saturated blues, greens, and true purples (blue + red). Such regions of the avian color space were suggested to be unoccupied because these colors are challenging to create, rather than because they might function poorly for communication2. Our results support this hypothesis because hummingbird coloration densely occupies these regions of the avian color gamut (Fig. 2d), using plumage patches that generally play particularly important roles in hummingbird communication, such as throat and crown plumage patches (Supplementary Fig. 5)16,17. The greater color diversity uncovered by our study suggests that barbule structural coloration is the most versatile class of all plumage coloration mechanisms and poses the least constraints on the evolvability of plumage color diversity. Barbule structural colors evolve through changes in the size, shape, spacing, and refractive index of barbule melanosome nanostructures, but little is known about how changes in these parameters themselves evolve18.The UV/V + green region of avian color space remains mostly unoccupied (Fig. 2c, d). It is challenging to create colors with separate reflectance peaks within the wavelength sensitivities of non-adjacent color cones because the peaks must be highly saturated to avoid stimulating neighboring cones2. However, this idea does not explain why there are far more true purple (blue + red) than UV/V + green plumage colors. Notably, birds particularly fail to fill the more UV/V regions (those closer to the UV/V vertex) of UV/V + green color space, which might indicate that it is difficult to create spectra with uv/v wavelength peaks higher than those in the m wavelengths.The differences between our methods and those of Stoddard and Prum2 likely contribute in part to the larger gamut size when comparing species data but not overall data. While the number of species included in our study was comparable to that of Stoddard and Prum2 (114 vs 111 species, respectively), we measured almost twice as many plumage patches as they did (+1600 vs. 965 patches). To prevent erroneous distortion to iridescent colors we did not average the three measurements per patch. Both studies measured six standard patches for all species and additional patches if necessary to capture other plumage color variation. The larger number of plumage patches we measured reflects how color diverse hummingbird plumages are. Our methods preserved the natural variation in hue due to iridescence and avoided the distorted flattening caused by averaging highly saturated peaks with slightly different peak hues. Although our methods are biased toward increasing variation, they are necessary to accurately capture the phenomenon of iridescent hummingbird coloration.There are multiple reasons why the hummingbird color gamut is so diverse. The size of the hummingbird color gamut, like the achieved color gamut of any clade, constitutes a combination of the history of selection on color function, the clade’s evolved capacities for color production, the age of the clade, and the number of species. Hummingbirds excel at all these criteria. The 336 species of extant hummingbirds have radiated rapidly over the last 22 million years19. Hummingbird plumage color diversity has evolved through a long history of persistent sexual and social selection on plumage coloration. Hummingbirds have polygynous breeding systems characterized by female only parental care, female mate choice, and often elaborate male courtship displays. Intersexual selection in hummingbirds has contributed to elaborate radiation in brilliant plumage coloration as well as vocalizations and non-vocal feather sounds14,16,20. Hummingbird plumage color evolution rates have even been shown to positively correlate with hummingbird speciation rates14. Furthermore, in some species, brilliant monomorphic plumage ornaments apparently function in aggressive, intra- and interspecific defense of floral resources21 and appear to be associated with socioecological features related to resource competition19. Our finding that crown and throat patches, which flash brilliantly when the head of the bird is oriented toward the observer, are more diverse in coloration than other plumage regions highlights the role of plumage coloration in direct inter-individual communication and social interactions.The mechanistic properties of hummingbird barbule structural color further explain the exceptional diversity of hummingbird plumage coloration. Hummingbird barbule structural coloration is among the most complex plumage coloration mechanisms, comprised of stacks of hollow, air-filled melanosomes, surrounded by a thin superficial, solid keratin cortex as well as sometimes superficial, miniature melanin platelets which lie just beneath this cortex9,10,11,12,13. Complex nanostructures allow for independent tuning of multiple components, and, hence, greater achievable color diversity12,18,22. Barbule structural color permits the production of any peak-reflected wavelength by varying the thickness of melanosome arrays, which can produce a diversity of single-peak spectra-hues, such as the unusual diversity of greens, blues, and blue + greens seen in hummingbirds (Fig. 2b). Hummingbird melanosomes are among the most unusual in birds in being both disc-shaped and air-filled9,10,11,12,13,23. The air in the center of hummingbird melanosomes approaches the maximum possible biological difference in refractive index (air = 1.0, melanin = ~1.7), which results in the efficient production of brilliant colors with the fewest layers of melanosomes, such that resulting spectra are narrow and near saturation13,24. Such spectra can thereby create colors that extend further in color space (Fig. 2a–c).Barbule structural color also allows for the production of plumage spectra with multiple saturated peaks, creating saturated color combinations that are not as commonly produced via other plumage coloration mechanisms. However, researchers have yet to identify exactly how hummingbird multipeak spectra are produced12,13, emphasizing the need for further analyses of the optics of hummingbird feathers. Many hummingbird melanosome arrays are non-ideal– i.e., the products of the thicknesses and refractive indices of the melanin and air cavity layers are not equal25. Non-ideal thin films can create more highly saturated, pure tone colors of the primary peak while also introducing additional, harmonic spectral peaks at shorter wavelengths25, which allows for complex reflectance spectra with multiple bright peaks within the avian visible spectrum. Also, melanosome arrays with a large average layer thickness ( >~300 nm) can create colors with fundamental interference peaks in the infrared and multiple, harmonic peaks in the avian visible range (300–700 nm). The presence of minute, superficial melanin platelets below the cortex in hummingbird barbules is also correlated with secondary, lower wavelength reflectance peaks, but the precise optical mechanism remains to be established12. These different nanostructural elements all contribute to distinctive multipeak reflectance spectra that can stimulate non-adjacent color cone combinations, which Stoddard and Prum2 identified as particularly difficult to accomplish: UV/V-purple (uv/v + s + l wavelengths; Schistes geoffroyi cheek, Fig. 4g); true purple (s + l wavelengths; Atthis ellioti gorget, Fig. 4h); UV/V-green (uv/v + m; Schistes geoffroyi crown, Fig. 4a); and UV/V-red (uv/v + l; Heliangelus viola, Fig. 4b). With multipeak spectra the potential for creating new and different colors is greatly expanded, allowing for a more versatile evolution of novel colors.Unexpectedly, the hummingbird plumage color gamut is larger in volume when modeled with the VS-type (34.2%) than with the UVS-type (29.6%) visual system. This apparently unique result contrasts notably with both Stoddard and Prum’s2 and our revised estimate of the color gamut of all birds combined– VS gamut = 40.5%; UVS gamut = 47.3%. Multiple previous analyses have shown that the UVS cone-type visual system does a more efficient job of discriminating the colors of natural objects because of the broader separation between the peak spectral sensitivities of the uv and s (blue) cone types2,26,27. Because the UVS-type visual system produces an even greater increase in color volume for a diverse plant color data set over the VS-type visual system, Stoddard and Prum2 rejected the hypothesis that the UVS-type visual system had specifically evolved to expand the diversity of avian color stimuli.However, our observations that the hummingbird plumage gamut is substantially greater in volume with the VS-visual system than with the more efficient UVS-visual system strongly suggests another hypothesis: Hummingbird plumage may have specifically evolved to be more diverse within the hummingbird VS-type color visual system via selection for highly saturated plumage colors. Given diversity in hue, the way to achieve greater color gamut volume, i.e., greater plumage color diversity, is through highly chromatic color vectors that extend toward the limits of the color space. The two visual systems map variation in wavelength to different maximum potential chroma—i.e., wavelengths with color vectors that extend toward the edges, faces, and vertices of the tetrahedron6. Color vectors that extend towards the vertices, i.e., plumage that best corresponds to a singular cone type’s peak sensitivity, have the highest maximum potential chroma because vertices are the regions furthest away from the tetrahedron’s center. Thus, hummingbird plumages may have specifically evolved to have maximum chroma within their own VS-visual system via peaks that correspond most closely to the peak sensitivities of the VS- rather than the UVS-visual system. For example, when comparing the UVS and VS plumage color gamuts for hummingbirds, it is notable that hummingbird coloration extends much further into the UV/V regions of color space for the VS-visual system (Supplementary Fig. 2). While in the VS system these color points map toward the v vertex, in the UVS-visual system they map towards the uv-s edge and the uv-s-l face. Such color vectors that contribute to expanded color volume of the VS gamut could have evolved by sexual or social selection for highly saturated plumage colors that are near in hue to the specific sensitivity peaks of hummingbird receptor cone types. Such selection could note preferences within some hummingbird species for hues with maximally possible chroma, not merely for maximal chroma of a given hue.Hummingbirds have tetrachromatic color vision with substantial sensitivity in the near ultraviolet28,29. Recently, Stoddard et al.30 used a series of elegant experiments with hummingbird feeders and LED lights to demonstrate for the first time that hummingbirds can distinguish non-spectral colors distributed throughout the tetrachromatic color space. However, the presence of this remarkably proficient four-color vision in hummingbirds poses an interesting evolutionary conundrum. Recent phylogenetic analyses have established that hummingbirds and swifts are phylogenetically embedded within the nocturnal caprimulgiforms31,32. The most parsimonious hypothesis is that the immediate ancestors of swifts and hummingbirds were extensively nocturnal for approximately 8 million years before they re-evolved diurnal ecology and behavior31. Given that an evolutionary history of nocturnality can lead to the degradation or loss of opsin genes33,34, it should be a high priority to establish what effect that ancestral nocturnality may have had on the molecular physiology and anatomy of the hummingbird color visual system.Our attempt to document the color diversity of an avian family has revealed that current estimates of the total avian color gamut are likely inaccurately low. Similar studies sampling from other color-diverse families, such as sunbirds (Nectariniidae), parrots (Psittacidae), tanagers (Thraupidae), birds of paradise (Paradiseidae), manakins (Pipridae), and starlings (Sturnidae), most of which have already been studied for their plumage coloration35,36,37,38,39, would help us obtain a better estimate of the true avian color gamut. More

Foote, M. The evolution of morphological diversity. Annu. Rev. Ecol. Syst. 28, 129–152 (1997).Article 

Google Scholar 
Félix, M. A. Phenotypic evolution with and beyond genome evolution. Curr. Top. Dev. Biol. 119, 291–347 (2016).PubMed 
Article 
CAS 

Google Scholar 
Carroll, S. B. Evo-devo and an expanding evolutionary synthesis: A genetic theory of morphological evolution. Cell 134, 25–36 (2008).CAS 
PubMed 
Article 

Google Scholar 
Harvey, P. & Pagel, M. The Comparative Method in Evolutionary Biology. (Oxford University Press, 1991).Huxley, J. S. & Teissier, G. Terminology of relative growth. Nature 137, 780–781 (1936).ADS 
Article 

Google Scholar 
Klingenberg, C. P. Size, shape, and form: Concepts of allometry in geometric morphometrics. Dev. Genes Evol. 226, 113–137 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Russell, E. S. Form and Function: A Contribution to the History of Animal Morphology. (John Murray, 1916).Goswami, A. & Polly, P. D. Methods for studying morphological integration and modularity. Paleontol. Soc. Pap. 16, 213–243 (2010).Article 

Google Scholar 
Vidal-García, M., Byrne, P. G., Roberts, J. D. & Keogh, J. S. The role of phylogeny and ecology in shaping morphology in 21 genera and 127 species of Australo-Papuan myobatrachid frogs. J. Evol. Biol. 27, 181–192 (2014).PubMed 
Article 

Google Scholar 
Erwin, D. H. Disparity: Morphological pattern and developmental context. Palaeontology 50, 57–73 (2007).Article 

Google Scholar 
Fišer, C., Robinson, C. T. & Malard, F. Cryptic species as a window into the paradigm shift of the species concept. Mol. Ecol. 27, 613–635 (2018).PubMed 
Article 

Google Scholar 
Wilson, D. E. & Mittermeier, R. A. Handbook of the Mammals of the World: Volume 8: Insectivores. vol. 8 (Lynx Edicions, 2018).Jacquet, F. et al. Phylogeography and evolutionary history of the Crocidura olivieri complex (Mammalia, Soricomorpha): From a forest origin to broad ecological expansion across Africa. BMC Evol. Biol. 15, 71. https://doi.org/10.1186/s12862-015-0344-y (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Ceríaco, L. M. P. et al. Description of a new endemic species of shrew (Mammalia, Soricomorpha) from PrÍncipe Island (Gulf of Guinea). Mammalia 79, 325–341 (2015).Article 

Google Scholar 
Nicolas, V. et al. Multilocus phylogeny of the Crocidura poensis species complex (Mammalia, Eulipotyphla): Influences of the palaeoclimate on its diversification and evolution. J. Biogeogr. 46, 871–883 (2019).Article 

Google Scholar 
Konečný, A., Hutterer, R., Meheretu, Y. & Bryja, J. Two new species of Crocidura (Mammalia: Soricidae) from Ethiopia and updates on the Ethiopian shrew fauna. J. Vertebr. Biol. 69, 20064.1. https://doi.org/10.25225/jvb.20064 (2020).Article 

Google Scholar 
Couvreur, T. L. P. et al. Tectonics, climate and the diversification of the tropical African terrestrial flora and fauna. Biol. Rev. 96, 16–51 (2021).PubMed 
Article 

Google Scholar 
Mayr, E. & O’Hara, R. J. The biogeographic evidence supporting the Pleistocene forest refuge hypothesis. Evolution 40, 55–67 (1986).PubMed 
Article 

Google Scholar 
Wiens, J. J. & Graham, C. H. Niche conservatism: Integrating evolution, ecology, and conservation biology. Annu. Rev. Ecol. Evol. Syst. 36, 519–539 (2005).Article 

Google Scholar 
Smith, T. B., Wayne, R. K., Girman, D. J. & Bruford, M. W. A role for ecotones in generating rainforest biodiversity. Science 276, 1855–1857 (1997).CAS 
Article 

Google Scholar 
Needham, A. E. & Hardy, A. C. The form-transformation of the abdomen of the female pea-crab, Pinnotheres pisum Leach. Proc. R Soc. Lond. Ser. B Biol. Sci. 137, 115–136 (1950).ADS 
CAS 

Google Scholar 
Hanken, J. & Hall, B. K. The Skull, Volume 3: Functional and Evolutionary Mechanisms. (University of Chicago Press, 1993).Hautier, L., Lebrun, R. & Cox, P. G. Patterns of covariation in the masticatory apparatus of hystricognathous rodents: Implications for evolution and diversification. J. Morphol. 273, 1319–1337 (2012).PubMed 
Article 

Google Scholar 
Aristide, L. et al. Multiple factors behind early diversification of skull morphology in the continental radiation of New World monkeys. Evolution 72, 2697–2711 (2018).PubMed 
Article 

Google Scholar 
Hardin, G. The competitive exclusion principle. Science 131, 1292–1297 (1960).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Denys, C. et al. Shrews (Mammalia, Eulipotyphla) from a biodiversity hotspot, Mount Nimba (West Africa), with a field identification key to species. Zoosystema 43, 729–757 (2021).Article 

Google Scholar 
Estevo, C. A., Nagy-Reis, M. B. & Nichols, J. D. When habitat matters: Habitat preferences can modulate co-occurrence patterns of similar sympatric species. PLoS One 12, e0179489. https://doi.org/10.1371/journal.pone.0179489 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Spaeth, P. A. Morphological convergence and coexistence in three sympatric North American species of Microtus (Rodentia: Arvicolinae). J. Biogeogr. 36, 350–361 (2009).Article 

Google Scholar 
Adams, D. C., Berns, C. M., Kozak, K. H. & Wiens, J. J. Are rates of species diversification correlated with rates of morphological evolution?. Proc. R. Soc. B Biol. Sci. 276, 2729–2738 (2009).Article 

Google Scholar 
Caumul, R. & Polly, P. D. Phylogenetic and environmental components of morphological variation: Skull, mandible, and molar shape in marmots (marmota, Rodentia). Evolution 59, 2460–2472 (2005).PubMed 
Article 

Google Scholar 
Da Silva, F. O. et al. The ecological origins of snakes as revealed by skull evolution. Nat. Commun. 9, 376. https://doi.org/10.1038/s41467-017-02788-3 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Hirano, T., Kameda, Y., Kimura, K. & Chiba, S. Substantial incongruence among the morphology, taxonomy, and molecular phylogeny of the land snails Aegista, Landouria, Trishoplita, and Pseudobuliminus (Pulmonata: Bradybaenidae) occurring in East Asia. Mol. Phylogenet. Evol. 70, 171–181 (2014).PubMed 
Article 

Google Scholar 
Ge, D., Yao, L., Xia, L., Zhang, Z. & Yang, Q. Geometric morphometric analysis of skull morphology reveals loss of phylogenetic signal at the generic level in extant lagomorphs (Mammalia: Lagomorpha). Contrib. Zool. 84, 267–284 (2015).Article 

Google Scholar 
Zou, Z. & Zhang, J. Morphological and molecular convergences in mammalian phylogenetics. Nat. Commun. 7, 12758. https://doi.org/10.1038/ncomms12758 (2016).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Ananjeva, N. B. Current state of the problems in the phylogeny of squamate reptiles (Squamata, Reptilia). Biol. Bull. Rev. 9, 119–128 (2019).Article 

Google Scholar 
Revell, L. J., Harmon, L. J. & Collar, D. C. Phylogenetic signal, evolutionary process, and rate. Syst. Biol. 57, 591–601 (2008).PubMed 
Article 

Google Scholar 
Klingenberg, C. P. & Marugán-Lobón, J. Evolutionary covariation in geometric morphometric data: Analyzing integration, modularity, and allometry in a phylogenetic context. Syst. Biol. 62, 591–610 (2013).PubMed 
Article 

Google Scholar 
Cardini, A. & Polly, P. D. Larger mammals have longer faces because of size-related constraints on skull form. Nat. Commun. 4, 2458. https://doi.org/10.1038/ncomms3458 (2013).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Esquerré, D., Sherratt, E. & Keogh, J. S. Evolution of extreme ontogenetic allometric diversity and heterochrony in pythons, a clade of giant and dwarf snakes. Evolution 71, 2829–2844 (2017).PubMed 
Article 

Google Scholar 
Marroig, G. & Cheverud, J. M. Size as a line of least evolutionary resistance: Diet and adaptive morphological radiation in New World monkeys. Evolution 59, 1128–1142 (2005).PubMed 
Article 

Google Scholar 
Cornette, R., Tresset, A., Houssin, C., Pascal, M. & Herrel, A. Does bite force provide a competitive advantage in shrews? The case of the greater white-toothed shrew. Biol. J. Linn. Soc. 114, 795–807 (2015).Article 

Google Scholar 
Rodgers, G. M., Downing, B. & Morrell, L. J. Prey body size mediates the predation risk associated with being “odd”. Behav. Ecol. 26, 242–246 (2015).Article 

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

Google Scholar 
Verschuren, D. Decadal and century-scale climate variability in tropical Africa during the past 2000 years. In Past Climate Variability Through Europe and Africa (eds. Battarbee, R. W., Gasse, F. & Stickley, C. E.) 139–158 (Springer Netherlands, 2004). https://doi.org/10.1007/978-1-4020-2121-3_8.Smith, T. B., Schneider, C. J. & Holder, K. Refugial isolation versus ecological gradients. Genetica 112, 383–398 (2001).PubMed 
Article 

Google Scholar 
Brown, W. L. Jr. & Wilson, E. O. Character displacement. Syst. Biol. 5, 49–64 (1956).
Google Scholar 
Vogel, P. et al. Genetic identity of the critically endangered Wimmer’s shrew Crocidura wimmeri. Biol. J. Linn. Soc. 111, 224–229 (2014).Article 

Google Scholar 
Esselstyn, J. A. et al. Fourteen new, endemic species of shrew (genus Crocidura) from Sulawesi reveal a spectacular island radiation. Bull. Am. Mus. Nat. Hist. 454, 1–108 (2021).Article 

Google Scholar 
Evin, A., Bonhomme, V. & Claude, J. Optimizing digitalization effort in morphometrics. Biol. Methods Protoc. 5, bpaa023. https://doi.org/10.1093/biomethods/bpaa023 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Blomberg, S. P., Garland, T. & Ives, A. R. Testing for phylogenetic signal in comparative data: Behavioral traits are more labile. Evolution 57, 717–745 (2003).PubMed 
Article 

Google Scholar 
Adams, D. C. A generalized K statistic for estimating phylogenetic signal from shape and other high-dimensional multivariate data. Syst. Biol. 63, 685–697 (2014).PubMed 
Article 

Google Scholar 
Kembel, S. W. et al. Picante: R tools for integrating phylogenies and ecology. Bioinformatics 26, 1463–1464 (2010).CAS 
PubMed 
Article 

Google Scholar 
Revell, L. J. phytools: Phylogenetic Tools for Comparative Biology (and Other Things). (2021).Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
Pebesma, E. Simple features for R: Standardized support for spatial vector data. R J. 10, 439 (2018).Article 

Google Scholar 
Oksanen, J. et al. vegan: Community Ecology Package. (2020).Dray, S., Legendre, P. & Peres-Neto, P. R. Spatial modelling: A comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecol. Model. 196, 483–493 (2006).Article 

Google Scholar 
Borcard, D., Gillet, F. & Legendre, P. Numerical Ecology with R. (Springer, 2018).Dray, S. et al. adespatial: Multivariate Multiscale Spatial Analysis. (2021).Collyer, M. & Adams, D. RRPP: Linear Model Evaluation with Randomized Residuals in a Permutation Procedure. (2021).Kassambara, A. rstatix: Pipe-Friendly Framework for Basic Statistical Tests. (2021).Borcard, D., Legendre, P. & Drapeau, P. Partialling out the spatial component of ecological variation. Ecology 73, 1045–1055 (1992).Article 

Google Scholar 
Rohlf, F. J. & Corti, M. Use of two-block partial least-squares to study covariation in shape. Syst. Biol. 49, 740–753 (2000).CAS 
PubMed 
Article 

Google Scholar 
Schlager, S., Jefferis, G. & Ian, D. Morpho: Calculations and Visualisations Related to Geometric Morphometrics. (2020). More

Study areaThe NCA encompasses an area of 8292 km2 and in 2020 had approximately 87,000 inhabitants23, who are primarily dependent on livestock for their livelihoods. It is a multiple-use area where people coexist with wildlife and livestock, and practise pastoralism with transhumance, characterised by seasonal movements of livestock for accessing resources such as grazing areas and water. The NCA comprises eleven administrative wards: Alailelai, Endulen, Eyasi, Laitole, Kakesio, Misigiyo, Ngorongoro, Naiyobi, Nainokanoka, Ngoile and Olbalbal (Fig. 1). The NCA was chosen for our study as it is known to be hyperendemic for anthrax4,17,20. In addition, informal consultations we held prior to the study, as well as tailored data collection at the community and household level, indicated that local communities have a good understanding of the disease in humans and animals, and of practices around carcass and livestock management that increase risks, particularly in certain locations and periods of the year24.Figure 1Locations of participatory mapping. Map showing the 11 administrative wards of the Ngorongoro Conservation Area in northern Tanzania and the locations where participatory mapping sessions took place (red dots). The maps were produced in QGIS 2.18.2 using data from the National Bureau of Statistics, Tanzania (http://www.nbs.go.tz/).Full size imageEthics approval and consent to participateThe study received approval from the National Institute for Medical Research, Tanzania, with reference number NIMRJHQ/R.8a/Vol. IX/2660; the Tanzania Commission for Science and Technology (numbers 2016-94-NA-2016-88 (O. R. Aminu), 2016-95-NA-2016-45 (T. L. Forde) and 2018-377-NA-2016-45 (T. Lembo)); Kilimanjaro Christian Medical University College Ethics Review Committee (certificate No. 2050); and the University of Glasgow College of Medical Veterinary & Life Sciences Ethics Committee (application number 200150152). Approval and permission to access communities and participants were also obtained from relevant local authorities. Written informed consent was obtained from all participants involved in the study. All data collected were analysed anonymously, ensuring the confidentiality of participants. All research activities were performed in accordance with relevant guidelines and regulations.Participatory mappingA participatory mapping approach based on methodology previously tested in East Africa25 was employed to define areas of anthrax risk for animals in the NCA based on community knowledge. Georeferenced maps of the NCA were produced using data from Google and DigitalGlobe (2016). The maps used datum Arc 1960/UTM zone 36S and grid intervals of 1000 km and were produced at 1:10,000 and 1:50,000 scales, in order to provide participants with a choice. Ten participatory mapping focus groups were held at ward administrative level (Fig. 1) in order to identify areas in the NCA that communities perceive as posing a high risk of anthrax. One mapping exercise was held in each ward. Ngoile and Olbalbal wards were covered at the same time and treated as one, as they had only recently (in 2015) been split from one ward (Olbalbal). Each session had between ten and thirteen participants, who consisted of village and ward administrators, animal health professionals (including community animal health workers and livestock field officers), community leaders, and selected community members. These participants represented members of the community concerned with animal health and owning livestock and, as such, were likely to hold in-depth knowledge relating to community experience of animal health and disease, including anthrax. Participants were recruited by consulting with animal health professionals as well as village and ward administrators, who gave permission to conduct the mapping sessions.The mapping sessions were conducted in Swahili and translated into English by an interpreter. Participants’ general knowledge of the area was first verified by testing whether they could correctly identify popular locations such as health centres, places of worship, markets and schools. Subsequently, participants discussed among themselves and came to a consensus about areas they considered to be at high risk of anthrax. Specifically, we asked them to identify locations they perceived as areas where they considered their animals to be at risk of being exposed to anthrax. These areas were drawn on the maps provided (Fig. 2). While they did not locate areas where the animals had succumbed to disease, we also asked for generic information on locations where anthrax outbreaks had occurred in the past to define areas that could be targeted for active surveillance of cases. In order to improve the fidelity of the data, participants defined risk areas in relation to their own locality (ward) and locations where their animals access resources. Therefore, the areas were not defined by administrative boundaries, as communities may access locations outside their wards, for instance for grazing or watering. The resulting maps were scanned, digitised and analysed as detailed in the following sections. Further detail on the participatory mapping process is provided in the Supplementary Methods (Additional File 1).Figure 2Participatory mapping of anthrax risk areas in the Ngorongoro Conservation Area. Images show (A) the set-up of a mapping session, (B) participants engaged during a session and (C) an example of a 1:50,000 scale map annotated by participants. The map was created with QGIS opensource mapping software. The basemap used was a scanned and geo-referenced full colour 1:50,000 scale topographic map produced by the Surveys & Mapping Division, Ministry of Lands, Housing & Human Settlements, Dar es Salaam, Tanzania. The grid is based on the Arc1960 UTM 36S projection and datum. The map was exported from QGIS in Acrobat Pdf format to enable it to be printed at suitable sizes for using in the fieldwork and to be manually annotated during the participatory mapping.Full size imageDigitisation of maps and generation of random pointsScanned maps were saved as PDF files and converted to high resolution TIFF files for digitisation in QGIS 2.18.2-Las Palmas free OpenSource software26. All maps were georeferenced with geographical coordinates during production and reference points were available to enable the precise mapping of all locations. The digitization was carried out using the QGIS digitizing tools and by creating polygon layers of the defined risk areas.Sourcing data on the environmental predictors of anthraxAvailable soil and environmental data (250 m grid) for Tanzania were obtained from various sources (Table 1). From the available data, we selected the following seven variables which have previously been shown to contribute to or explain the risk of anthrax based on the biology of B. anthracis (Table 1).Table 1 Environmental factors with potential to influence anthrax occurrence.Full size tableCation exchange capacity (CEC)Measured in cmol/kg, CEC is the total capacity of the soil to retain exchangeable cations such as Ca2+, Mg2+ etc. It is an inherent soil characteristic and is difficult to alter significantly. It influences the soil’s ability to hold on to essential nutrients and provides a buffer against soil acidification27. CEC has been reported to be positively correlated with anthrax risk. In addition, CEC is a proxy for calcium content, which may contribute to anthrax risk in a pH-dependent manner as explained below19,22.Predicted topsoil pH (pH)Soil pH below 6.0 (acidic soil) is thought to inhibit the viability of spores19 thus a positive effect of higher pH on the risk of anthrax is expected. It has been suggested that the exosporium of B. anthracis is negatively charged in soils with neutral to slightly alkaline pH. This negative charge attracts positively charged cations in soil, mainly calcium, enabling the spores to be firmly attached to soil particles and calcium to be maintained within the spore core, thereby promoting the viability of B. anthracis19,28.Distance to inland water bodies (DOWS)Both the distance from water and proximity to water may increase anthrax risk. Distance to inland water may indicate the degree to which an area is dry/arid. Anthrax outbreaks have been shown to occur in areas with very dry conditions19. Although anthrax occurrence has also been associated with high soil moisture, this relates more to the spore germination in the environment (a mechanism that is disputed) and the concentration of spores in moist humus that amount to an infectious dose18,29. Spores will survive much longer in soils with low moisture content19. Low moisture may also be associated with low vegetation which results in animals grazing close to the soil, increasing the risk of ingesting soil with spores. Hampson et al. reported that anthrax outbreaks occurred close to water sources in the Serengeti ecosystem of Tanzania in periods of heavy rainfall20, and Steenkamp et al. found that close proximity to water bodies was key to the transmission of B. anthracis spores in Kruger National Park, South Africa22. Water is an important resource for livestock and a large number of animals may congregate at water sources during dry seasons. The close proximity of a water source to a risk area may increase the chance of infection, particularly during periods of high precipitation which might unearth buried spores.Average enhanced vegetation index (EVI)Vegetation density may influence the likelihood of an animal ingesting soil or inhaling dust that may be contaminated with spores. Grazing animals are more likely to encounter bacteria in soil with low vegetation density20, although there is a possibility that spores can be washed onto higher vegetation by the action of water19. Vegetation index may also reflect the moisture content of soil. Arid/dry conditions favour the formation and resistance of spores in the environment, thus lower vegetation may be associated with the occurrence of anthrax.Average daytime land surface temperature (LSTD)Anthrax has been more commonly reported to occur in regions with warmer climates worldwide. Minett observed that under generally favourable conditions and at 32 °C to 37 °C, sporulation of B. anthracis occurs readily but vegetative cells are more likely to disintegrate at temperatures below 21 °C30. Another hypothesis for the association of high temperature with anthrax occurrence is altered host immune response to disease due to stress caused by elevated temperatures19. In addition, elevated temperatures are usually associated with arid areas where vegetation is low, limiting access to adequate nutrition, which in turn affects immunity. Similarly, in hotter climates where infectious diseases occur more often, host interactions with other pathogens may modulate immune response to anthrax31. In this case, a lower infectious and lethal dose of spores would be sufficient to cause infection and death, respectively19. Contact with and ingestion of soil, spores and abrasive pasture is also higher with low vegetation in hot and arid areas19,32. In boreal regions such as in northern Canada, where anthrax occurs in wood bison, and Siberia, the disease is more commonly reported in the summer19. We therefore hypothesised a positive effect of LSTD on the risk of anthrax.SlopeSpores of B. anthracis are hypothesized to persist more easily in flat landscapes that are characterised by shallow slopes19, as it is thought that wind and water may disperse spores more easily along areas with a higher slope gradient, thereby decreasing the density of spores to levels that may be insufficient to cause infection in a susceptible host. Therefore, we expected a negative relationship between slope and the risk of anthrax.Predicted topsoil organic carbon content (SOC)Organic matter (g/kg) may aid spore persistence by providing mechanical support. The negatively charged exosporium of spores is attracted to the positive charges on hummus-rich soil, thus anthrax is thought to persist in soil rich in organic matter18. Based on available evidence, we expected a positive effect of SOC on the risk of anthrax.Creating the datasetThe annotated and digitised maps yielded polygons of high-risk areas within the NCA (Fig. 3). After digitization, 5000 random points were generated33 to cover the 8292 km2 area of the NCA. This enabled us to obtain distinct points allowed by the 250 m grid resolution of the environmental variables. Points falling within the defined risk areas were selected to represent risk areas while those falling outside represented low-risk areas. Measures of the environmental characteristics associated with individual points were obtained with the ‘add Raster data to points’ feature in QGIS.Figure 3Ngorongoro Conservation Area map showing (A) defined risk areas (in red) and (B) distance to settlements. For analysis, 5000 random points were generated throughout the area; points falling within 4.26 km of human settlements (the average distance herds are moved from settlements in a day as determined through interviews of resident livestock owners) were retained for analysis (n = 2173, shown in blue in 3a). The maps were created in QGIS 2.18.2 using data from the National Bureau of Statistics, Tanzania (http://www.nbs.go.tz/).Full size imageIn order to focus on areas of greatest risk to humans and livestock and to exclude locations that are not accessible, only points within a certain range of distance from settlements were included (Fig. 3). On average, herders in the NCA move their livestock 4.26 km away from settlements for grazing and watering during the day (unpublished data obtained through a cross-sectional survey of 209 households). Thus, only points falling within this distance from settlements were selected, providing us with data on areas where infection is most likely to occur. Data on locations of settlements were obtained from satellite imagery and included permanent residences as well as temporary settlements (e.g. seasonal camps set up after long distance movement away from permanent settlements, typically in the dry season, in search of pasture and water). These data were collated from the Center for International Earth Science Information Network (CIESIN).After adjusting for accessibility of resource locations using the average distance moved by livestock, 2173 points were retained for analysis, of which 239 (11%) fell within high-risk areas.Data analysisAll statistical analyses were carried out in R (v 4.1.0) within the RStudio environment34. The aims of the statistical analysis were to infer the relationship between anthrax risk areas as determined through participatory mapping and the environmental factors identified in Table 1, and to use this relationship to make spatial predictions of anthrax risk across the study area. We achieved both aims by modelling the binary risk status (high or low) of the randomly generated points as a function of their environmental characteristics in a Bayesian spatial logit-binomial generalised linear mixed-effects model (GLMM), implemented in the package glmmfields35. Spatial autocorrelation (residual non-independence between nearby points) was accounted for by including spatial random effects in the GLMM. We chose relatively non-informative priors for the intercept and the covariates, using Student’s t-distributions centred at 0 and wide variances (intercept: df = 3, location = 0, scale = 10; betas: df = 3, location = 0, scale = 3). For the spatial Gaussian Process and the observation process scale parameters, we adopted the default glmmfields settings and used half-t priors (both gp_theta and gp_sigma: df = 3, location = 0, scale = 5), and 12 knots. To achieve convergence, the models were run for 5000 iterations35.First, univariable models were fitted to estimate unadjusted associations between each environmental factor (CEC, pH, DOWS, EVI, LSTD, slope, and SOC; Table 1; Supplementary Table S1) and high- and low-risk areas. Second, we constructed multivariable models by fitting multiple environmental variables (Supplementary Table S2). Three variables, SOC, slope and EVI showed a strongly right-skewed distribution and were therefore log-transformed prior to GLMM analysis to prevent excessive influence of outliers. All predictor variables were centred to zero mean and scaled to unit standard deviation for analysis, and odds ratios were rescaled back to the original units for ease of interpretation. Prior to fitting the multivariable GLMM, the presence of collinearity among the predictor variables—which were all continuous—was assessed using variance inflation factors (VIFs)36, calculated with the car package and illustrated using scatter plots (Supplementary Fig. S1)36. Three predictor variables showed a VIF greater than 3 (LSTD, ln EVI and pH with VIFs of 6.8, 4.2 and 3.5, respectively). Removal of LSTD and ln EVI reduced all VIFs to below 3, therefore these two variables were excluded from the multivariable regression analysis37.The model performance was assessed by calculating the area under the receiver operating characteristic curve. The predicted probability of being an anthrax high-risk area was determined and depicted on a map of the NCA using a regular grid of points generated throughout the NCA with one point sampled every 500 m.Consent for publicationPermission to publish was granted by the National Institute for Medical Research, Tanzania. More