Philippa Kaur

Seymour, C. L. et al. Caught on camera: The impacts of urban domestic cats on wild prey in an African city and neighbouring protected areas. Glob. Ecol. Conserv. 23, e01198 (2020).Article 

Google Scholar 
Mori, E. et al. License to Kill? Domestic Cats Affect a Wide Range of Native Fauna in a Highly Biodiverse Mediterranean Country. Front. Ecol. Evol. 7, 477 (2019).Kays, R. et al. The small home ranges and large local ecological impacts of pet cats. Anim. Conserv. 23, 516–523 (2020).Loss, S. R., Will, T. & Marra, P. P. The impact of free-ranging domestic cats on wildlife of the United States. Nat. Commun. 4, 1396 (2013).ADS 
Article 

Google Scholar 
Van Heezik, Y., Smyth, A., Adams, A. & Gordon, J. Do domestic cats impose an unsustain386 able harvest on urban bird populations?. Biol. Conserv. 143, 121–130 (2010).Article 

Google Scholar 
Woods, M., McDonald, R. A. & Harris, S. Predation of wildlife by domestic cats Felis catus in Great Britain. Mammal Rev. 33, 174–188 (2003).Article 

Google Scholar 
Li, Y. et al. Estimates of wildlife killed by free-ranging cats in China. Biol. Conserv. 253, 108929 (2021).Article 

Google Scholar 
Barratt, D. G. Home range size, habitat utilisation and movement patterns of suburban and farm cats Felis catus. Ecography 20, 271–280 (1997).Article 

Google Scholar 
Moseby, K. E., Stott, J. & Crisp, H. Movement patterns of feral predators in an arid environment–implications for control through poison baiting. English. Wildl. Res. 36, 422–435 (2009).Article 

Google Scholar 
Hall, C. M. et al. Factors determining the home ranges of pet cats: A meta-analysis. Biol. Conserv. 203, 313–320 (2016).Article 

Google Scholar 
Castañeda, I. et al. Trophic patterns and home-range size of two generalist urban carnivores: A review. J. Zool. 307, 79–92 (2019).Article 

Google Scholar 
Hebblewhite, M. & Haydon, D. T. Distinguishing technology from biology: A critical review of the use of GPS telemetry data in ecology. Philos. Trans. R. Soc. B Biol. Sci. 365, 2303–2312 (2010).Article 

Google Scholar 
Allen, A. M. et al. Scaling up movements: From individual space use to population patterns. Ecosphere 7, e01524 (2016).
Google Scholar 
Trouwborst, A., McCormack, P. C. & Martínez Camacho, E. Domestic cats and their impacts on biodiversity: A blind spot in the application of nature conservation law. People Nat. 2, 235–250 (2020).Article 

Google Scholar 
Sims, V., Evans, K. L., Newson, S. E., Tratalos, J. A. & Gaston, K. J. Avian assemblage structure and domestic cat densities in urban environments. Divers. Distrib. 14, 387–399 (2008).Article 

Google Scholar 
Lepczyk, C. A., Mertig, A. G. & Liu, J. Landowners and cat predation across rural-to-urban landscapes. Biol. Conserv. 115, 191–201 (2004).Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing R Foundation for Statistical Computing (Vienna, Austria, 2021).Heggøy, O. & Shimmings, P. Huskattens predasjon på fugler i Norge. En vurdering basert på en litteraturgjennomgang tech. rep. 36 (2018).Morgan, S. et al. Urban cat (Felis catus) movement and predation activity associated with a wetland reserve in New Zealand. Wildl. Res. 36, 574–580 (2009).Calver, M., Grayson, J., Lilith, M. & Dickman, C. Applying the precautionary principle to the issue of impacts by pet cats on urban wildlife. Biol. Conserv. 144, 1895–1901 (2011).Article 

Google Scholar 
Crowley, S., Cecchetti, M. & Mcdonald, R. Diverse perspectives of cat owners indicate bar riers to and opportunities for managing cat predation of wildlife. Front. Ecol. Environ. 18, 544–549 (2020).Treves, A., Krofel, M., Ohrens, O. & van Eeden, L. M. Predator control needs a standard of unbiased randomized experiments with cross-over design. Front. Ecol. Evol. 7, 462 (2019).Ferreira, G. A., Machado, J. C., Nakano-Oliveira, E., Andriolo, A. & Genaro, G. The effect of castration on home range size and activity patterns of domestic cats living in a natural area in a protected area on a Brazilian island. Appl. Anim. Behav. Sci. 230, 105049 (2020).Bengsen, A. J. et al. Feral cat home-range size varies predictably with landscape productivity and population density. J. Zool. 298, 112–120 (2016).Article 

Google Scholar 
López-Jara, M. J. et al. Free-roaming domestic cats near conservation areas in Chile: Spatial movements, human care and risks for wildlife. Perspect. Ecol. Conserv. 19, 387–398 (2021).Gillies, C. & Clout, M. The prey of domestic cats (Felis catus) in two suburbs of Auckland City, New Zealand. J. Zool. 259, 309–315 (2003).Article 

Google Scholar 
Pirie, T. J., Thomas, R. L. & Fellowes, M. D. E. Pet cats (Felis catus) from urban boundaries use different habitats, have larger home ranges and kill more prey than cats from the suburbs. Landsc. Urban Plan. 220, 104338 (2022).Article 

Google Scholar 
Vucetich, J. A., Hebblewhite, M., Smith, D. W. & Peterson, R. O. Predicting prey population dynamics from kill rate, predation rate and predator-prey ratios in three wolf-ungulate systems. J. Anim. Ecol. 80, 1236–1245 (2011).Article 

Google Scholar 
Kennedy, M., Phillips, B. E. N. L., Legge, S., Murphy, S. A. & Faulkner, R. A. Do dingoes suppress the activity of feral cats in northern Australia?. Austral Ecol. 37, 134–139 (2012).Article 

Google Scholar 
Crooks, K. R. & Soule, M. E. Mesopredator release and avifaunal extinctions in a fragmented system. English. Nature 400, 563–566 (1999).ADS 
CAS 
Article 

Google Scholar 
Ferreira, J. P., Leita, O. I., Santos-Reis, M. & Revilla, E. Human-related factors regulate the spatial ecology of domestic cats in sensitive areas for conservation. PLOS ONE 6, e25970 (2011).ADS 
CAS 
Article 

Google Scholar 
Brook, L. A., Johnson, C. N. & Ritchie, E. G. Effects of predator control on behaviour of an apex predator and indirect consequences for mesopredator suppression. J. Appl. Ecol. 49, 1278–1286 (2012).Article 

Google Scholar 
Laundre, J. W., Hernandez, L. & Altendorf, K. B. Wolves, elk, and bison: Reestablishing the “landscape of fear’’ in Yellowstone National Park, USA. English. Can. J. Zool. 79, 1401–1409 (2001).Article 

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

Google Scholar 
Milleret, C. et al. GPS collars have an apparent positive effect on the survival of a large carnivore. Biol. Lett. 17, 20210128 (2021).Cecchetti, M., Crowley, S. L., Goodwin, C. E. D. & McDonald, R. A. Provision of high meat content food and object play reduce predation of wild animals by domestic cats Felis catus. Curr. Biol. 31, 1107-1111.e5 (2021).CAS 
Article 

Google Scholar 
Linklater, W., Farnworth, M., van Heezik, Y., Stafford, K. & Macdonald, E. Prioritizing cat owner behaviors for a campaign to reduce wildlife depredation. Conserv. Sci. Pract. 1, 1:e29 (2019).Selinske, M. J. et al. Identifying and prioritizing human behaviors that benefit biodiversity. Conserv. Sci. Pract. 2, e249 (2020).
Google Scholar 
McDonald, J. L., Maclean, M., Evans, M. R. & Hodgson, D. J. Reconciling actual and perceived rates of predation by domestic cats. Ecol. Evol. 5, 2745–2753 (2015).Article 

Google Scholar 
Bischof, R. et al. Estimating and forecasting spatial population dynamics of apex predators using transnational genetic monitoring. Proc. Natl. Acad. Sci. 117, 30531–30538 (2020).CAS 
Article 

Google Scholar 
Bischof, R., Gjevestad, J. G. O., Ordiz, A., Eldegard, K. & Milleret, C. High frequency GPS bursts and path-level analysis reveal linear feature tracking by red foxes. Sci. Rep. 9, 8849 (2019).ADS 
Article 

Google Scholar 
Gupte, P. R. et al. A guide to pre-processing high-throughput animal tracking data. J. Anim. Ecol. 91, 287–307 (2022).Article 

Google Scholar 
Morris, G. & Conner, L. Assessment of accuracy, fix success rate, and use of estimated horizontal position error (EHPE) to filter inaccurate data collected by a common commercially available GPS logger. PLoS ONE 12, e0189020 (2017).Article 

Google Scholar 
Clapp, J. G., Holbrook, J. D. & Thompson, D. J. GPSeqClus: An R package for sequential clustering of animal location data for model building, model application and field site investigations. Methods Ecol. Evol. 12, 787–793 (2021).Article 

Google Scholar 
Nielson, M., R., Sawyer, H. & McDonald, T. L. BBMM: Brownian Bridge Movement Model R Package Version 3.0 (2013).Horne, J. S., Garton, E. O., Krone, S. M. & Lewis, J. S. Analyzing animal movements using Brownian bridges. Ecology 88, 2354–2363 (2007).Article 

Google Scholar 
Sawyer, H., Kauffman, M. J., Nielson, R. M. & Horne, J. S. Identifying and prioritizing ungulate migration routes for landscape-level conservation. Ecol. Appl. 19, 2016–2025 (2009).Article 

Google Scholar 
Fischer, J. W., Walter, W. D. & Avery, M. L. Brownian bridge movement models to characterize birds’ home ranges. Condor 115, 298–305 (2013).Article 

Google Scholar 
Seidler, R., Long, R., Berger, J., Bergen, S. & Beckmann, J. Identifying impediments to long-distance mammal migrations. Conserv. Biol. 29 (2014).Collins, G. Seasonal distribution and routes of pronghorn in the Northern Great Basin. West. N. Am. Nat. 76, 101–112 (2016).Article 

Google Scholar  More

Abbott RJ, Brennan AC (2014) Altitudinal gradients, plant hybrid zones and evolutionary novelty. Philos Trans R Soc B Biol Sci 369:6–9Article 

Google Scholar 
Avise JC (2000) Phylogeography: the history and formation of species. Harvard University Press, Cambridge, MABook 

Google Scholar 
Baldassarre DT, White TA, Karubian J, Webster MS (2014) Genomic and morphological analysis of a semipermeable avian hybrid zone suggests asymmetrical introgression of a sexual signal. Evolution 68:2644–2657PubMed 
Article 
PubMed Central 

Google Scholar 
Barr KR, Dharmarajan G, Rhodes OE, Lance R, Leberg PL (2007) Novel microsatellite loci for the study of the black-capped vireo (Vireo atricapillus). Mol Ecol Notes 7:1067–1069CAS 
Article 

Google Scholar 
Barton NH, Gale KS (1993) Hybrid zones and the evolutionary process. In: Harrison RG (ed.) Hybrid Zones and the Evolutionary Process. Oxford University Press, New York, NY
Google Scholar 
Barton NH, Hewitt GM (1989) Adaption, speciation and hybrid zones. Nature 341:497–503CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Billerman SM, Murphy MA, Carling MD (2016) Changing climate mediates sapsucker (Aves: Sphyrapicus) hybrid zone movement. Ecol Evol 6:7976–7990PubMed 
PubMed Central 
Article 

Google Scholar 
Bell RC, Irian CG (2019) Phenotypic and genetic divergence in reed frogs across a mosaic hybrid zone on São Tomé Island. Biol J Linn Soc 128:672–680Article 

Google Scholar 
Bensch S, Price T, Kohn J (1997) Isolation and characterization of microsatellite loci in a Phylloscopus warbler. Mol Ecol 6:91–92CAS 
PubMed 
Article 

Google Scholar 
Bradbury IR, Bowman S, Borza T, Snelgrove PVR, Hutchings JA, Berg PR et al. (2014) Long distance linkage disequilibrium and limited hybridization suggest cryptic speciation in Atlantic cod. PLoS ONE 9:e106330Article 
CAS 

Google Scholar 
Brelsford A, Irwin DE (2009) Incipient speciation despite little assortative mating: the yellow-rumped warbler hybrid zone. Evolution 63:3050–3060PubMed 
Article 

Google Scholar 
Burg TM, Croxall JP (2004) Global population structure and taxonomy of the wandering albatross species complex. Mol Ecol 13:2345–2355CAS 
PubMed 
Article 

Google Scholar 
Carling MD, Zuckerberg B (2011) Spatio-temporal changes in the genetic structure of the Passerina bunting hybrid zone. Mol Ecol 20:1166–1175PubMed 
Article 

Google Scholar 
Carling MD, Thomassen HA (2012) The role of environmental heterogeneity in maintaining reproductive isolation between hybridizing Passerina (Aves: Cardinalidae) buntings. Int J Ecol 2012:295463Article 

Google Scholar 
Carpenter AM, Graham BA, Spellman GM, Klicka J, Burg TM (2021) Genetic, bioacoustic and morphological analyses reveal cryptic speciation in the warbling vireo complex (Vireo gilvus: Vireonidae: Passeriformes). Zool J Linn Soc zlab036 https://doi.org/10.1093/zoolinnean/zlab036Cicero C, Johnson NK (1998) Molecular phylogeny and ecological diversification in a clade of New World songbirds (genus Vireo). Mol Ecol 7:1359–1370CAS 
PubMed 
Article 

Google Scholar 
Chenuil A, Cahill AE, Délémontey N, Du Salliant du Luc E, Fanton H (2019) Problems and questions posed by cryptic species. A framework to guide future studies. Assessing to conserving biodiversity. History, philosophy and theory of the life sciences, Vol. 24. Springer. Daubenmire, Cham
Google Scholar 
Cheviron ZA, Brumfield RT (2012) Genomic insights into adaptation to high-altitude environments. Heredity 108:354–361CAS 
PubMed 
Article 

Google Scholar 
Coyne JA, Orr HA (2004) Speciation. Sinauer and Associates, Sunderland, Massachusetts
Google Scholar 
Culumber ZW, Shepard DB, Colemans SW, Rosenthal GG, Tobler M (2012) Physiological adaptation along environmental gradients and replicated hybrid zone structure in swordtails (Teleostei: Xiphophorus). J Evol Biol 25:1800–1814CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Dubay SG, Witt CC (2014) Differential high-altitude adaptation and restricted gene flow across a mid-elevation hybrid zone in Andean tit-tyrant flycatchers. Mol Ecol 23:3551–3565PubMed 
Article 
PubMed Central 

Google Scholar 
Garroway CJ, Bowman J, Cascaden TJ, Holloway GL, Mahan CG, Malcolm JR et al. (2010) Climate change induced hybridization in flying squirrels. Glob Chang Biol 16:113–121Article 

Google Scholar 
Grabenstein KC, Taylor SA (2018) Breaking barriers: Causes, consequences, and experimental utility of human-mediated hybridization. Trends Ecol Evol 33:198–212PubMed 
Article 

Google Scholar 
Graham BA, Cicero C, Strickland D, Woods JG, Coneybeare H, Dohms KM et al. (2021) Cryptic genetic diversity and cytonuclear discordance characterize contact among Canada jay (Perisoreus canadensis) morphotypes in western North America. Biol J Linn Soc 132:725–740Article 

Google Scholar 
Hammer Ø, Harper DA, Ryan PD (2001) Paleontological statistics software package for education and data analysis. Palaeontol Electron 4:9Haselhorst MSH, Parchman TL, Buerkle CA (2019) Genetic evidence for species cohesion, substructure and hybrids in spruce. Mol Ecol 28:2029–2045PubMed 
Article 

Google Scholar 
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A (2005) Very high resolution interpolated climate surfaces for global land areas. Int J Climatol 25:1965–1978Article 

Google Scholar 
Hawley DM (2005) Isolation and characterization of eight microsatellite loci from the house finch (Carpodactus mexicanus). Mol Ecol Notes 5:443–445CAS 
Article 

Google Scholar 
Hebert PDN, Stoeckle MY, Zemlak TS, Francis CM (2004) Identification of birds through DNA barcodes. PLoS Biol 2:e312PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Hewitt GM (1988) Hybrid zones-natural laboratories for evolutionary studies. Trends Ecol Evol 3:158–167CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Hewitt GM (2001) Speciation, hybrid zones and phylogeography—or seeing genes in space and time. Mol Ecol 10:537–549CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A (2005) Very high resolution interpolated climate surfaces for global land areas. Int J Climatol 25:1965–1978Article 

Google Scholar 
Hindley JA, Graham BA, Pulgarin-R PC, Burg TM (2018) The influence of latitude, geographic distance, and habitat discontinuities on genetic variation in a high latitude montane species. Sci Rep. 8:11846CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hinojosa JC, Koubínová D, Szenteczki MA, Pitteloud C, Dincă V, Alvarez N et al. (2019) A mirage of cryptic species: Genomics uncover striking mitonuclear discordance in the butterfly Thymelicus sylvestris. Mol Ecol 28:3857–3868PubMed 
Article 
PubMed Central 

Google Scholar 
Hubisz MJ, Falush D, Stephens M, Pritchard JK (2009) Inferring weak population structure with the assistance of sample group information. Mol Ecol Res 9:1322–1332Article 

Google Scholar 
Irwin DE (2020) Assortative mating in hybrid zones is remarkably ineffective in promoting speciation. Evolution 195:E150–E167
Google Scholar 
Johnson NK (1995) Speciation in vireos. I. Macrogeographic patterns of allozymic variation in the Vireo solitarius complex in the contiguous United States. Condor 97:903–919Article 

Google Scholar 
Johnson NK, Cicero C (2004) New mitochondrial DNA data affirm the importance of Pleistocene speciation in North American birds. Evolution 58:1122–1130PubMed 
Article 
PubMed Central 

Google Scholar 
Larson EL, Tinghitella RM, Taylor SA (2019) Insect hybridization and climate change. Front Ecol Evol 7:348Article 

Google Scholar 
Legendre P, Fortin M-J (2010) Comparison of the Mantel test and alternative approaches for detecting complex multivariate relationships in the spatial analysis of genetic data. Mol Ecol Resour 10:831–844PubMed 
Article 

Google Scholar 
Lovell SF, Lein MR, Rogers SM (2021) Cryptic speciation in the warbling vireo (Vireo gilvus). Ornithology 138:ukaa071Article 

Google Scholar 
MacDonald ZG, Dupuis JR, Davis CS, Acorn JH, Nielsen SE, Sperling FAH (2020) Gene flow and climate-associated genetic variation in a vagile habitat specialist. Mol Ecol 29:3889–3906PubMed 
Article 

Google Scholar 
Manthey JD, Klicka J, Spellman GM (2011) Cryptic diversity in a widespread North American songbird: phylogeography of the brown creeper (Certhia americana). Mol Phylogenet Evol 58:502–512PubMed 
Article 

Google Scholar 
Marchetti K, Price T, Richman A (1995) Correlates of wing morphology with foraging behaviour and migration distance in the genus Phylloscopus. J Av Biol 26:177–181Article 

Google Scholar 
Martin H, Touzet P, Dufay M, Gode C, Schmitt E, Lahiani E et al. (2017) Lineages of Silene nutans developed rapid, strong, asymmetric postzygotic reproductive isolation in allopatry. Evolution 71:1519–1531CAS 
PubMed 
Article 

Google Scholar 
Martinez JG, Soler JJ, Soler M, Moller AP, Burke T (1999) Comparative population structure and gene flow of a brood parasite, the great spotted cuckoo (Clamator glandarius) and its primary host, the magpie (Pica pica). Evolution 53:269–278CAS 
PubMed 
PubMed Central 

Google Scholar 
Mettler RD, Spellman GM (2009) A hybrid zone revisited: Molecular and morphological analysis of the maintenance, movement, and evolution of a Great Plains avian (Cardinalidae: Pheucticus) hybrid zone. Mol Ecol 18:3256–3267CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Meirmans PG, Van Tienderen PH (2004) GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms. Mol Ecol Notes 4:792–794Article 

Google Scholar 
Nowakowski JK, Szulc J, Remisiewicz M (2014) The further the flight, the longer the wing: relationship between wing length and migratory distance in Old World reed and bush warblers (Acrocephalidae and Locustellidae). Ornis Fennica 91:178–186
Google Scholar 
Pavolova A, Amos JN, Joseph L, Loynes K, Austin JJ, Keogh JS et al. (2013) Perched at the mito-nuclear crossroads: divergent mitochondrial lineages correlate with environment in the face of ongoing nuclear gene flow in an Australian bird. Evol 67:3412–3428Article 
CAS 

Google Scholar 
Piertney SB, Marquiss M, Summers R (1998) Characterization of tetranucleotide microsatellite markers in the Scottish crossbill (Loxia scotica). Mol Ecol 7:1261–1263CAS 
PubMed 
Article 

Google Scholar 
Porras-Hurtado L, Ruiz Y, Santos C, Phillips C, Carracedo A, Lareu MV (2013) An overview of STRUCTURE: Applications, parameter settings, and supporting software. Front Genet 4:98PubMed 
PubMed Central 
Article 
CAS 

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

Google Scholar 
Reding DM, Castañeda-Rico S, Shirazi S, Hofman CA, Cancellare IA, Lance SL et al. (2021) Mitochondrial genomes of the United States distribution of gray fox (Urocyon cinereoargenteus) reveal a major phylogeographic break at the Great Plains suture zone. Front Ecol Evol. https://doi.org/10.3389/fevo.2021.666800.Richardson DS, Jury FL, Dawson DA, Salgueiro P, Komdeur J, Burke T (2003) Fifty Seychelles warbler (Acrocephalus sechellensis) microsatellite loci polymorphic in Sylviidae species and their cross-species amplification in other passerine birds. Mol Ecol 9:2225–2230Article 

Google Scholar 
Riordan EC, Gugger PF, Ortego J, Smith C, Gaddis K, Thompson P et al. (2016) Association of genetic and phenotypic variability with geography and climate in three southern California oaks. Am J Bot 103:73–85PubMed 
Article 

Google Scholar 
Rush AC, Cannings RJ, Irwin DE (2009) Analysis of multilocus DNA reveals hybridization in a contact zone between Empidonax flycatchers. J Avian Biol 40:614–624Article 

Google Scholar 
Sartor CC, Cushman SA, Wan HY, Kretschmer R, Pereira JA, Bou N et al. (2021) The role of the environment in the spatial dynamics of an extensive hybrid zone between two neotropical cats. J Evol Biol 34:614–627PubMed 
Article 

Google Scholar 
Schukman JM, Lira-Noriega A, Townsend Peterson A (2011) Multiscalar ecological characterization of Say’s and eastern phoebes and their zone of contact in the Great Plains. Condor 113:372–384Article 

Google Scholar 
Seehausen O, Takimoto G, Roy D, Jokela J (2008) Speciation reversal and biodiversity dynamics with hybridization in changing environments. Mol Ecol 17:30–44PubMed 
Article 
PubMed Central 

Google Scholar 
Semenchuk GP (1992) The Atlas of Breeding Birds of Alberta. Fed. of Alberta Naturalists, Edmonton, p 243
Google Scholar 
Peakall R, Smouse PE (2012) GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research–an update. Bioinformatics 28:2537–2539CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Sorenson MD, Ast JC, Dimcheff DE, Yuri T, Mindell DP (1999) Primers for a PCR-based approach to mitochondrial genome sequencing in birds and other vertebrates. Mol Phylogent Evol 12:105–114CAS 
Article 

Google Scholar 
Spellman GM, Klicka J (2007) Phylogeography of the white-breasted nuthatch (Sitta carolinensis): diversification in North American pine and oak woodlands. Mol Ecol 16:1729–1740CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Stenzler LM, Fitzpatrick JW (2002) Isolation of microsatellite loci in the Florida scrub jay Aphelocoma coerulescens. Mol Ecol Notes 2:547–550CAS 
Article 

Google Scholar 
Swenson NG (2006) GIS-based niche models reveal unifying climatic mechanisms that maintain location of avian hybrid zones in a North America suture zone. J Evol Biol. 19:717–725CAS 
PubMed 
Article 

Google Scholar 
Swenson NG, Howard DJ (2005) Clustering of contact zones, hybrid zones, and phylogeographic breaks in North America. Am Nat 166:581–591PubMed 
Article 

Google Scholar 
Tarr CL, Fleischer RC (1998) Primers for polymorphic GT microsatellites isolated from the Mariana crow, Corvus kubaryi. Mol Ecol 7:253–255CAS 
PubMed 
Article 

Google Scholar 
Tarroso P, Pereira RJ, Martínez-Freiría F, Godinho R, Brito JC (2014) Hybridization at an ecotone: Ecological and genetic barriers between three Iberian vipers. Mol Ecol 23:1108–1123CAS 
PubMed 
Article 

Google Scholar 
Taylor SA, Larson EL, Harrison RG (2015) Hybrid zones: windows on climate change. Trends Ecol Evol 30:398–406PubMed 
PubMed Central 
Article 

Google Scholar 
Toews DPL, Mandic M, Richards JG, Irwin DE (2014) Migration, mitochondria and the yellow-rumped warbler. Evolution 68:241–255CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Toews DPL, Campagna L, Taylor SA, Balakrishnan CN, Baldassarre DT, Deane-Coe PE et al. (2016) Genomic approaches to understanding population divergence and speciation in birds. Auk 133:13–30Article 

Google Scholar 
Toews DPL, Irwin DE (2008) Cryptic speciation in a Holarctic passerine revealed by genetic and bioacoustic analyses. Mol Ecol 17:2691–2705CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
van Els P, Cicero C, Klicka J (2012) High latitudes and high genetic diversity: Phylogeography of a widespread boreal bird, the gray jay (Perisoreus canadensis). Mol Phylogenet Evol 63:456–465PubMed 
Article 
PubMed Central 

Google Scholar 
Voelker G, Rohwer S (1998) Contrasts in scheduling of molt and migration in eastern and western warbling vireos. Auk 155:142–155Article 

Google Scholar 
Walsh J, Billerman SM, Rohwer VG, Butcher BG, Lovette IJ (2020) Genomic and plumage variation across the controversial Baltimore and Bullock’s oriole hybrid zone. Auk 137:1–15Article 

Google Scholar 
Walsh J, Rowe RJ, Olsen BJ, Shriver WG, Kovach AI (2016) Genotype-environment associations support a mosaic hybrid zone between two tidal marsh birds. Ecol Evol 6:279–294PubMed 
Article 
PubMed Central 

Google Scholar 
Walsh P, Metzger D, Higuchi R (1991) Chelex 100 as a medium for simple extraction of DNA for PCR-based typing from forensic material. Biotechniques 10:506–513CAS 
PubMed 
PubMed Central 

Google Scholar 
Weir JT, Schluter D (2004) Ice sheets promote speciation in boreal birds. Proc R Soc B 271:1881–1887PubMed 
PubMed Central 
Article 

Google Scholar 
Williams JW (2003) Variations in tree cover in North America since the last glacial maximum. Glob Planet Change 35:1–23Article 

Google Scholar 
Williams DA, Berg EC, Hale AM, Hughes CR (2004) Characterization of microsatellites for parentage studies of white-throated magpie-jays (Calocitta formosa) and brown jays (Cyanocorax morio). Mol Ecol Notes 4:509–511CAS 
Article 

Google Scholar 
Zwartjes PW (2001) Genetic structuring among migratory populations of the black-whiskered vireo, with a comparison to the red-eyed vireo. Condor 103:439–448Article 

Google Scholar  More

FennelThe main effect of fertilization (F) significantly impacted all measured parameters of fennel. Intercropping (I) pattern affected all parameters except plant height and 1000-seed weight. Significant I × F interactions occurred for umbel number, seed yield, essential oil content (EO), EO yield, oil content, and oil yield (Table 2).Table 2 Analysis of variance for the effect of cropping pattern and fertilization on evaluated traits in fennel.Full size tablePlant heightThe tallest fennel plants (125.8 cm) occurred in the sole cropping (Fs), while the shortest plants (98.6 cm) occurred in 1F:2FG. Across intercropping patterns, the Fs treatment had 28%, 14%, 11% taller fennel plants than 1F:2FG, 2F:2FG, and 2F:4FG, respectively (Fig. 1A). Compared to the unfertilized control, HA and BFS increased fennel plant height by 10% and 13%, respectively (Fig. 1B).Figure 1Means comparison for the main effects of cropping patterns [Fs (fennel sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on plant height (A), and fertilization [C (control), HA (humic acid), BFS (biofertilizers)] on plant height (B) and 1000-seed weight (C) of fennel. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size image1000-seed weightCompared to the unfertilized control (3.9 g per 1000 seeds), BFS and HA increased the 1000-seed weight of fennel by 24.1% and 14.5% (4.9 and 4.5 g per 1000 seeds), respectively (Fig. 1C).Umbel numberThe fennel sole cropping fertilized with HA produced the most umbels of fennel (51.5), while 2F:4FG without fertilization produced the least (32). Averaged across fertilizer types within each intercropping system, 1F:2FG, 2F:2FG, and 2F:4FG had 21.1%, 16.3%, and 26.7% fewer umbels than fennel sole cropping, respectively. Across intercropping systems, HA and BFS increased the umbel number by 17.8% and 16.5% compared with the unfertilized control, respectively (Fig. 2A).Figure 2Means comparison for the interaction effect of fertilization [C (control), HA (humic acid), BFS (biofertilizers)] and different cropping patterns [Fs (fennel sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on umbel number (A) and seed yield (B) of fennel. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size imageSeed yieldThe different intercropping patterns had lower fennel seed yields than fennel sole cropping. Sole cropping fertilized with BFS and HA produced the highest fennel seed yields (2233 and 2209 kg ha–1, respectively), followed by unfertilized sole cropping (1960 kg ha–1). The lowest seed yields occurred in the unfertilized controls in 1F:2FG (933 kg ha–1) and 2F:4FG (1033 kg ha–1). Averaged across fertilization treatments, fennel seed yield in 1F:2FG, 2F:2FG, and 2F:4FG decreased by 41.7, 26.8, and 36.3%, respectively, compared to fennel sole cropping (Fs). Averaged across intercropping patterns, HA and BFS increased fennel seed yield by 33.3% and 39.5% compared with the unfertilized control, respectively (Fig. 2B).Essential oil content and yieldThe different intercropping patterns produced higher EO contents of fennel than fennel sole cropping. The highest absolute EO content of fennel (4.22%) occurred in 2F:2FG fertilized with BFS, although this did not statistically differ from the 2F:2FG fertilized with HA (4.04%) or 2F:4FG fertilized with HA or BFS (3.8% and 4.00%, respectively) (Fig. 3A). The lowest EO contents occurred in the unfertilized control (2.38%), HA (2.55%), and BFS (2.57%) in the Fs system. Averaged across fertilization treatments, the EO content of fennel in 1F:2FG, 2F:2FG, and 2F:4FG increased by 36%, 52%, and 44% compared to fennel sole cropping, respectively. Within each intercropping pattern, and with the exception of Fs, the HA and BFS treatments had higher EO contents of fennel, none of which significantly differed, increasing by 25% and 29%, respectively (Fig. 3A).Figure 3Means comparison for the interaction effect of fertilization [C (control), HA (humic acid), BFS (biofertilizers)] and different cropping patterns [Fs (fennel sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on essential oil content (A), essential oil yield (B), oil content (C), and oil yield (D) of fennel. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size imageMaximum EO yields of fennel occurred with HA or BFS applied in 2F:2FG (65.2 and 66.6 kg ha–1) and 2F:4FG (60.7 and 65.5 kg ha–1), respectively, while the lowest EO yields occurred in the unfertilized control in 1F:2FG (27.2 kg ha–1), 2F:2FG (33.2 kg ha–1), and 2F:4FG (32.1 kg ha–1). Averaged across intercropping patterns, the EO yield of fennel increased by 66.1% and 74.7% with HA and BFS, respectively (Fig. 3B).Fennel essential oil compositionGC–FID and GC–MS analyses identified 14 components in the fennel EO (representing 97.4–99.9% of the total composition) (Table 3), with the main constituents being trans-anethole (78.3–84.85%), estragole (3.02–7.17%), fenchone (4.14–7.52%), and limonene (3.15–4.88%). The highest percentage of (E)-anethole, estragole, and fenchone occurred in 2F:2FG with BFS. The highest limonene content occurred in 2F:4 FG with HA. The relative contents of trans-anethole, fenchone, and limonene increased by 3.9%, 16.6%, and 8.4% compared with fennel sole cropping. Notably, the contents of most compounds increased with HA and BFS. Compared to the unfertilized control, trans-anethole, fenchone, and limonene contents increased by 2.9%, 21.5%, and 7.9% with BFS and 2.3%, 22.4%, and 11.9% with HA, respectively (Table 3).Table 3 Proportion of fennel essential oil constituents under different cropping patterns and fertilization.Full size tableFennel oil content and yieldAmong the studied treatments, the highest fennel oil content occurred with HA or BFS application in 1F:2FG (16.3% and 16.6%) and 2F:2FG (16.3% and 17.4%), respectively. The lowest fennel oil contents occurred in the unfertilized control, HA, and BFS treatments (12.5%, 12.8%, and 12.9%, respectively) under fennel sole cropping, and the unfertilized control in 2F:4FG (12.6%). Averaged across fertilizer treatments, fennel oil content in 1F:2FG, 2F:2FG, and 2F:4FG increased by 22.8%, 26.0%, and 12.6% compared with fennel sole cropping, respectively. Across intercropping patterns, HA and BFS increased fennel oil content by 13.5% and 16.5%, respectively (Fig. 3C).The maximum oil yield of fennel (318.6 kg ha–1) occurred in 2F:2FG fertilized with BFS, while the lowest oil yield (129.3 kg ha–1) occurred in 1F:2FG without fertilization. Across intercropping patterns, HA and BFS increased fennel oil yield by 50.8% and 62.6%, respectively (Fig. 3D).Oil compoundsGC–FID and GC–MS analyses identified nine constituents that represented 94.3–97.9% of the total fennel oil composition. The main oil constituents were oleic acid (39.2–48.3%), linoleic acid (17.1–24.8%), stearic acid (10.9–15.4%), lauric acid (10.1–14.00%), and arachidic acid (2.2–3.4%). The highest oleic and linoleic acid contents occurred in 2F:4FG and 2F:2FG fertilized with BFS, respectively. Across fertilizer treatments, oleic and linoleic acid contents increased by 6% and 21%, respectively, under different intercropping patterns compared with fennel sole cropping. Across systems, HA and BFS enhanced oleic acid content by 1.8% and 8% and linoleic acid by 7.9% and 8.2%, respectively, compared with the unfertilized control. The highest percentage of stearic and lauric acids occurred in the unfertilized control of fennel sole cropping. Conversely, the lowest stearic and lauric acid contents occurred in 2F:2FG and 2F:4FG fertilized with BFS, 16.1% and 14.2% higher than fennel sole cropping, respectively. Finally, HA and BFS decreased stearic acid content by an average of 5.4% and 7.2%, respectively (Table 4).Table 4 Proportion of fennel oil constituents under different cropping patterns and fertilization.Full size tablePhenolic compoundsThe main phenolic compounds of fennel were chlorogenic acid (10.4–15.3 ppm), quercetin (7.0–17.2 ppm), and cinnamic acid (4.1–8.9 ppm). The highest chlorogenic acid and quercetin contents occurred in 2F:2FG fertilized with BFS and HA, respectively, while the lowest contents occurred in the fennel sole cropping system without fertilizer. Averaged across the three intercropping patterns, the chlorogenic acid and quercetin contents were 18.5% and 80.1% higher than the fennel sole cropping system. The chlorogenic acid and quercetin contents increased by 13% and 17% with BFS and 22% and 15% with HA, respectively (Table 5).Table 5 Concentration of phenolic compounds in fennel under different cropping patterns and fertilization.Full size tableFenugreekThe main effects of intercropping (I) pattern (C) and fertilizer (F) were significant for all parameters analyzed in fenugreek. Significant I × F interactions occurred for plant height, pod number per plant, seed yield, oil content, and oil yield of fenugreek (Table 6).Table 6 Analysis of variance for the effects of cropping patterns and fertilization on evaluated traits in fenugreek.Full size tablePlant heightThe 2F:2FG intercropping system fertilized with BFS produced the tallest fenugreek plants (63 cm), followed by 1F:2FG with BFS (53.3 cm) and 2F:4FG with BFS (56 cm), and 2F:2FG with HA (55 cm). The unfertilized control produced the shortest fenugreek plants (42 cm) in the sole cropping. Most fertilizer treatments across different intercropping patterns produced taller fenugreek plants than their sole cropping counterparts. Across fertilizer treatments, 1F:2FG, 2F:2FG, and 2F:4FG produced 16.2%, 26.8%, and 14.6% taller fenugreek plants than sole cropping, respectively. Across cropping patterns, BFS and HA increased fenugreek plant height by 5.7% and 15.2% compared with the unfertilized control, respectively (Fig. 4A).Figure 4Means comparison for the interaction effect of fertilization [C (control), HA (humic acid), BFS (biofertilizers)] and different cropping patterns [FGs (fenugreek sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on plant height (A) and pod number per plant (B) of fenugreek. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size imagePod number per plantThe fenugreek sole cropping with BFS and HA and 2F:4FG with BFS produced the most pods per plant (21.3, 20.3, and 20, respectively), while the unfertilized controls in 1F:2FG, 2F:2FG, and 2F:4FG produced the least (11.6, 12, and 13.3, respectively). Across fertilization treatments, 1F:2FG, 2F:2FG, and 2F:4FG had 30.1%, 25.6%, and 14.3% fewer pods per plant, respectively, than the fenugreek sole cropping system. Across cropping systems, HA and BFS increased pod number per plant in fenugreek by 25% and 33%, respectively, relative to the corresponding sole cropping (Fig. 4B).Seed number per podAcross fertilization treatments, fenugreek sole cropping produced the most seeds per pod (7.09), followed by 2F:4FG (6.02), 2F:2FG (4.93), and 1F:2FG (4.41) (Fig. 5A). In relative terms, sole cropping produced 60.5%, 43.9%, and 17.6% more seeds per pod than 1F:2FG, 2F:2FG, and 2F:4FG (Fig. 5A). Across cropping patterns, BFS and HA increased seed number per pod by 8.1% and 17.4% compared with the unfertilized control, respectively (Fig. 5B).Figure 5Means comparison for the main effects of cropping patterns [FGs (fenugreek sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on seed number per pod (A) and 1000-seed weight (C), and fertilization [C (control), HA (humic acid), BFS (biofertilizers)] on seed number per pod (B) and 1000-seed weight (D) of fennel. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size image1000-seed weightAmong different cropping patterns, sole cropping and 1F:2FG produced the highest (10.45 g) and lowest (8.34 g) fenugreek seed weights, respectively. In relative terms, fenugreek sole cropping produced 25.3%, 21.8%, and 12.4% higher seed weights than 1F:2FG, 2F:2FG, and 2F:4FG, respectively (Fig. 5C). Across cropping patterns, BFS and HA increased fenugreek seed weight by 3.7% and 5.7% compared with the control, respectively (Fig. 5D).Seed yieldMeans comparisons showed that sole cropping produced higher fenugreek seed yields than intercropping patterns. Sole cropping with BFS (1240 kg ha–1) and HA (1217 kg ha–1) produced the highest seed yields followed by the unfertilized control (Fig. 6A). The unfertilized control in 1F:2FG (437 kg ha–1) and 2F:2FG (467 kg ha–1) produced the lowest fenugreek seed yields. In all cases, and within each cropping pattern, BFS and HS produced higher fenugreek seed yields than the unfertilized control. As a result, BFS and HA increased fenugreek seed yield by 25.2% and 31.5% compared with the unfertilized control, respectively (Fig. 6A).Figure 6Means comparison for the interaction effects of fertilization [C (control), HA (humic acid), BFS (biofertilizers)] and different cropping patterns [FGs (fenugreek sole cropping), 1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] on seed yield (A), oil content (B), and oil yield (C) of fenugreek. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size imageOil content and yieldThe 2F:2FG cropping pattern with BFS produced the highest fenugreek oil content (8.3%), while the unfertilized control in sole cropping produced the lowest (5.9%). Across fertilizer treatments, 1F:2FG, 2F:2 FG, and 2F:4 FG produced 11.7%, 18.5%, and 15.7% higher fenugreek oil contents than sole cropping, respectively. In the 2F:2FG and 2F:4FG cropping patterns, BFS produced higher oil content (%) than HA. As a result, across cropping patterns, HA and BFS increased fenugreek oil content by 12.3% and 19.4%, respectively (Fig. 6B).Sole cropping with HA and BFS and 2F:2FG with BFS produced the highest fenugreek oil yields (77.1, 80.0, and 74.4 kg ha–1, respectively), while the unfertilized controls in 1F:2FG and 2F:4FG produced the lowest (27.51 and 29.8 kg ha–1, respectively). The 1F:2FG, 2F:2FG, and 2F:4FG cropping patterns produced 45.9%, 20.7%, and 41.5% lower fenugreek oil yields than fenugreek sole cropping, respectively. Moreover, except for sole cropping, BFS produced the highest fenugreek oil yield, followed by HA and the unfertilized control (Fig. 6C).Oil compoundsGC–FID and GC–MS analyses identified seven constituents (representing 91.09–99.27% of the total composition) in fenugreek oil. The main oil constituents were linoleic acid (26.1–37.1%), linolenic acid (16.9–22.4%), oleic acid (15.1–21.2%), palmitic acid (11.2–17.1%), lauric acid (5.0–12.3%), and myristic acid (3.1–6.4%). The highest linoleic and oleic acid percentages occurred in 1F:2FG and 2F:4FG with BFS. The 1F:2FG cropping pattern with BFS also had the highest linolenic acid percentage. The fenugreek sole cropping system without fertilization (control) had the lowest content of these three compounds. The intercropping patterns had 17%, 18.2%, and 17.1% higher oleic, linoleic, and linolenic acid contents than fenugreek sole cropping. In addition, HA and BFS increased oleic acid content by 15.6% and 8.8%, linoleic acid content by 12.8% and 7%, and linolenic acid content by 7.5% and 12.9%, respectively. Fenugreek sole cropping without fertilization produced the highest lauric acid and palmitic contents, 29.33% and 22.81% higher than the intercropping patterns (Table 7).Table 7 Proportion of fenugreek oil constituents under different cropping patterns and fertilization.Full size tablePhenolic compoundsThe main phenolic compounds in fenugreek were chlorogenic acid (2.01–5.49 ppm), caffeic acid (2.42–4.93 ppm), quercetin (1.98–4.45 ppm), comaric (1.09–2.43 ppm), apigenin (1.97–2.99 ppm), and gallic acid (1.76–2.92 ppm). The 2F:2FG cropping pattern with HA produced the highest quercetin and gallic acid contents, and 2F:4FG with HA produced the highest chlorogenic and caffeic acid contents. The 2F:2FG and 2F:4FG cropping patterns with BFS produced the highest comaric and apigenin contents, respectively. In contrast, fenugreek sole cropping without fertilization produced the lowest contents of the abovementioned compounds (Table 8).Table 8 Proportion of fenugreek concentration of phenolic compounds under different cropping patterns and fertilization.Full size tableLand equivalent ratio (LER)The 2F:4FG and 2F:2FG intercropping patterns treated with BFS had the highest partial LERs for fennel (0.82) and fenugreek (0.72), respectively. In addition, 2F:2FG with BFS and 1F:2FG without fertilization produced the highest (1.42) and lowest (0.86) total LERs, respectively (Fig. 7).Figure 7Partial and total land equivalent ratio (LER) for seed yields of different fennel and fenugreek intercropping patterns [1F:2FG, 2F:2FG, 2F:4FG (ratios of fennel and fenugreek in the intercropping patterns)] and fertilization [C (Control), HA (humic acid), BFS (biofertilizers)]. Different lower-case letters above the bars indicate significant (p ≤ 0.05) differences.Full size image More

IPCC Climate Change 2013: The Physical Science Basis (eds Stocker, T. F. et al.) (Cambridge Univ. Press, 2013).Sindelar, A. J. et al. Winter oilseed production for biofuel in the US Corn Belt: opportunities and limitations. GCB Bioenergy 9, 508–524 (2017).CAS 

Google Scholar 
Stöckle, C. O. et al. Evaluating opportunities for an increased role of winter crops as adaptation to climate change in dryland cropping systems of the U.S. Inland Pacific Northwest. Clim. Change 146, 247–261 (2018).
Google Scholar 
Williams, C. M., Henry, H. A. L. & Sinclair, B. J. Cold truths: how winter drives responses of terrestrial organisms to climate change. Biol. Rev. 90, 214–235 (2015).
Google Scholar 
Seifert, C. A., Azzari, G. & Lobell, D. B. Satellite detection of cover crops and their effects on crop yield in the Midwestern United States. Environ. Res. Lett. 13, 064033 (2018).
Google Scholar 
Marcillo, G. S. & Miguez, F. E. Corn yield response to winter cover crops: an updated meta-analysis. J. Soil Water Conserv. 72, 226–239 (2017).
Google Scholar 
Zhu, L., Ives, A. R., Zhang, C., Guo, Y. & Radeloff, V. C. Climate change causes functionally colder winters for snow cover-dependent organisms. Nat. Clim. Change 9, 886–893 (2019).
Google Scholar 
Mankin, J. S. & Diffenbaugh, N. S. Influence of temperature and precipitation variability on near-term snow trends. Clim. Dynam. 45, 1099–1116 (2015).
Google Scholar 
Zhu, L., Radeloff, V. C. & Ives, A. R. Characterizing global patterns of frozen ground with and without snow cover using microwave and MODIS satellite data products. Remote Sens. Environ. 191, 168–178 (2017).
Google Scholar 
Huning, L. S. & AghaKouchak, A. Global snow drought hot spots and characteristics. Proc. Natl Acad. Sci. USA 117, 19753–19759 (2020).CAS 

Google Scholar 
Qin, Y. et al. Agricultural risks from changing snowmelt. Nat. Clim. Change 10, 459–465 (2020).
Google Scholar 
Trnka, M. et al. Adverse weather conditions for European wheat production will become more frequent with climate change. Nat. Clim. Change 4, 637–643 (2014).
Google Scholar 
Li, D., Wrzesien, M. L., Durand, M., Adam, J. & Lettenmaier, D. P. How much runoff originates as snow in the western United States, and how will that change in the future? Geophys. Res. Lett. 44, 6163–6172 (2017).
Google Scholar 
Biemans, H. et al. Importance of snow and glacier meltwater for agriculture on the Indo-Gangetic Plain. Nat. Sustain. 2, 594–601 (2019).
Google Scholar 
Acevedo, E., Silva, P. & Silva, H. in Bread Wheat: Improvement and Production (eds Curtis, B. C. et al.) 39–70 (FAO Plant Production and Protection, 2002).Baker, J. T., Pinter, P. J., Reginato, R. J. & Kanemasu, E. T. Effects of temperature on leaf appearance in spring and winter wheat cultivars. Agron. J. 78, 605–613 (1986).
Google Scholar 
Tack, J., Barkley, A. & Nalley, L. L. Effect of warming temperatures on US wheat yields. Proc. Natl Acad. Sci. USA 112, 6931–6936 (2015).CAS 

Google Scholar 
Müller, C. et al. Global gridded crop model evaluation: benchmarking, skills, deficiencies and implications. Geosci. Model Dev. 10, 1403–1422 (2017).
Google Scholar 
Talukder, A. S. M. H. M., McDonald, G. K. & Gill, G. S. Effect of short-term heat stress prior to flowering and early grain set on the grain yield of wheat. Field Crops Res. 160, 54–63 (2014).
Google Scholar 
Farooq, M., Bramley, H., Palta, J. A. & Siddique, K. H. M. Heat stress in wheat during reproductive and grain-filling phases. Crit. Rev. Plant Sci. 30, 491–507 (2011).Cuadra, S. V., Kimball, B. A., Boote, K. J., Suyker, A. E. & Pickering, N. Energy balance in the DSSAT-CSM-CROPGRO model. Agric. For. Meteorol. 297, 108241 (2021).
Google Scholar 
Harder, P., Helgason, W. D. & Pomeroy, J. W. Modeling the snowpack energy balance during melt under exposed crop stubble. J. Hydrometeorol. 19, 1191–1214 (2018).
Google Scholar 
Barlow, K. M., Christy, B. P., O’Leary, G. J., Riffkin, P. A. & Nuttall, J. G. Simulating the impact of extreme heat and frost events on wheat crop production: a review. Field Crops Res. 171, 109–119 (2015).
Google Scholar 
Wang, W. et al. Evaluation of air–soil temperature relationships simulated by land surface models during winter across the permafrost region. Cryosphere 10, 1721–1737 (2016).
Google Scholar 
Seifert, C. A. & Lobell, D. B. Response of double cropping suitability to climate change in the United States. Environ. Res. Lett. 10, 024002 (2015).
Google Scholar 
Pullens, J. W. M. et al. Risk factors for European winter oilseed rape production under climate change. Agric. For. Meteorol. 272–273, 30–39 (2019).
Google Scholar 
Chopra, R. et al. Identification and stacking of crucial traits required for the domestication of pennycress. Nat. Food 1, 84–91 (2020).
Google Scholar 
Crews, T. E., Carton, W. & Olsson, L. Is the future of agriculture perennial? Imperatives and opportunities to reinvent agriculture by shifting from annual monocultures to perennial polycultures. Glob. Sustain. 1, e11 (2018).Harkness, C. et al. Adverse weather conditions for UK wheat production under climate change. Agric. Meteorol. 282–283, 107862 (2020).
Google Scholar 
Schierhorn, F., Hofmann, M., Gagalyuk, T., Ostapchuk, I. & Müller, D. Machine learning reveals complex effects of climatic means and weather extremes on wheat yields during different plant developmental stages. Clim. Change 169, 39 (2021).Michel, S. et al. Improving and maintaining winter hardiness and frost tolerance in bread wheat by genomic selection. Front. Plant Sci. 10, 1195 (2019).
Google Scholar 
Mahfoozi, S., Limin, A. E. & Fowler, D. B. Influence of vernalization and photoperiod responses on cold hardiness in winter cereals. Crop Sci. 41, 1006–1011 (2001).
Google Scholar 
Dutra, E. et al. An improved snow scheme for the ECMWF land surface model: description and offline validation. J. Hydrometeorol. 11, 899–916 (2010).
Google Scholar 
Ge, Y. & Gong, G. Land surface insulation response to snow depth variability. J. Geophys. Res. Atmos. 115, 8107 (2010).
Google Scholar 
Hunt, J. R. et al. Early sowing systems can boost Australian wheat yields despite recent climate change. Nat. Clim. Change 9, 244–247 (2019).
Google Scholar 
Sloat, L. L. et al. Climate adaptation by crop migration. Nat. Commun. 11, 1243 (2020) .Ainsworth, E. A. & Long, S. P. 30 years of free-air carbon dioxide enrichment (FACE): what have we learned about future crop productivity and its potential for adaptation? Glob. Change Biol. 27, 27–49 (2021).
Google Scholar 
Shimoda, S. et al. Effects of snow compaction ‘yuki-fumi’ on soil frost depth and volunteer potato control in potato–wheat rotation system in Hokkaido. Plant Prod. Sci. 24, 186–197 (2021).CAS 

Google Scholar 
Luojus, K. et al. GlobSnow v3.0 Northern Hemisphere snow water equivalent dataset. Sci. Data 8, 163 (2021)..IMS Daily Northern Hemisphere Snow and Ice Analysis at 1 km, 4 km, and 24 km Resolutions Version 1 (NSIDC, 2008).Jing, Q. et al. Assessing the options to improve regional wheat yield in Eastern Canada using the CSM–CERES–wheat model. Agron. J. 109, 510–523 (2017).
Google Scholar 
Vogel, F. A. & Bange, G. A. Understanding USDA Crop Forecasts (USDA, 1999).Daly, C. et al. Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int. J. Climatol. 28, 2031–2064 (2008).
Google Scholar 
Brown, R. D. & Brasnett, B. Daily Snow Depth Analysis Data Version 1 (Canadian Meteorological Centre, 2010).Brasnett, B. A global analysis of snow depth for numerical weather prediction. J. Appl. Meteorol. Climatol. 38, 726–740 (1999).
Google Scholar 
Toure, A. M., Reichle, R. H., Forman, B. A., Getirana, A. & De Lannoy, G. J. M. Assimilation of MODIS snow cover fraction observations into the NASA catchment land surface model. Remote Sens. 10, 316 (2018).
Google Scholar 
Snauffer, A. M., Hsieh, W. W. & Cannon, A. J. Comparison of gridded snow water equivalent products with in situ measurements in British Columbia, Canada. J. Hydrol. 541, 714–726 (2016).
Google Scholar 
Census of Agriculture (USDA National Agricultural Statistics Service, 2017).Skinner, D. Z. & Mackey, B. Freezing tolerance of winter wheat plants frozen in saturated soil. Field Crops Res. 113, 335–341 (2009).
Google Scholar 
Lollato, R. P. et al. Climate-risk assessment for winter wheat using long-term weather data. Agron. J. 112, 2132–2151 (2020).
Google Scholar 
Siebers, M. H. et al. Heat waves imposed during early pod development in soybean (Glycine max) cause significant yield loss despite a rapid recovery from oxidative stress. Glob. Change Biol. 21, 3114–3125 (2015).
Google Scholar 
Çakir, R. Effect of water stress at different development stages on vegetative and reproductive growth of corn. Field Crops Res. 89, 1–16 (2004).
Google Scholar 
Lobell, D. B. et al. The critical role of extreme heat for maize production in the United States. Nat. Clim. Change 3, 497–501 (2013).
Google Scholar 
Chen, M., Griffis, T. J., Baker, J., Wood, J. D. & Xiao, K. Simulating crop phenology in the Community Land Model and its impact on energy and carbon fluxes. J. Geophys. Res. Biogeosci. 120, 310–325 (2015).CAS 

Google Scholar 
Larson, K. M. & Small, E. E. Daily Snow Depth and SWE from GPS Signal-to-Noise Ratios Version 1 (NSIDC, 2017).Sturm, M. et al. Estimating snow water equivalent using snow depth data and climate classes. J. Hydrometeorol. 11, 1380–1394 (2010).
Google Scholar 
McCabe, G. J. & Wolock, D. M. Recent declines in western U.S. snowpack in the context of twentieth-century climate variability. Earth Interact. 13, 1–15 (2009).
Google Scholar 
Wu, X. et al. Uneven winter snow influence on tree growth across temperate China. Glob. Change Biol. 25, 144–154 (2019).
Google Scholar 
Qiao, S. et al. Robust negative impacts of climate change on African agriculture. Environ. Res. Lett. 5, 014010 (2010).
Google Scholar 
Lobell, D. B., Schlenker, W. & Costa-Roberts, J. Climate trends and global crop production since 1980. Science 333, 616–620 (2011).CAS 

Google Scholar 
Xie, Y., Gibbs, H. K. & Lark, T. J. Landsat-based Irrigation Dataset (LANID): 30 m resolution maps of irrigation distribution, frequency, and change for the US, 1997–2017. Earth Syst. Sci. Data 13, 5689–5710 (2021).
Google Scholar 
Mueller, N. D. et al. Closing yield gaps through nutrient and water management. Nature 490, 254–257 (2012).CAS 

Google Scholar 
Friedman, J., Hastie, T. & Tibshirani, R. Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33, 1–22 (2010).
Google Scholar 
Elliott, J. et al. The global gridded crop model intercomparison: data and modeling protocols for phase 1 (v1.0). Geosci. Model Dev. 8, 261–277 (2015).
Google Scholar 
Li, X., Shen, Z., Harri, A. & Coble, K. H. Comparing survey-based and programme-based yield data: implications for the U.S. Agricultural Risk Coverage-County programme. Geneva Pap. Risk Insur. Issues Pract. 45, 184–202 (2020).
Google Scholar 
Hawkins, E., Osborne, T. M., Ho, C. K. & Challinor, A. J. Calibration and bias correction of climate projections for crop modelling: an idealised case study over Europe. Agric. Meteorol. 170, 19–31 (2013).
Google Scholar 
Ho, C. K., Stephenson, D. B., Collins, M., Ferro, C. A. T. & Brown, S. J. Calibration strategies: a source of additional uncertainty in climate change projections. Bull. Am. Meteorol. Soc. 93, 21–26 (2012).
Google Scholar  More

Previous studies have used number of species threats6,7, countryside species-area relationship1,3,17, and potentially disappeared fraction of species4 to quantify biodiversity loss. We introduce the non-normalised Species Threat Abatement and Restoration (nSTAR) metric as the quantifiable representation of biodiversity loss in our analysis, a unit-less, species-centred metric which relies on detailed information curated in the IUCN Red List of Threatened Species11. On its own, this metric can be used to support production-based accounting of the extinction risk of species and identify the most significant threats at a specific location to inform direct interventions26. However, once manipulated into a structure that allows it to be appended to a multi-region input–output (MRIO) table, an environmentally-extended MRIO can be created. This unlocks the power of consumption-based accounting of this extinction risk, connecting the direct environmental impact with the consumption which ultimately induces it.IUCN Red List of Threatened SpeciesThe IUCN Red List version 2020–211 provided information on extinction risk for over 122,000 species and details of the threats acting on those species, including the threat classification, scope, timing, and severity. The species scope was limited to comprehensively assessed terrestrial species, ensuring that only species which have been assessed across all countries were included, and thus eliminating any geographical bias introduced by incomplete assessments27. Species with an extinction risk category of Near Threatened (NT), Vulnerable (VU), Endangered (EN), or Critically Endangered (CR) were included. Three species were excluded to avoid double counting where two different extinction risk categories were provided for the same species, leaving 5295 amphibian, mammal, and bird species in scope.The information contained in the IUCN Red List regarding the threats facing each species is crucial, since many of these threats are attributable to economic activity28,29. Specialist assessors are required to assign one or more of 118 different threat classes to each species’ record, with additional documentation of the severity, scope and timing of each threat recommended, based on the impact of that threat on the species’ population30. To connect this threat information to economic sectors, a key requirement for input–output analysis, background information on threat classes was sourced from the IUCN Threats Classification Scheme version 3.229. Each threat was assessed for connection to each of the 6357 economic sectors classified in the UN Statistics Division Central Product Classification Standard31, based on the likelihood that activity associated with each sector directly contributes to the threat being assessed. As an example, the economic sectors associated with rice cultivation were allocated to the threats grouped under IUCN Threat Class 2.1—Annual & perennial non-timber crops. A total of 55 out of 118 threats were allocated to at least one economic sector, with higher-level threat classes excluded from this allocation if information was available for the associated lower-level threat classes to avoid double counting. Species threats driven by activity that cannot be attributed to an economic sector, such as invasive species, were not allocated to any sectors and as a result, the extinction-risk footprint does not necessarily represent the full magnitude of extinction risk for each species. While not all threats were allocated to an economic sector, all economic sectors were allocated to at least one threat. Further details on the connection of economic sectors to threats are available in Supplementary Note S5, which includes a link to the detailed 6357 × 118 binary concordance matrix used to execute these sector-threat allocations.The IUCN Red List also requires inclusion of a range map and habitat classification, which were combined with remote sensed land cover and elevation data to generate a high-resolution area of habitat (AOH) map for each in-scope species32,33. These maps, reapplied from Strassburg et al.34, were used to calculate the percentage of each species’ AOH present in each country.Quantifying biodiversity loss: the nSTAR metricThis detailed information from the IUCN Red List was used to calculate the nSTAR metric, which quantifies each threat’s impact, rather than just its presence, on each species. Adapted from the newly developed Species Threat Abatement and Restoration metric (STAR)26 by removing the normalisation step, the nSTAR metric, which has no units, was calculated for each species in two stages.First, a numeric representation of each species’ extinction risk category (Wi) was determined, following the equal steps methodology introduced by Butchart et al.35. Extinction risk categories of Data Deficient (DD) and Least Concern (LC) were assigned Wi = 0, Near Threatened (NT) was assigned Wi = 1, Vulnerable (VU) was assigned Wi = 2, Endangered (EN) was assigned Wi = 3, and Critically Endangered (CR) was assigned Wi = 4.Next, a Threat Impact score (TSij) for each threat (j) acting on a species (i) was determined based on the scope and severity information recorded for that threat, according to the values set out in Table 1, which are adapted from those proposed by Garnett et al.36. Reapplying the methodology of the STAR metric, where no value was recorded for the scope or severity of a threat, the median possible value for these were used, and only threats noted as Ongoing or Future were included. Further details on these methodological choices and sensitivity analyses to support them are available in Mair et al.26.Table 1 Numeric representation of threat information.Full size tableThe numeric nSTAR value for each species-threat combination (ij) was calculated by multiplying the value representing the species’ extinction risk category (Wi) by the Threat Impact score (TSij) for that threat:$${text{nSTAR}}_{ij} = W_{i} *TS_{ij}$$
(1)
The total nSTAR for species (i) can be calculated by multiplying the extinction risk category value (Wi) for that species by the sum of all Threat Impact scores for the species:$${text{nSTAR}}_{i} = W_{i} *(TS_{i1} + TS_{i2} + TS_{i3} + cdots + TS_{ij} )$$
(2)
Once calculated according to Eq. (1), the nSTARij value for each species-threat combination was allocated to economic sectors using the 6357 × 118 sector-threat concordance (available in Supplementary Note S5), which was normalised based on the economic size of each sector. Finally these nSTAR values, derived for each species-sector combination, were allocated to each country based on the country’s share of the AOH for that species, calculated from the intersection of the species’ AOH map with each country’s borders34.The nSTAR metric introduced here differs from the STAR metric from which it is adapted in that the normalisation step executed at this point in the STAR methodology is omitted. This ensures that the nSTAR metric is both additive and independent across all three dimensions of species, country, and economic sector, a necessary condition for use in input–output analysis. The STAR metric normalises the total value calculated in Eq. (2) to ensure that the total STAR value for any species is equal to Wi * 100, resulting in all species with the same extinction risk category being allocated the same STAR value regardless of the number of threats acting on them26. This normalisation facilitates the aggregation of the STAR metric by species taxonomy however it is problematic when aggregating the STAR metric by threat, since the STAR value attributed to each species-threat combination will be dependent not only on the characteristics of that threat, but also on the number and characteristics of other threats acting on the species. This dependence on more than one variable in the calculation of the STAR value for each species-threat combination means that it is not suitable for aggregation by threat and, by extension, economic sectors once the threat to sector allocation has been carried out.In order to provide a metric which can be aggregated and disaggregated across species, sector, and country hierarchies the nSTAR methodology excludes this normalisation step. Consistent with the STAR methodology, the nSTAR metric is calculated using numeric values only and therefore has no unit of measure26.Input–output analysisOnce calculated, the nSTAR metric was partnered with the global supply-chain data available in the 2013 Eora MRIO, chosen for its extensive coverage of 190 regions (189 countries and one ‘rest of world’ region) and between 26 and 1022 economic sectors in each country, depending on the level of detail in each country’s publicly available National Accounts12.A satellite block, or Q matrix, was created using the nSTAR values for 5295 species across 6357 economic sectors for 190 regions. This satellite block was then aggregated to match the sectoral structure of the Eora MRIO, a total of 14,839 country-sector combinations. A process flow diagram to illustrate the stages of data manipulation required to convert the IUCN Red List data to a satellite block ready for use with the Eora MRIO is included in Supplementary Fig. S5.The Eora MRIO provided the intermediate transaction matrix T, the final demand matrix Y, and the value-added matrix V. Consumption based footprints were calculated by connecting the nSTAR value captured in the satellite block Q to the final demand matrix Y following Leontief’s methodology9,10. Central to this methodology is the Leontief Inverse L, a concise mathematical representation of the interdependencies across all economic sectors, which is expressed as:$${mathbf{L}} = left( {{mathbf{I}}{-}{mathbf{A}}} right)^{{ – {1}}}$$
(3)
where I is an identity matrix with dimensions equal to the those of the intermediate transaction matrix T, and A is the direct requirements matrix, derived from the T matrix in a number of stages. First the total output vector x is calculated, then diagonalised, and the inverse calculated to derive ({widehat{mathbf{X}}}^{-1}), which returns the direct requirements matrix A when multiplied by T.Next the satellite block was converted into an intensity matrix q by multiplying Q by ({widehat{mathbf{X}}}^{-1}) to calculate the nSTAR value attributable to each dollar of total output produced by each sector. Once the q, L and Y matrices are available, the consumption extinction-risk footprint for a sector k (fk) can be calculated using Eq. (4):$${mathbf{f}}_{k} = {mathbf{q}}*{mathbf{L}}*{mathbf{Y}}_{k}$$
(4)
where Yk represents the final demand for that sector. Consumption extinction-risk footprint values were generated for each species-sector-country combination, a total of more than 78 million datapoints.Further matrix manipulation was used to calculate the country-level imported, exported, and domestic extinction-risk footprints. For each country the final demand matrix, Y, was separated into two matrices, Ydom, representing demand from that country for the economic sectors in that country, and Yoth, representing demand from all other countries for the economic sectors in that country. Next, the intensity matrix, q, was separated into two matrices, qdom, representing the nSTAR intensity for each of the species within that country’s borders, and qoth, representing the nSTAR intensity for all remaining species. The three sub-footprints for each country were calculated using Eqs. (5), (6) & (7). A simplified illustration of this methodology is included in Supplementary Fig. S3.$${mathbf{f}}_{dom} = {mathbf{q}}_{dom} *{mathbf{L}}*{mathbf{Y}}_{dom}$$
(5)
$${mathbf{f}}_{exp} = {mathbf{q}}_{dom} *{mathbf{L}}*{mathbf{Y}}_{oth}$$
(6)
$${mathbf{f}}_{imp} = {mathbf{q}}_{oth} *{mathbf{L}}*{mathbf{Y}}_{dom}$$
(7)
Imported, exported, and domestic extinction-risk footprints were calculated for 188 countries.LimitationsWhile very powerful in unravelling the intricacies of the global economy, there are limitations to the effectiveness of input–output analysis. Since it relies on National Accounts data, only activity which can be directly connected into reported economic activity is captured. This means that any activities that are not transacted within the boundaries of the formal economy, such as subsistence hunting and illegal logging, will be excluded unless they have been incorporated into the relevant country’s National Accounts data. The exclusion of threats due to their timing or non-economic classification (such as geological events, disease, and invasive species) resulted in a zero nSTAR value for 519 species, leaving 4776 species with a material nSTAR value. In addition, any limitations in the sector categorisations, their spatial and technological homogeneity, or assumptions included in the allocation of economic activity to sectors within the National Accounts data in each country will be propagated through to the footprint calculations. These limitations are common to consumption-based analyses5,6,7,17,25 and we do not further address them here.Further limitations exist with the use of the scope and severity data for each threat captured in the IUCN Red List, since this does not take into account interaction between threats, or between the severity and scope of an individual threat36. As a result, the impact from a single threat acting on a species may be overstated, and higher nSTAR values attributed to that species than would otherwise be warranted. In addition, any variations in the location, scope, or severity of threats acting across a species’ distribution range are not captured and thus the impact of different economic sectors may be over or under-represented26.There is a temporal displacement between the economic activity and the extinction risk used in this analysis. The extinction risk category assigned to each species is due to the cumulative sum of current and historical impacts acting on it, while the value of economic interactions used to trace this extinction risk through the global economy is based on one year of activity. This is typical of related approaches1,6, and may not introduce much uncertainty given that current economic activity is higher than at any time in history37. Nevertheless, there is no doubt that some current extinction risk is due to past economic activity and development of methods to incorporate this temporal dimension would be a valuable research avenue. More

Butzer, K. W. Early Hydraulic Civilization in Egypt: a Study in Cultural Ecology (University of Chicago Press, Chicago, 1976).Said, R. The River Nile: Geology, Hydrology and Utilization (Pergamon Press, Oxford, 1993).Zeder, M. A. Domestication and early agriculture in the Mediterranean Basin: origins, diffusion, and impact. Proc. Natl Acad. Sci. USA 105, 11597–11604 (2008).CAS 
Article 

Google Scholar 
Shirai, N. The Archaeology of the First Farmer-Herders in Egypt: New Insights into the Fayum Epipalaeolithic and Neolithic (Uni. Leiden Press, the Netherlands, 2010).Garcea, E. A. A. Multi-stage dispersal of Southwest Asian domestic livestock and the path of pastoralism in the Middle Nile Valley. Quat. Int. 412, 54–64 (2016).Article 

Google Scholar 
Wilson, P. Prehistoric settlement in the western Delta: a regional and local view from Sais (Sa el-Hagar). J. Egypt. Archaeol. 92, 75–126 (2006).Article 

Google Scholar 
Van Geel, B. Non-Pollen Palynomorphs. Smol J. P., Birks H. J. B., Last W. M., Bradley R. S., Alverson K. (eds) Tracking Environmental Change Using Lake Sediments. Developments in Paleoenvironmental Research, 3 (Springer, Dordrecht, 2002).Van Geel, B., Hallewas, J. P. & Pals, J. P. A Late Holocene deposit under the Westfriese Zeedijk near Nkhuizen (Prov. of N-Holland, The Netherlands): palaeoecological and archaeological aspects. Rev. Palaeobot. Palyno 38, 269–335 (1983).Article 

Google Scholar 
Van Geel, B. A paleoecological study of Holocene peat bog sections in Germany and the Netherlands. Rev. Palaeobot. Palyno 25, 1–120 (1978).Article 

Google Scholar 
Marinova, E. & Atanassova, J. Anthropogenic impact on vegetation and environment during the bronze age in the area of Lake Durankulak, NE Bulgaria: pollen, microscopic charcoal, non-pollen palynomorphs and plant macrofossils. Rev. Palaeobot. Palyno. 141, 165–178 (2006).Article 

Google Scholar 
Van Geel, B. et al. Diversity and ecology of tropical African fungal spores from a 25,000-year palaeoenvironmental record in southeastern Kenya. Rev. Palaeobot. Palynol. 164, 174–190 (2011).Article 

Google Scholar 
Gelorini, V., Verbeken, A., van Geel, B. B., Cocquyt, C. & Verschuren, D. Modern non-pollen palynomorphs from East African lake sediments. Rev. Palaeobot. Palyno 164, 143–173 (2011).Article 

Google Scholar 
Hillbrand, M., Geel, B. V., Hasenfratz, A., Hadorn, P. & Haas, J. N. Non-pollen palynomorphs show human-and livestock-induced eutrophication of Lake Nussbaumersee (Thurgau, Switzerland) since Neolithic times (3840 BC). Holocene 24, 559–568 (2014).Article 

Google Scholar 
Stanley, J. D. & Warne, A. G. Sea level and initiation of Predynastic culture in the Nile delta. Nature 363, 435–438 (1993).Article 

Google Scholar 
Pennington, B. T., Sturt, F., Wilson, P., Rowland, J. & Brown, A. G. The fluvial evolution of the Holocene Nile Delta. Quarter. Sci. Rev 170, 212–231 (2017).Article 

Google Scholar 
Negm, A. M., Saavedra O., & El-Adawy A. In The Handbook of Environmental Chemistry, 55 (Springer, 2017).Viste, E. & Sorteberg, A. The effect of moisture transport variability on Ethiopian summer precipitation. Int. J. Climatol. 33, 3106–3123 (2013).Article 

Google Scholar 
Revel, M., Colin, C., Bernasconi, S., Combourieu-Nebout, N. & Mascle, J. 21,000 years of Ethiopian African monsoon variability recorded in sediments of the western Nile deep-sea fan. Reg. Environ. Change 14, 1685–1696 (2014).Article 

Google Scholar 
Wijmstra, T. A., Smit, A., Van der Hammen, T. & Van Geel, B. Vegetational succession, fungal spores and short-term cycles in pollen diagrams from the Wietmarscher Moor. Acta Botanica Neerlandica 20, 401–410 (1971).Article 

Google Scholar 
Wilson, P. In The Nile Delta as a centre of cultural interactions between Upper Egypt and the Southern Levant in the 4th millennium BC, 299–318 (Poznań Archaeological Museum, Poznan, 2014).Zong, Y. Q. et al. Fire and flood management of coastal swamp enabled first rice paddy cultivation in east China. Nature 449, 459–462 (2007).CAS 
Article 

Google Scholar 
Yang, S. et al. Modern pollen assemblages from cultivated rice fields and rice pollen morphology: application to a study of ancient land use and agriculture in the Pearl River delta, China. The Holocene 22, 1393–1404 (2012).Article 

Google Scholar 
He, K. et al. Middle-Holocene sea-level fluctuations interrupted the developing Hemudu Culture in the lower Yangtze River. China. Quarter. Sci. Rev. 188, 90–103 (2018).Article 

Google Scholar 
Edwards, K. J., Whittington, G., Robinson, M. & Richter, D. Palaeoenvironments, the archaeological record and cereal pollen detection at Clickimin, Shetland, Scotland. J. Archaeo. Sci. 32, 1741–1756 (2005).Article 

Google Scholar 
Andersen, S. T. Identification of Wild Grass and Cereal Pollen [fossil Pollen, Annulus Diameter, Surface Sculpturing], Aarbog, 69–92 (Danmarks Geologiske Undersoegelse, 1979).Tweddle, J. C., Edwards, K. J. & Fieller, N. R. Multivariate statistical and other approaches for the separation of cereal from wild Poaceae pollen using a large Holocene dataset. Veg. Hist. Archaeobot. 14, 15–30 (2005).Article 

Google Scholar 
Joly, C., Barille, L., Barreau, M., Mancheron, A. & Visset, L. Grain and annulus diameter as criteria for distinguishing pollen grains of cereals from wild grasses. Rev. Palaeobot. Palynol. 146, 221–233 (2007).Article 

Google Scholar 
Salgado-Labouriau, M. L. & Rinaldi, M. Palynology of Gramineae of the Venezuelan mountains. Grana Palynologica 29, 119–128 (1990).Article 

Google Scholar 
Josefsson, T., Ramqvist, P. H. & Rnberg, G. The history of early cereal cultivation in northernmost Fennoscandia as indicated by palynological research. Veg. Hist. Archaeobot. 23, 821–840 (2014).Article 

Google Scholar 
Zhao, X. S. et al. Climate-driven early agricultural origins and development in the Nile Delta. Egypt. J. Archaeo. Sci. 136, 105498 (2021).Article 

Google Scholar 
Willcox, G. The distribution, natural habitats and availability of wild cereals in relation to their domestication in the near east: multiple events, multiple centres. Veg. Hist. Archaeobot. 14, 534–541 (2005).Article 

Google Scholar 
Riemer, H. Barbara e. barich. People, water and grain: the beginnings of domestication in the Sahara and the Nile Valley, Roma 1998. Archol. Inf. 24, 117–119 (2014).
Google Scholar 
Arranz-Otaegui, A., Colledge, S., Zapata, L., Teira-Mayolini, L. C. & Juan, J. Regional diversity on the timing for the initial appearance of cereal cultivation and domestication in Southwest Asia. Proc. Natl Acad. Sci. USA 113, 201612797 (2016).Article 

Google Scholar 
Zohary, D., Hopf, M. & Weiss, E. Domestication of plants in the Old World (Oxford University Press, Oxford, 2012).Kvavadze, E. & Bitadze, N. L. Special issue: fresh insights into the palaeoecological and palaeoclimatological value of quaternary non-pollen palynomorphs || Fibres of Linum (flax), Gossypium (cotton) and animal wool as non-pollen palynomorphs in the Late Bronze Age burials of Saphar-Kharaba, southern Georgia. Veg. Hist. Archaeobot. 19, 479–494 (2010).Article 

Google Scholar 
Karg, S. New research on the cultural history of the useful plant Linum usitatissimum L. (flax), a resource for food and textiles for 8,000 years. Veg. Hist. Archaeobot. 20, 507–508 (2011).Article 

Google Scholar 
Zhao, X. S. et al. Holocene climate change and its influence on early agriculture in the Nile Delta, Egypt. Palaeogeogr. Palaeoclimatol. Palaeoecol. 547, 109702 (2020).Article 

Google Scholar 
Reimer, P. et al. The IntCal20 northern hemisphere radiocarbon age calibration curve (0–55 cal kBP). Radiocarbon 62, 725–757 (2020).CAS 
Article 

Google Scholar 
Blaauw, M. & Christen, J. A. Flexible paleoclimate age-depth models using an autoregressive gamma process. Bayesian Analysis. 6, 457–474 (2011).Article 

Google Scholar 
Moore, P. D., Webb, J. A. & Collison, M. E. Pollen analysis (Blackwell Scientific Publications, Oxford, UK, 1991).Kholeif, S. E. A. & Mudie, P. J. Palynological records of climate and oceanic conditions in the Late Pleistocene and Holocene of the Nile Cone, Southeastern Mediterranean, Egypt. Palynology 33, 1–24 (2009).Article 

Google Scholar 
Leroy, S. A. G. Palynological evidence of Azolla nilotica Dec. in recent Holocene of the eastern Nile Delta and palaeoenvironment. Veg. Hist. Archaeobot. 1, 43–52 (1992).Article 

Google Scholar 
Kholeif, S. E. A. Holocene paleoenvironmental change in inner continental shelf sediments, Southeastern Mediterranean, Egypt. J. Afr. Earth. Sci. 57, 143–153 (2010).CAS 
Article 

Google Scholar  More

The MiDAS global consortium was established in 2018 to coordinate the sampling and collection of metadata from WWTPs across the globe (Supplementary Data 1). Samples were obtained in duplicates from 740 WWTPs in 425 cities, 31 countries on six continents (Fig. 1a). The majority of the WWTPs were configured with the activated sludge process (69.7%) (Fig. 1b), and these were the main focus of the subsequent analyses. Nevertheless, WWTPs based on biofilters, moving bed bioreactors (MBBR), membrane bioreactors (MBR), and granular sludge were also sampled to cover the microbial diversity in other types of WWTPs. The activated sludge plants were designed for carbon removal only (C; 22.1%), carbon removal with nitrification (C,N; 9.5%), carbon removal with nitrification and denitrification (C,N,DN; 40.9%), and carbon removal with nitrogen removal and enhanced biological phosphorus removal, EBPR (C,N,DN,P; 21.7%) (Fig. 1c). The first type represents the simplest design whereas the latter represents the most advanced process type with varying oxic and anoxic stages or compartments.Fig. 1: Sampling of WWTPs across the world.a Geographical distribution of WWTPs included in the study and their process configuration. b Distribution of plant types. MBBR moving bed bioreactor, MBR membrane bioreactor. c Distribution of process types for the activated sludge plants. C carbon removal, C,N carbon removal with nitrification, C,N,DN carbon removal with nitrification and denitrification, C,N,DN,P carbon removal with nitrogen removal and enhanced biological phosphorus removal (EBPR). The values next to the bars are the number of WWTPs in each group.Full size imageMiDAS 4: a global 16S rRNA gene catalogue and taxonomy for WWTPsMicrobial community profiling at high taxonomic resolution (genus- and species-level) using 16S rRNA gene amplicon sequencing requires a reference database with high-identity reference sequences (≥99% sequence identity) for the majority of the bacteria in the samples and a complete seven-rank taxonomy (domain to species) for all reference sequences16,20. To create such a database for bacteria in WWTPs globally, we applied synthetic long-read full-length 16S rRNA gene sequencing20,21 on samples from all WWTPs included in this study.More than 5.2 million full-length 16S rRNA gene sequences were obtained after quality filtering and primer trimming. The sequences were processed with AutoTax20 to yield 80,557 full-length 16S rRNA gene amplicon sequence variant (FL-ASVs). These reference sequences were added to our previous MiDAS 3 database16, providing a combined database (MiDAS 4) with a total of 90,164 unique, chimera-free FL-ASV reference sequences. The absence of detectable chimeric sequences is a unique feature of the database and is achieved due to the attachment of unique molecular identifiers (UMIs) to each end of the original template molecules before any PCR amplification steps21. This allows filtering of true biological sequences from chimera already in the synthetic long-read assembly20,21. The novelty of the FL-ASVs were determined based on the percent identity shared with their closest relatives in the SILVA 138 SSURef NR99 database and the threshold for each taxonomic rank proposed by Yarza et al.22. Out of all FL-ASVs, 88% had relatives above the genus-level threshold (≥94.5% identity) and 56% above the species-level threshold (≥98.7% identity) (Fig. 2 and Table 1).Fig. 2: Novel sequences and de novo taxa defined in the MiDAS 4 reference database.The phylogenetic trees are based on a multiple alignment of all MiDAS 4 reference sequences, which were first aligned against the global SILVA 138 alignment using the SINA aligner, and subsequently pruned according to the ssuref:bacteria positional variability by parsimony filter in ARB to remove hypervariable regions. The eight phyla with most FL-ASVs are highlighted in different colours. Sequence novelty was determined by the percent identity between each FL-ASV and their closest relative in the SILVA_138_SSURef_Nr99 database according to Usearch mapping and the taxonomic thresholds proposed by Yarza et al.22 shown in Table 1. Taxonomy novelty was defined based on the assignment of de novo taxa by AutoTax20.Full size imageTable 1 Novel sequences and de novo taxa observed in the MiDAS 4 reference database.Full size tableMiDAS 4 provides placeholder names for many environmental taxaAlthough only a small percentage of the reference sequences in MiDAS 4 represented new putative taxa at higher ranks (phylum, class, or order) according to the sequence identity thresholds proposed by Yarza et al.22, a large number of sequences lacked lower-rank taxonomic classifications and was assigned de novo placeholder names by AutoTax20 (Fig. 2 and Table 1). In total, de novo taxonomic names were generated by AutoTax for 26 phyla (30.6% of observed), 83 classes (37.2% of observed), 297 orders (46.8% of observed), and more than 8000 genera (86.3% of observed). Without the de novo taxonomy we would not be able to discuss these taxa across studies to unveil their potential role in wastewater treatment systems.Phylum-specific phylogenetic trees were created to determine if the FL-ASV reference sequences that were assigned to de novo phyla were actual phyla or simply artifacts related to the naive sequence identity-based assignment of de novo placeholder taxonomies (Supplementary Fig. 1a). The majority (65 FL-ASVs) created deep branches from within the Alphaproteobacteria together with 16S rRNA gene sequences from mitochondria, suggesting they represented divergent mitochondrial genes rather than true novel phyla. We also observed several FL-ASVs assigned to de novo phyla that branched from the classes Parcubacteria (3 FL-ASVs) and Microgenomatis (22 FL-ASVs) within the Patescibacteria phylum. These two classes were originally proposed as superphyla due to an unusually high rate of evolution of their 16S rRNA genes23,24. It is, therefore, likely that these de novo phyla are also artefacts due to the simple taxonomy assignment approach, which does not take different evolutionary rates into account20. Most of the class- and order-level novelty was found within the Patescibacteria, Proteobacteria, Firmicutes, Planctomycetota and Verrucomicrobiota. (Supplementary Fig. 1b). At the family- and genus-level, we also observed many de novo taxa affiliated to Bacteroidota, Bdellovibrionota and Chloroflexi.MiDAS 4 provides a common taxonomy for the fieldThe performance of the MiDAS 4 database was evaluated based on an independent amplicon dataset from the Global Water Microbiome Consortium (GWMC) project2, which covers ~1200 samples from 269 WWTPs. The raw GWMC amplicon data of the 16S rRNA gene V4 region was resolved into ASVs, and the percent identity to their best hits in MiDAS 4 and other reference databases was calculated (Fig. 3). The MiDAS 4 database had high-identity hits (≥99% identity) for 72.0 ± 9.5% (mean ± SD) of GWMC ASVs with ≥0.01% relative abundance, compared to 57.9 ± 8.5% for the SILVA 138 SSURef NR99 database, which was the best of the universal reference databases (Fig. 3). The relative abundance cutoff selects taxa that likely have a quantitative impact on the ecosystem while filtering out the rare biosphere which includes many bacteria introduced with the influent wastewaters25. Similar analyses of ASVs obtained from the samples included in this study showed, not surprisingly, even better performance with high-identity hits for 90.7 ± 7.9% of V1–V3 ASVs and 90.0 ± 6.6% of V4 ASVs with ≥0.01% relative abundance, compared to 60.6 ± 11.9% and 73.9 ± 10.3% for SILVA (Supplementary Fig. 2a). Although the sampling of WWTPs was focused towards activated sludge plants, the MiDAS 4 database also includes high-identity references for most ASVs in other plant types (granules, biofilters, etc.) (Supplementary Fig. 2b). This suggests that most taxa were shared across plant types, although often present in other relative abundances.Fig. 3: Database evaluation based on amplicon data from the Global Water Microbiome Consortium project.Raw amplicon data from the Global Water Microbiome Consortium project2 was processed to resolve ASVs of the 16S rRNA gene V4 region. The ASVs for each of the samples were filtered based on their relative abundance (only ASVs with ≥0.01% relative abundance were kept) before the analyses. The percentage of the microbial community represented by the remaining ASVs after the filtering was 88.35 ± 2.98% (mean ± SD) across samples. High-identity (≥99%) hits were determined by the stringent mapping of ASVs to each reference database. Classification of ASVs was done using the SINTAX classifier. The violin and box plots represent the distribution of percent of ASVs with high-identity hits or genus/species-level classifications for each database across n = 1165 biologically independent samples. Box plots indicate median (middle line), 25th, 75th percentile (box) and the min and max values after removing outliers based on 1.5x interquartile range (whiskers). Outliers have been removed from the box plots to ease visualisation. Different colours are used to distinguish the different databases.Full size imageUsing MiDAS 4 with the SINTAX classifier, it was possible to obtain genus-level classifications for 75.0 ± 6.9% of the GWMC ASVs with ≥0.01% relative abundance (Fig. 3). In comparison, SILVA 138 SSURef NR99, which was the best of the universal reference databases, could only classify 31.4 ± 4.2% of the ASVs to genus-level. When MiDAS 4 was used to classify amplicons from this study, we obtained genus-level classification for 92.0 ± 4.0% of V1–V3 ASVs and 84.8 ± 3.6% of V4 ASVs (Supplementary Fig. 2a). This is close to the theoretical limit set by the phylogenetic signal provided by each amplicon region analyzed20. Improved classifications were also observed for archaeal V4 ASVs (93.3 ± 10.6% for MiDAS 4 vs 69.3 ± 21.3% for SILVA), although no additional archaeal reference sequences were added to the MiDAS database in this study.MiDAS 4 was also able to assign species-level classifications to 40.8 ± 7.1% of the GWMC ASVs. In contrast, the 16S rRNA gene reference database obtained from GTDB SSU r89, which is the only universal reference database that contains a comprehensive species-level taxonomy, only classified 9.9 ± 2.0% of the ASVs (Fig. 3). For the ASVs created in this study, MiDAS 4 provided a species-level classification for 68.4 ± 6.1% of the V1–V3 and 48.5 ± 6.0% of the V4 ASVs (Supplementary Fig. 2a).Based on the large number of WWTPs sampled, their diversity, and the independent evaluation based on the GWMC dataset2, we expect that the MiDAS 4 reference database essentially covers the large majority of bacteria in WWTPs worldwide. Therefore, the MiDAS 4 taxonomy should act as a shared vocabulary for wastewater treatment microbiologists, providing opportunities for cross-study comparisons and ecological studies at high taxonomic resolution.Comparison of the V1–V3 and V4 primer sets for community profiling of WWTPsBefore investigating what factors shape the activated sludge microbiota, we compared short-read amplicon data created for all activated sludge samples belonging to the four main process types (C; C,N; C,N,DN and C,N,DN,P) collected in the Global MiDAS project using two commonly used primer sets that target the V1–V3 or V4 variable region of the 16S rRNA gene. The V1–V3 primers were chosen because the corresponding region of the 16S rRNA gene provides the highest taxonomic resolution of common short-read amplicons20,26, and these primers have previously shown great correspondence with metagenomic data and quantitative fluorescence in situ hybridisation (FISH) results for wastewater treatment systems17. The V4 region has a lower phylogenetic signal, but the primers used for amplification have better theoretical coverage of the bacterial diversity in the SILVA database20,26.The majority of genera (62%) showed less than twofold difference in relative abundances between the two primer sets, and the rest were preferentially detected with either the V1–V3 or the V4 primer (19% for both) (Fig. 4). We observed that several genera of known importance detected in high abundance by V1–V3 were hardly observed by V4, including Acidovorax, Rhodoferax, Ca. Villigracilis, Sphaerotilus and Leptothrix. Similarly, we observed genera abundant with V4 but strongly underestimated by V1–V3, such as Acinetobacter and Prosthecobacter. A complete list of differentially detected genera (Supplementary Data 2) serves as a valuable tool in combination with in silico primer evaluation for deciding which primer pair to use for targeted studies of specific taxa.Fig. 4: Comparison of relative genus abundance based on V1–V3 and V4 region 16S rRNA gene amplicon data.a Mean relative abundance was calculated based on 709 activated sludge samples. Genera present at ≥0.001% relative abundance in V1–V3 and/or V4 datasets are considered. Genera with less than twofold difference in relative abundance between the two primer sets are shown with gray circles, and those that are overrepresented by at least twofold with one of the primer sets are shown in red (V4) and blue (V1–V3). The twofold difference is an arbitrary choice; however, it relates to the uncertainty we usually encounter in amplicon data. Genus names are shown for all taxa present at a minimum of 0.1% mean relative abundance (excluding those with de novo names). b Heatmaps of the most abundant genera with more than twofold relative abundance difference between the two primer sets.Full size imageBecause the V1–V3 primers provide better classification rates at the genus- and species-level (Supplementary Fig. 2a), we primarily focused on this dataset for the following analyses. It should be noted that the V1–V3 primer set performs poorly on anammox bacteria27,28 and does not target archaea at all. To determine the importance of these groups, we estimated their relative read abundance using the V4 amplicon data. Ca. Brocadia and Ca. Anammoximicrobium were the only anammox genera detected, and the latter was never more than 0.6% abundant. Ca. Brocadia was observed in MBBR reactors and granular sludge in anammox reactors with relative read abundances reaching 29%, but it was below 0.1% relative abundance in all but two of the activated sludge samples investigated. For archaea, the relative read abundance was generally low (median = 0.18%), but for a few WWTPs high (up to 11.7%), so archaea should not be neglected in these cases.Process and environmental factors affecting the activated sludge microbiotaAlpha diversity analysis revealed that the rarefied (10,000 read per sample) richness and diversity in activated sludge plants were most strongly affected by process type, industrial load and continent (Supplementary Fig. 3 and Supplementary Note 1). The richness and diversity increased with the complexity of the treatment process, as found in other studies, reflecting the increased number of niches29. In contrast, it decreased with high industrial loads, presumably because industrial wastewater often is less complex and therefore promotes the growth of fewer specialised species7. The effect of continents is presumably caused by the necessary unbalanced sampling of WWTPs and confounded by the effects of plant types and industrial loads.Distance decay relationship (DDR) analyses were used to determine the effect of geographic distance on the microbial community similarity of activated sludge plants with the four main process types (Supplementary Fig. 4 and Supplementary Note 2). We found that distance decay was only effective within shorter geographical distances (2500 km) at the ASV-level, but higher similarities with OTUs clustered at 97% and even more at the genus-level. This suggests that many ASVs are geographically restricted and functionally redundant in the activated sludge microbiota, so different strains or species from the same genus across the world may provide similar functions.To gain a deeper understanding of the factors that shape the activated sludge microbiota, we examined the genus-level taxonomic beta-diversity using principal coordinate analysis (PCoA) and permutational multivariate analysis of variance (PERMANOVA) analyses (Fig. 5 and Supplementary Note 3). We have chosen taxonomic diversity instead of phylogenetic diversity (UniFrac) because many of the important traits are categorical (yes/no) and only conserved at lower taxonomic ranks (genus/species). The analysis was made at the genus-level due to the high classification rate achieved with MiDAS 4 and because genera were less affected by DDR compared to ASVs. We found that the overall microbial community was most strongly affected by continent and temperature in the WWTPs. However, process type, industrial load and the climate zone also had significant impacts. The percentage of total variation explained by each parameter was generally low, indicating that the global WWTPs microbiota represents a continuous distribution rather than distinct states, as observed for the human gut microbiota30.Fig. 5: Effects of process and environmental factors on the activated sludge microbial community structure. Principal coordinate analyses of Bray–Curtis and Soerensen beta-diversity for genera based on V1–V3 amplicon data. Samples are coloured based on metadata.The fraction of variation in the microbial community explained by each variable in isolation was determined by PERMANOVA (Adonis R2-values). Exact P values 0.1% relative abundance in 80% (strict core), 50% (general core) and 20% (loose core) of all activated sludge plants (Fig. 6a).Fig. 6: Identification of core and conditionally rare or abundant taxa based on V1–V3 amplicon data.a Identification of strict, general and loose core genera based on how often a given genus was observed at a relative abundance above 0.1% in WWTPs. b Identification of conditionally rare or abundant (CRAT) genera based on whether a given genus was observed at a relative abundance above 1% in at least one WWTP. The cumulative genus abundance is based on all ASVs classified at the genus-level. All core genera were removed before identification of the CRAT genera. c, d Number of genera and species, respectively, and their abundance in different process types across the global WWTPs. Values for genera and species are divided into strict core, general core, loose core, CRAT, other taxa and unclassified ASVs. The relative abundance of different groups was calculated based on the mean relative abundance of individual genera or species across samples. C carbon removal, C,N carbon removal with nitrification, C,N,DN carbon removal with nitrification and denitrification, C,N,DN,P carbon removal with nitrogen removal and enhanced biological phosphorus removal (EBPR).Full size imageIn addition to the core taxa, we also identified conditionally rare or abundant taxa (CRAT)32 (Fig. 6b). These are taxa typically present in low abundance but occasionally become prevalent, including taxa related to process disturbances, such as bacteria causing activated sludge foaming or those associated with the degradation of specific residues in industrial wastewater. CRAT have only been studied in a single WWTP treating brewery wastewater, despite their potential effect on performance32,33. CRAT are here defined as taxa which are not part of the core, but present in at least one WWTP with a relative abundance above 1%.Core taxa and CRAT were identified for both the V1–V3 and V4 amplicon data to ensure that critical taxa were not missed due to primer bias. We identified 250 core genera (15 strict, 65 general and 170 loose) and 715 CRAT genera (Supplementary Data 4). The strict core genera (Fig. 7) mainly contained genera with versatile metabolisms found in several environments, including Flavobacterium, Novosphingobium and Haliangium. The general core (Fig. 7) included many known bacteria associated with nitrification (Nitrosomonas and Nitrospira), polyphosphate accumulation (Tetrasphaera, Ca. Accumulibacter) and glycogen accumulation (Ca. Competibacter). The loose core contained well-known filamentous bacteria (Ca. Microthrix, Ca. Promineofilum, Ca. Sarcinithrix, Gordonia, Kouleothrix and Thiothrix), but also Nitrotoga, a less common nitrifier in WWTPs.Fig. 7: Percent relative abundance of strict and general core taxa across process types.The taxonomy for the core genera indicates phylum and genus. For general core species, genus names are also provided. De novo taxa in the core are highlighted in red. C carbon removal, C,N carbon removal with nitrification, C,N,DN carbon removal with nitrification and denitrification, C,N,DN,P carbon removal with nitrogen removal and enhanced biological phosphorus removal (EBPR).Full size imageBecause MiDAS 4 allowed for species-level classification, we also identified core and CRAT species based on the same criteria as for genera (Supplementary Fig. 7 and Supplementary Data 4). This revealed 113 core species (0 strict, 9 general and 104 loose). The general core species (Fig. 7) included Nitrospira defluvii and Tetrasphaera midas_s_5, a common nitrifier and PAO, respectively. Arcobacter midas_s_2255, a potential pathogen commonly abundant in the influent wastewater, was also part of the general core34. The loose core contained additional species associated with nitrification (Nitrosomonas midas_s_139 and Nitrospira nitrosa), polyphosphate accumulation (Ca. Accumulibacter phosphatis, Dechloromonas midas_s_173, Tetrasphaera midas_s_45), as well as known filamentous species (Ca. Microthrix parvicella and midas_s_2 (recently named Ca. M. subdominans35), Ca. Villigracilis midas_s_471 and midas_s_9223, Leptothrix midas_s_884). In addition to the core species, we identified 1417 CRAT species. As CRAT taxa are generally found in low abundance and the current study does not include time series or influent data, we cannot say anything conclusive about their general implications for the ecosystem. However, they may be present due to short-term mass immigration25 or specific operational conditions36 and in both cases, potentially affect the plant operation. They should therefore be considered important target for further investigations together with the core taxa.Many core taxa and CRAT can only be identified with MiDAS 4The core taxa and CRAT included a large proportion of MiDAS 4 de novo taxa. At the genus-level, 106/250 (42%) of the core genera and 500/715 (70%) of the CRAT genera had MiDAS placeholder names. At the species-level, the proportion was even higher. Here placeholder names were assigned to 101/113 (89%) of the core species and 1352/1417 (95%) CRAT species. This highlights the importance of a comprehensive taxonomy that includes the uncultured environmental taxa.The core and CRAT taxa cover the majority of the global activated sludge microbiotaAlthough the core taxa and CRAT represent a small fraction of the total diversity observed in the MiDAS 4 reference database, they accounted for the majority of the observed global activated sludge microbiota (Fig. 6c, d). Accumulated read abundance estimates ranged from 57–68% for the core genera and 11–13% for the CRAT, and combined they accounted for 68–79% of total read abundance in the WWTPs depending on process types. The core taxa represented a larger proportion of the activated sludge microbiota for the more advanced process types, which likely reflects the requirement of more versatile bacteria associated with the alternating redox conditions in these types of WWTPs. The remaining fraction, 21–32%, consisted of 6–8% unclassified genera and genera present in very low abundance, presumably with minor importance for the plant performance. The species-level core taxa and CRAT represented 11–24% and 24–33% accumulated read abundance, respectively. Combined, they accounted for almost 50% of the observed microbiota.Global diversity within important functional guildsThe general change from simple to advanced WWTPs with nutrient removal and the transition to water resource recovery facilities (WRRFs) requires increased knowledge about the bacteria responsible for the removal and recovery of nutrients, so we examined the global diversity of well-described nitrifiers, denitrifiers, PAOs and GAOs (Fig. 8). GAOs were included because they may compete with the PAOs for nutrients and thereby interfere with the biological recovery of phosphorus37. Because MiDAS 4 provided species-level resolution for a large proportion of activated sludge microbiota, we also investigated the species-level diversity within genera affiliated with the functional guilds. A complete overview of species in all genera detected in this global study is provided in the MiDAS field guide (https://www.midasfieldguide.org/guide).Fig. 8: Global diversity of genera belonging to major functional groups.The percent relative abundance represents the mean abundance for each country considering only WWTPs with the relevant process types. Countries are grouped based on continent (shifting colour).Full size imageNitrosomonas and potential comammox Nitrospira were the only abundant (≥0.1% average relative abundance) genera found among ammonia-oxidising bacteria (AOBs), whereas both Nitrospira and Nitrotoga were abundant among the nitrite oxidisers (NOBs), with Nitrospira being the most abundant across all countries (Fig. 8). Nitrobacter was not detected, and Nitrosospira was detected in only a few plants in very low abundance (≤0.01% average relative abundance). At the species-level, each genus had 2–5 abundant species (Supplementary Fig. 8). The most abundant and widespread Nitrosomonas species was midas_s_139. However, midas_s_11707 and midas_s_11733 were dominating in a few countries. For Nitrospira, the most abundant species in nearly all countries was N. defluvii. ASVs classified as the comammox N. nitrosa38,39 was also common in many countries across the world. However, because the comammox trait is not phylogenetically conserved at the 16S rRNA gene level38,39, we cannot conclude that these ASVs represent true comammox bacteria. For Nitrotoga, only two species were detected with notable abundance, midas_s_181 and midas_s_9575. Ammonia-oxidising archaea (AOAs) were not detected with MiDAS 4 due to the lack of reference sequences, and because AOAs are not targeted by the V1–V3 primer pair. However, analyses of our V4 amplicon dataset classified with the SILVA database revealed a considerable relative read abundance of AOAs in Malaysia and the Philippines, but absence or low abundance of AOAs in other countries (Supplementary Fig. 9). Other studies have occasionally found AOAs across the world, but generally in lower abundance than AOBs40,41,42. To ensure detection of AOAs with MiDAS 4, we anticipate adding external reference sequences for AOAs in a future release of the database.Denitrifying bacteria are very common in advanced activated sludge plants, but are generally poorly described. Among the known genera, Rhodoferax, Zoogloea and Thauera were most abundant (Fig. 8). Zoogloea and Thauera are well-known floc formers, sometimes causing unwanted slime formation43. Rhodoferax was the most common denitrifier in Europe, whereas Thauera dominated in Asia. Many denitrifiers could not be classified at the species-level (Supplementary Fig. 10), likely due to highly conserved 16S rRNA genes. An exception was Zoogloea, where midas_s_1080 and Z. caeni and were the most abundant species worldwide.EBPR is performed by PAOs, with three genera recognised as important in full-scale WWTPs: Tetrasphaera, Dechloromonas and Ca. Accumulibacter13. According to relative read abundance, all three were found in EBPR plants globally, with Tetrasphaera as the most prevalent (Fig. 8). Dechloromonas was also abundant in nitrifying and denitrifying plants without EBPR, indicating a more diverse ecology. Four recognised GAOs were found globally: Ca. Competibacter, Defluviicoccus, Propionivibrio and Micropruina, with Ca. Competibacter being the most abundant (Fig. 8). Only a few species (2–6 species) in each genus were dominant across the world for both PAOs (Supplementary Fig. 11) and GAOs (Supplementary Fig. 12), except for Ca. Competibacter, which covered ~20 abundant but country-specific species. Among PAOs, the abundant species were Tetrasphaera midas_s_5, Dechloromonas midas_s_173, (recently named D. phosphorivorans) Ca. Accumulibacter midas_s_315, Ca. A. phosphatis and Ca. A. aalborgensis. Interestingly, some of the most abundant PAOs and GAOs were also abundant in the simple process design with C-removal, indicating more versatile metabolisms.Global diversity of filamentous bacteriaFilamentous bacteria are essential for creating strong activated sludge flocs. However, in large numbers, they can also lead to loose flocs and poor settling properties. This is known as bulking, a major operational problem in many WWTPs. Many can also form foam on top of process tanks due to hydrophobic surfaces. Presently, approximately 20 genera are known to contain filamentous species44, and among those, the most abundant are Ca. Microthrix, Leptothrix, Ca. Villigracilis, Trichococcus and Sphaerotilus (Fig. 9). They are all well-known from studies on mitigation of poor settling properties in WWTPs. Interestingly, Leptothrix, Sphaerotilus and Ca. Villigracilis belong to the genera where abundance-estimation depended strongly on primers, with V4 underestimating their abundance (Fig. 3). Ca. Microthrix and Leptothrix were strongly associated with continents, most common in Europe and less in Asia and North America (Fig. 9).Fig. 9: Global diversity of known filamentous organisms.The percent relative abundance represents the mean abundance for each country across all process types. Countries are grouped based on the continent (shifting colour).Full size imageMany of the filamentous bacteria were linked to specific process types (Supplementary Fig. 13), e.g. Ca. Microthrix were not observed in WWTPs with carbon removal only, and Ca. Amarolinea were only abundant in plants with nutrient removal. The number of abundant species within the genera were generally low, with one species in Trichococcus, two in Ca. Microthrix and approximately five in Leptothrix and Ca. Villigracilis (Supplementary Fig. 14). Only five abundant species were observed for Sphaerotilus. However, a substantial fraction of unclassified ASVs was also observed, demonstrating that certain species within this genus are poorly resolved based on the 16S rRNA gene. Ca. Promineofilum was also poorly resolved at the species-level (Supplementary Fig. 15).Conclusion and perspectivesWe present a worldwide collaborative effort to produce MiDAS 4, an ASV-resolved full-length 16S rRNA gene reference database, which covers more than 31,000 species and enables genus- to species-level resolution in microbial community profiling studies. MiDAS 4 covers the vast majority of WWTP bacteria globally and provides a strongly needed common taxonomy for the field, which provides the foundation for comprehensive linking of microbial taxa in the ecosystem with their functional traits. Presently, hundreds of studies are undertaken to combine engineering and microbial aspects of full-scale WWTPs. However, most ASVs or OTUs in these studies are classified at poor taxonomic resolution (family-level or above) due to the use of incomplete universal reference databases. Because many important functional traits are only conserved at high taxonomic resolution (genus- or species-level), this strongly hampers our ability to transfer taxa-specific knowledge from one study to another. This will change with MiDAS 4, and we expect that reprocessing of data from earlier studies may reveal new perspectives into wastewater treatment microbiology. Our online Global MiDAS Field Guide presents the data generated in this study and summarises present knowledge about all taxa. We encourage researchers within the field to contribute new knowledge to MiDAS using the contact link in the MiDAS website (https://www.midasfieldguide.org/guide/contact).The global microbiota of activated sludge plants has been predicted to harbour a massive diversity with up to one billion species2. However, most of these occur at very low abundance and are of little importance for the treatment process. By focusing only on the abundant taxa, we can see that this number is much smaller, i.e., ~1000 genera and 1500 species. We consider these taxa functionally the most important globally, representing a “most wanted list” for future studies. Some taxa are abundant in most WWTPs (core taxa), and others are occasionally abundant in fewer plants (CRAT). The CRAT have received little attention in the field of wastewater treatment, but they can be of profound importance for WWTP performance. Both groups have a high fraction of poorly characterised species. The high taxonomic resolution provided by MiDAS 4 enables us to identify samples where these important core taxa occur in high abundance. This provides an ideal starting point for obtaining high-quality metagenome-assembled genomes (MAGs), isolation of pure cultures, in addition to targeted culture-independent studies to uncover their physiological and ecological roles.Among the known functional guilds, such as nitrifiers or polyphosphate-accumulating organisms, the same genera were found worldwide, with only a few abundant species in each genus. There were differences in the community structure, and the abundance of dominant species was mainly shaped by process type, temperature, and in some cases, continent. This discovery sends an important message to the field: relatively few species are abundant worldwide, so research or operational results can reliably be transferred from one geographical region to another, stimulating the transition from WWTPs to more sustainable WRRFs.The relatively low number of uncharacterised abundant species also shows that it is within our reach to describe them all in terms of identity, physiology, ecology and dynamics, providing the necessary knowledge for informed process optimisation and management. The number of poorly described genera (i.e. those with only a MiDAS placeholder genus name) was 88 among the 250 core genera (35%) and more than 89% at the species-level, so there is still some work to do to link their identities and function. An important step in this direction is the visualisation of the populations. With the comprehensive set of FL-ASVs, it is possible to design highly specific FISH probes, and to critically evaluate the old probes. In the Danish WWTPs, we have successfully done this for groups in the Acidobacteriota42 based on the MiDAS 3 database18. Our recent retrieval of more than 1000 high-quality MAGs from Danish WWTPs with advanced process design is also an important step to link identity to function43. The HQ-MAGs can be linked directly to MiDAS 4 as they contain complete 16S rRNA genes. They cover 62% (156/250) of the core genera and 61% (69/113) of the core species identified in this study. These MAGs may also form the basis for further studies to link identity and function, e.g. by applying metatranscriptomics44 and other in situ techniques such as FISH combined with Raman45,46, guided by the “most wanted” list provided in this study. We expect that MiDAS 4 will have significant implications for future microbial ecology studies in wastewater treatment systems. More

Ellis, E. C., Beusen, A. H. W. & Goldewijk, K. K. Anthropogenic Biomes: 10,000 BCE to 2015 CE. Land. 9(5), 129 (2020).Article 

Google Scholar 
Seto, K. C., Guneralp, B. & Hutyra, L. R. Global forecasts of urban expansion to 2030 and direct impacts on biodiversity and carbon pools. PNAS. 109, 16083–8 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Doherty, T. S., Hays, G. C. & Driscoll, D. A. Human disturbance causes widespread disruption of animal movement. Nat. Ecol. Evol. 5, 513–519 (2021).PubMed 
Article 

Google Scholar 
Frid, A. & Dill, L. Human-caused disturbance stimuli as a form of predation risk. Conserv. Ecol. 6, 11 (2002).
Google Scholar 
Rodgers, J. A. & Schwikert, S. T. Buffer-zone distances to protect foraging and loafing waterbirds from disturbance by personal watercraft and outboard-powered boats. Conserv. Bio. 16, 216–224 (2002).Article 

Google Scholar 
Constantine, R., Brunton, D. H. & Dennis, T. Dolphin-watching tour boats change bottlenose dolphin (Tursiops truncates) behaviour. Biol. Conserv. 117, 299–307 (2004).Article 

Google Scholar 
Gill, J. A., Sutherland, W. J. & Watkinson, A. R. A method to quantify the effects of human disturbance on animal populations. J. Appl. Ecol. 33, 786–792 (1996).Article 

Google Scholar 
King, J. M. & Heinen, J. T. An assessment of the behaviors of overwintering manatees as influenced by interactions with tourists at two sites in central Florida. Biol. Conserv 117, 227–234 (2004).Article 

Google Scholar 
Stockwell, C. A., Bateman, G. C. & Berger, J. Conflicts in national parks: A case study of helicopters and bighorn sheep time budgets at the Grand Canyon. Biol. Conserv 56, 317–328 (1991).Article 

Google Scholar 
Diamond, J. M. The design of a nature reserve system for Indone-Asian New Guinea. In Conservation Biology: The Science of Scarcity and Cliversity (ed. Soule, M.) 485–503 (Sinauer, Sunderland, Massachusetts, 1986).Ceballos, G., García, A. & Ehrlich, P. R. The sixth extinction crisis loss of animal populations and species. J. Cosmol. 8, 1821–1831 (2010).
Google Scholar 
Kerr, J. T. & Deguise, I. Habitat loss and the limits to endangered species recovery. Ecol. Lett. 7, 1163–1169 (2004).Article 

Google Scholar 
Mbora, D. N. M. & McPeek, M. A. Host density and human activities mediate increased parasite prevalence and richness in primates threatened by habitat loss and fragmentation. J. Anim. Ecol. 78, 210–218 (2009).PubMed 
Article 

Google Scholar 
Low, T. The New Nature (Penguin Books Limited, 2003).
Google Scholar 
Baxter-Gilbert, J., Riley, J. L. & Measey, J. Fortune favors the bold toad: Urban-derived behavioral traits may provide advantages for invasive amphibian populations. Behav. Ecol. Sociobiol. 75, 130 (2021).Article 

Google Scholar 
Coleman, J. L. & Barclay, R. M. R. Prey availability and foraging activity of grassland bats in relation to urbanization. J. Mammal. 94, 1111–1122 (2013).Article 

Google Scholar 
Castellano, M. J. & Valone, T. J. Effects of livestock removal and perennial grass recovery on the lizards of a desertified arid grassland. J. Arid Environ. 66, 87e95 (2006).Article 

Google Scholar 
Huey, R. B. Temperature, physiology, and the ecology of reptiles. In Biology of the Reptilia (eds. Gans, C., & Pough, F.H.) Vol. 12. (Academic Press, London, 1982).White, D. et al. Assessing risks to biodiversity from future landscape change. Conserv. Biol. 11, 349360 (1997).Article 

Google Scholar 
Carpio, A. J., Oteros, J., Tortosa, F. S. & Guerrero-Casado, J. Land use and biodiversity patterns of the herpetofauna: The role of olive groves. Acta Oecol. 70, 103–111 (2016).Article 

Google Scholar 
Geyle, H. M., Tingley, R., Amey, A. P. & Chapple, D. G. Reptiles on the brink: Identifying the Australian terrestrial snake and lizard species most at risk of extinction. Pac. Conserv. Biol. 27, 3–12 (2021).Article 

Google Scholar 
Doherty, T. S. et al. Reptile responses to anthropogenic habitat modification: A global meta-analysis. Glob. Ecol. Biogeogr. 29(7), 1265–1279 (2020).Article 

Google Scholar 
Hu, Y., Doherty, T. S. & Jessop, T. S. How influential are squamate reptile traits in explaining population responses to environmental disturbances?. Wildl. Res. 47(3), 249–259 (2020).Article 

Google Scholar 
Poole, G. & Berman, C. An ecological perspective on in-stream temperature: natural heat dynamics and mechanisms of human-caused thermal degradation. Environ. Manag. 27, 787–802 (2001).CAS 
Article 

Google Scholar 
Tang, X. et al. Human activities enhance radiation forcing through surface albedo associated with vegetation in beijing. Remote Sens. 12(5), 837 (2020).Article 

Google Scholar 
Barna, A., Masum, A. K. M., Hossain, M. E., Bahadur, E.H. & Alam, M. S. A study on human activity recognition using gyroscope, accelerometer, temperature and humidity data. In 2019 International Conference on Electrical, Computer and Communication Engineering (ECCE), pp. 1–6 (2019).Moore, M. & Seigel, R. A. No place to nest or bask: Effects of human disturbance on the nesting and basking habits of yellow-blotched map turtles (Graptemys flavimaculata). J. Biol. Conserv. 130(3), 386–393 (2006).Article 

Google Scholar 
Bonnet, X., Naulleau, G. & Shine, R. The dangers of leaving home: Dispersal and mortality in snakes. Biol. Conserv. 89(1), 39–50 (1999).Article 

Google Scholar 
Haxton, T. Road mortality of Snapping Turtles, Chelydra serpentina, in central Ontario during their nesting period. Can. Field-Nat. 114(1), 106–110 (2000).
Google Scholar 
Koenig, J., Shine, R. & Shea, G. L. The ecology of an Australian reptile icon: How do blue-tongued lizards (Tiliqua scincoides) survive in suburbia?. Wildl. Res. 28(3), 214–227 (2001).Article 

Google Scholar 
Uetz, P. How many Reptile species?. Herpetol. Rev. 31, 13–15 (2000).
Google Scholar 
Todd, R. L., Steven, P., Rowland, G. & Oldham, G. Herpetological observations from field expeditions to North Karnataka and Southwest Maharashtra, India. Herpetol. Bull. 112, 17–37 (2010).
Google Scholar 
Sathish Kumar, V. M. The conservation of Indian Reptiles: An approach with molecular aspects. Reptile Rap. 14, 2–8 (2012).
Google Scholar 
Berryman, A. A. & Hawkins, B. A. The refuge as an integrating concept in ecology and evolution. Oikos. 115, 92–196 (2006).Article 

Google Scholar 
Webb, J. K., Pringle, R. M. & Shine, R. How do nocturnal snakes select diurnal retreat sites?. Copeia 2004, 919–925 (2004).Article 

Google Scholar 
Skinner, M. & Miller, N. Aggregation and social interaction in garter snakes (Thamnophis sirtalis sirtalis). Behav. Ecol. Sociobiol. 74, 51 (2020).Article 

Google Scholar 
Aubret, F. & Shine, R. Causes and consequences of aggregation by neonatal tiger snakes (Notechis scutatus, Elapidae). Austral Ecol. 34(2), 210–217 (2009).Article 

Google Scholar 
Myres, B. & Eells, M. Thermal aggregation in Boa constrictor. Herpetologica 24(1), 61–66 (1968).
Google Scholar 
Parrish, J. K. & Edelstein-keshet, L. Coinplexity, pattern, and evolutionary trade-offs in animal aggregation. Science 284, 99–101 (1999).CAS 
PubMed 
Article 

Google Scholar 
Trevesa, A. Theory and method in studies of vigilance and aggregation. Anim. Behav. 60, 711–722 (2000).Article 

Google Scholar 
Greene, H. W. Snakes (University of California Press, 1997).Book 

Google Scholar 
Huey, R. B., Peterson, C. R., Arnold, S. J. & Porter, W. P. Hot rocks and not-so-hot rocks: Retreat-site selection by garter snakes and its thermal consequences. Ecology 70, 931–944 (1989).Article 

Google Scholar 
Christian, K. & Weavers, B. Analysis of activity and energetics of the lizard Varanus rosenbergi. Copeia 1994, 289–295 (1994).Article 

Google Scholar 
Autumn, K. & de Nardo, D. F. Behavioural thermoregulation increases growth rate in nocturnal lizard. J. Herpetol. 29, 157–162 (1995).Article 

Google Scholar 
Milne, T., Bull, C. M. & Hutchinson, M. N. Use of burrows by the endangered pygmy blue-tongue lizard, Tiliqua adelaidensis (Scincidae). Wildl. Res. 30, 523–528 (2003).Article 

Google Scholar 
Sunday, J. M. et al. Thermal-safety margins and the necessity of thermoregulatory behavior across latitude and elevation. PNAS. 111, 5610–5615 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kearney, M., Shine, R. & Porter, W. P. The potential for behavioral thermoregulation to buffer ‘“cold-blooded”’ animals against climate warming. PNAS 106, 3835–3840 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Stevenson, D. J., Dyer, K. J. & Willis-Stevenson, B. A. Survey and monitoring of the eastern indigo snake in georgia. Southeast. Nat. 2(3), 393–408 (2003).Article 

Google Scholar 
Zappalorti, R. T. & Reinert, H. K. Artificial refugia as a habitat-improvement strategy for snake conservation. Contrib. Herpetol. 11, 369–375 (1994).
Google Scholar 
Griffith, B., Scott, J. M., Carpenter, J. W. & Reed, C. Translocation as a species conservation tool: Status and strategy. Science 245, 477–480 (1989).CAS 
PubMed 
Article 

Google Scholar 
Mullin, S. J. Snakes Ecology and Conservation (eds. Stephen, J. M. & Richard, A. S.). (Cornell University Press, 2011).Lei, J., Booth, D. T. & Dwyer, R. G. Spatial ecology of yellow-spotted goannas adjacent to a sea turtle nesting beach. Aust. J. Zool. 65, 77–86 (2017).Article 

Google Scholar 
Ermi, Z. Snakes of China. (Anhui Science and Technology Press, 2006).Schulz, K. D. A Monograph of the Colubrid Snakes of the Genus Elaphe Fitzinger (Czech Republic, Koeltz Scientific Books, 1996).
Google Scholar 
Pallas, P. S. Reise durch verschiedene Provinzen des Russischen Reiches, Vol. 2. 744 (Kaiserl. Akad. Wiss., St. Petersburg, 1773).Auffenberg, W., Arian, Q. N. & Kurshid, N. Preferred habitat, home range and movement patterns of Varanus bengalensis in southern Pakistan. Mertensiella 2, 7–28 (1991).
Google Scholar 
McDiarmid, R. W. Reptile Biodiversity: Standard Methods for Inventory and Monitoring. (University of California Press, 2002).Riley, J. L., Baxter-gilbert, J. H. & Litzgus, J. D. A comparison of three external transmitter attachment methods for snakes. Wildl. Soc. Bull. 41(1), 132–139 (2017).Article 

Google Scholar 
Meine, C., & Archibald, G. The Cranes: Status Survey and Conservation Action Plan (IUCN, 1996).Mori, A. & Toda, M. Body temperature of subtropical snakes at night: How cold is their blood?. Curr. Herpetol. 37(2), 151–157 (2018).Article 

Google Scholar 
Crane, M., Silva, I., Marshall, B. M. & Strine, C. T. Lots of movement, little progress: A review of reptile home range literature. PeerJ 9, e11742 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Calabrese, J. M., Fleming, C. H. & Gurarie, E. ctmm: An R package for analyzing animal relocation data as a continuous-time stochastic process. Methods Ecol. Evol. 7, 1124–1132 (2016).Article 

Google Scholar 
Fleming, C. H. & Calabrese, J. M. A new kernel density estimator for accurate home-range and species-range area estimation. Methods Ecol. Evol. 8, 571–579 (2017).Article 

Google Scholar 
Fleming, C. H. et al. From fine-scale foraging to home ranges: A semivariance approach to identifying movement modes across spatiotemporal scales. Am. Nat. 183, 154–167 (2014).Article 

Google Scholar 
Fleming, C. H., Noonan, M. J., Medici, E. P. & Calabrese, J. M. Overcoming the challenge of small effective sample sizes in home-range estimation. Methods Ecol. Evol. 10, 1679–1689 (2019).Article 

Google Scholar 
Uhlenbeck, G. E. & Ornstein, L. S. On the theory of the Brownian motion. Phys. Rev. 36, 823 (1930).CAS 
MATH 
Article 

Google Scholar 
Bürkner, P. C. brms: An R package for Bayesian multilevel models using Stan. J. Stat. Softw. 80, 1–28 (2017).Article 

Google Scholar 
Bhattacharyya, A. On a measure of divergence between two statistical populations defined by their probability distributions. News Bull. Calcutta Math. Soc. 35, 99–109 (1943).MathSciNet 
MATH 

Google Scholar 
Winner, K. et al. Statistical inference for home range overlap. Methods Ecol. Evol. 9, 1679–1691 (2018).Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna (2016). http://www.R-project.org. Accessed September 2022.Calenge, A. The package ‘“adehabitat”’ for the R software: Tool for the analysis of space and habitat use by animals. Ecol. Modell. 197, 516–519 (2006).Article 

Google Scholar 
Manley, B. F. J., McDonald, L. L. & Thomas, D. L. Resource Selection by Animals: Statistical Design and Analysis for Field Studies (Chapman and Hall, 1993).Book 

Google Scholar  More