Ecology

Thrash JC, Boyd A, Huggett MJ, Grote J, Carini P, Yoder RJ, et al. Phylogenomic evidence for a common ancestor of mitochondria and the SAR11 clade. Sci Rep. 2011;1:9.
Google Scholar 
Ferla MP, Thrash JC, Giovannoni SJ, Patrick WM. New rRNA gene-based phylogenies of the alphaproteobacteria provide perspective on major groups, mitochondrial ancestry and phylogenetic instability. PLoS One. 2013;8:e83383.
Google Scholar 
Giovannoni SJ. SAR11 bacteria: the most abundant plankton in the oceans. Annu Rev Mar Sci. 2017;9:231–55.
Google Scholar 
Zhao X, Schwartz CL, Pierson J, Giovannoni SJ, McIntosh RJ, Nicastro D. Three-dimensional structure of the ultraoligotrophic marine bacterium “Candidatus pelagibacter ubique”. Appl Environ Microbiol. 2017;83:807–16.
Google Scholar 
Giovannoni SJ, DeLong EF, Schmidt TM, Pace NR. Tangential flow filtration and preliminary phylogenetic analysis of marine picoplankton. Appl Environ Microbiol. 1990;56:4.
Google Scholar 
Morris RM, Rappé MS, Connon SA, Vergin KL, Siebold WA, Carlson CA, et al. SAR11 clade dominates ocean surface bacterioplankton communities. Nature. 2002;420:806–10.CAS 

Google Scholar 
Rappé MS, Connon SA, Vergin KL, Giovannoni SJ. Cultivation of the ubiquitous SAR11 marine bacterioplankton clade. Nature. 2002;418:630–3.
Google Scholar 
Grote J, Thrash JC, Huggett MJ, Landry ZC, Carini P, Giovannoni SJ, et al. Streamlining and core genome conservation among highly divergent members of the SAR11 clade. mBio. 2012;3:e00252–12.CAS 

Google Scholar 
Field KG, Gordon D, Wright T, Rappé M, Urback E, Vergin K, et al. Diversity and depth-specific distribution of SAR11 cluster rRNA genes from marine planktonic bacteria. Appl Environ Microbiol. 1997;63:63–70.CAS 

Google Scholar 
Suzuki MT, Beja O, Taylor LT, DeLong EF. Phylogenetic analysis of ribosomal RNA operons from uncultivated coastal marine bacterioplankton. Environ Microbiol. 2001;3:323–31.CAS 

Google Scholar 
Carlson CA, Morris R, Parsons R, Treusch AH, Giovannoni SJ, Vergin K. Seasonal dynamics of SAR11 populations in the euphotic and mesopelagic zones of the northwestern Sargasso Sea. ISME J. 2009;3:283–95.CAS 

Google Scholar 
Brown MV, Lauro FM, DeMaere MZ, Muir L, Wilkins D, Thomas T, et al. Global biogeography of SAR11 marine bacteria. Mol Syst Biol. 2012;8:595.
Google Scholar 
Haro‐Moreno JM, Rodriguez‐Valera F, Rosselli R, Martinez‐Hernandez F, Roda‐Garcia JJ, Gomez ML, et al. Ecogenomics of the SAR11 clade. Environ Microbiol. 2020;22:1748–63.
Google Scholar 
Carini P, White AE, Campbell EO, Giovannoni SJ. Methane production by phosphate-starved SAR11 chemoheterotrophic marine bacteria. Nat Commun. 2014;5:4346.CAS 

Google Scholar 
Sun J, Steindler L, Thrash JC, Halsey KH, Smith DP, Carter AE, et al. One carbon metabolism in SAR11 Pelagic marine bacteria. PLoS One. 2011;6:e23973.CAS 

Google Scholar 
Schwalbach MS, Tripp HJ, Steindler L, Smith DP, Giovannoni SJ. The presence of the glycolysis operon in SAR11 genomes is positively correlated with ocean productivity. Environ Microbiol. 2010;12:490–500.CAS 

Google Scholar 
Sun J, Todd JD, Thrash JC, Qian Y, Qian MC, Temperton B, et al. The abundant marine bacterium Pelagibacter simultaneously catabolizes dimethylsulfoniopropionate to the gases dimethyl sulfide and methanethiol. Nat Microbiol. 2016;1:16065.CAS 

Google Scholar 
Halsey KH, Giovannoni SJ, Graus M, Zhao Y, Landry Z, Thrash JC, et al. Biological cycling of volatile organic carbon by phytoplankton and bacterioplankton: VOC cycling by marine plankton. Limnol Oceanogr. 2017;62:2650–61.CAS 

Google Scholar 
Carlson CA, Giovannoni SJ, Hansell DA, Goldberg SJ, Parsons R, Vergin K. Interactions among dissolved organic carbon, microbial processes, and community structure in the mesopelagic zone of the northwestern Sargasso Sea. Limnol Oceanogr. 2004;49:1073–83.CAS 

Google Scholar 
Wagner S, Schubotz F, Kaiser K, Hallmann C, Waska H, Rossel PE, et al. Soothsaying DOM: a current perspective on the future of oceanic dissolved organic carbon. Front Mar Sci. 2020;7:341.
Google Scholar 
Quinn PK, Bates TS. The case against climate regulation via oceanic phytoplankton sulphur emissions. Nature. 2011;480:51–6.CAS 

Google Scholar 
Bolaños LM, Choi CJ, Worden AZ, Baetge N, Carlson CA, Giovannoni S. Seasonality of the microbial community composition in the North Atlantic. Front Mar Sci. 2021;8:624164.
Google Scholar 
Tucker SJ, Freel KC, Monaghan EA, Sullivan CES, Ramfelt O, Rii YM, et al. Spatial and temporal dynamics of SAR11 marine bacteria across a nearshore to offshore transect in the tropical Pacific Ocean. PeerJ. 2021;9:e12274.
Google Scholar 
Giovannoni SJ, Vergin KL. Seasonality in ocean microbial communities. Science. 2012;335:671–6.CAS 

Google Scholar 
Eren AM, Maignien L, Sul WJ, Murphy LG, Grim SL, Morrison HG, et al. Oligotyping: differentiating between closely related microbial taxa using 16S RRNA gene data. Methods Ecol Evol. 2013;4:1111–9.
Google Scholar 
Vergin K, Done B, Carlson C, Giovannoni S. Spatiotemporal distributions of rare bacterioplankton populations indicate adaptive strategies in the oligotrophic ocean. Aquat Microb Ecol. 2013;71:1–13.
Google Scholar 
Salter I, Galand PE, Fagervold SK, Lebaron P, Obernosterer I, Oliver MJ, et al. Seasonal dynamics of active SAR11 ecotypes in the oligotrophic Northwest Mediterranean Sea. ISME J. 2015;9:347–60.CAS 

Google Scholar 
Ortmann AC, Santos TTL. Spatial and temporal patterns in the Pelagibacteraceae across an estuarine gradient. FEMS Microbiol Ecol. 2016;92:fiw133.
Google Scholar 
Vergin KL, Beszteri B, Monier A, Cameron Thrash J, Temperton B, Treusch AH, et al. High-resolution SAR11 ecotype dynamics at the Bermuda Atlantic Time-series Study site by phylogenetic placement of pyrosequences. ISME J. 2013;7:1322–32.CAS 

Google Scholar 
Needham DM, Fichot EB, Wang E, Berdjeb L, Cram JA, Fichot CG, et al. Dynamics and interactions of highly resolved marine plankton via automated high-frequency sampling. ISME J. 2018;12:2417–32.CAS 

Google Scholar 
Benway HM, Lorenzoni L, White AE, Fiedler B, Levine NM, Nicholson DP, et al. Ocean time series observations of changing marine ecosystems: an era of integration, synthesis, and societal applications. Front Mar Sci. 2019;12:6–393.
Google Scholar 
Steinberg DK, Carlson CA, Bates NR, Johnson RJ, Michaels AF, Knap AH. Overview of the US JGOFS Bermuda Atlantic Time-series Study (BATS): a decade-scale look at ocean biology and biogeochemistry. Deep Sea Res Part II Top Stud Oceanogr. 2001;48:1405–47.CAS 

Google Scholar 
Southward AJ, Langmead O, Hardman-Mountford NJ, Aiken J, Boalch GT, Dando PR, et al. Long-term oceanographic and ecological research in the Western English Channel. In: Advances in marine biology. Elsevier. 2005;47:1–105.Gilbert JA, Field D, Swift P, Newbold L, Oliver A, Smyth T, et al. The seasonal structure of microbial communities in the Western English Channel. Environ Microbiol. 2009;11:3132–9.CAS 

Google Scholar 
Gilbert JA, Steele JA, Caporaso JG, Steinbrück L, Reeder J, Temperton B, et al. Defining seasonal marine microbial community dynamics. ISME J. 2012;6:298–308.CAS 

Google Scholar 
Caporaso JG, Paszkiewicz K, Field D, Knight R, Gilbert JA. The Western English Channel contains a persistent microbial seed bank. ISME J. 2012;6:1089–93.CAS 

Google Scholar 
Warwick-Dugdale J, Solonenko N, Moore K, Chittick L, Gregory AC, Allen MJ, et al. Long-read viral metagenomics captures abundant and microdiverse viral populations and their niche-defining genomic islands. PeerJ. 2019;7:e6800.
Google Scholar 
Vergin KL, Done B, Carlson CA, Giovannoni SJ. Spatiotemporal distributions of rare bacterioplankton populations indicate adaptive strategies in the oligotrophic ocean. Aquat Microb Ecol. 2013;71:1–3.
Google Scholar 
Choi CJ, Jimenez V, Needham DM, Poirier C, Bachy C, Alexander H, et al. Seasonal and geographical transitions in eukaryotic phytoplankton community structure in the Atlantic and Pacific Oceans. Front Microbiol. 2020;11:542372.
Google Scholar 
Bolaños LM, Karp-Boss L, Choi CJ, Worden AZ, Graff JR, Haëntjens N, et al. Small phytoplankton dominate western North Atlantic biomass. ISME J. 2020;14:1663–74.
Google Scholar 
Matsen FA, Kodner RB, Armbrust E. pplacer: linear time maximum-likelihood and Bayesian phylogenetic placement of sequences onto a fixed reference tree. BMC Bioinform. 2010;11:1–6.
Google Scholar 
Treusch AH, Vergin KL, Finlay LA, Donatz MG, Burton RM, Carlson CA, et al. Seasonality and vertical structure of microbial communities in an ocean gyre. ISME J. 2009;3:1148–63.
Google Scholar 
Giovannoni SJ, Rappe MS, Vergin KL, Adair NL. 16S rRNA genes reveal stratified open ocean bacterioplankton populations related to the Green Non-Sulfur bacteria. Proc Natl Acad Sci. 1996;93:7979–84.CAS 

Google Scholar 
Morris RM, Vergin KL, Cho J-C, Rappé MS, Carlson CA, Giovannoni SJ. Temporal and spatial response of bacterioplankton lineages to annual convective overturn at the Bermuda Atlantic Time-series Study site. Limnol Oceanogr. 2005;50:1687–96.CAS 

Google Scholar 
Daims H, Brühl A, Amann R, Schleifer K-H, Wagner M. The domain-specific probe EUB338 is insufficient for the detection of all bacteria: development and evaluation of a more comprehensive probe set. Syst Appl Microbiol. 1999;22:434–44.CAS 

Google Scholar 
Lane DJ. Nucleic acid techniques in bacterial systematics. In: Nucleic acid techniques in bacterial systematics. New York: Wiley; p. 115–75.Callahan BJ, McMurdie PJ, Holmes SP. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 2017;11:2639–43.
Google Scholar 
McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One. 2013;8:e61217.CAS 

Google Scholar 
Eren AM, Borisy GG, Huse SM, Mark Welch JL. Oligotyping analysis of the human oral microbiome. Proc Natl Acad Sci. 2014;111:E2875–84.CAS 

Google Scholar 
Buchholz HH, Michelsen ML, Bolaños LM, Browne E, Allen MJ, Temperton B. Efficient dilution-to-extinction isolation of novel virus–host model systems for fastidious heterotrophic bacteria. ISME J. 2021;15:1585–98.CAS 

Google Scholar 
R Core Team. R: a language and environment for statistical computing. Vienna, Austria; https://www.R-project.org/Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. Package “vegan”.Wickham H. ggplot2: ggplot2. Wiley Interdiscip Rev Comput Stat. 2011;3:180–5.
Google Scholar 
Wang W, Yan J. Shape-restricted regression splines with R package splines2. J Data Sci. 2021;19:498–517.
Google Scholar 
Auladell A, Sánchez P, Sánchez O, Gasol JM, Ferrera I. Long-term seasonal and interannual variability of marine aerobic anoxygenic photoheterotrophic bacteria. ISME J. 2019;13:1975–87.CAS 

Google Scholar 
Ahdesmaki M, Fokianos K, Strimmer K, Ahdesmaki MM. Package ‘GeneCycle’ 2015.Roesch A, Schmidbauer H and Roesch MA. Package ‘WaveletComp.’ 2014.Lomas MW, Bates NR, Johnson RJ, Knap AH, Steinberg DK, Carlson CA. Two decades and counting: 24-years of sustained open ocean biogeochemical measurements in the Sargasso Sea. Deep Sea Res Part II Top Stud Oceanogr. 2013;93:16–32.CAS 

Google Scholar 
Lomas MW, Bates NR, Johnson RJ, Steinberg DK, Tanioka T. Adaptive carbon export response to warming in the Sargasso Sea. Nature Commun. 2022;13:1–0.
Google Scholar 
Sargeant SL, Murrell JC, Nightingale PD, Dixon JL. Basin-scale variability of microbial methanol uptake in the Atlantic Ocean. Biogeosciences. 2018;15:5155–67.CAS 

Google Scholar 
Smyth TJ, Allen I, Atkinson A, Bruun JT, Harmer RA, Pingree RD, et al. Ocean net heat flux influences seasonal to interannual patterns of plankton abundance. PLoS One. 2014;9:e98709.
Google Scholar 
Van de Peer Y. A quantitative map of nucleotide substitution rates in bacterial rRNA. Nucleic Acids Res. 1996;24:3381–91.
Google Scholar 
Baker GC, Smith JJ, Cowan DA. Review and re-analysis of domain-specific 16S primers. J Microbiol Methods. 2003;55:541–55.CAS 

Google Scholar 
Vasileiadis S, Puglisi E, Arena M, Cappa F, Cocconcelli PS, Trevisan M. Soil bacterial diversity screening using single 16S rRNA gene V regions coupled with multi-million read generating sequencing technologies. PLoS ONE. 2012;7:e42671.CAS 

Google Scholar 
Stingl U, Tripp HJ, Giovannoni SJ. Improvements of high-throughput culturing yielded novel SAR11 strains and other abundant marine bacteria from the Oregon coast and the Bermuda Atlantic Time-series study site. ISME J. 2007;1:361–71.CAS 

Google Scholar 
Delmont TO, Kiefl E, Kilinc O, Esen OC, Uysal I, Rappé MS, et al. Single-amino acid variants reveal evolutionary processes that shape the biogeography of a global SAR11 subclade. eLife. 2019;8:e46497.
Google Scholar 
Lévy M, Jahn O, Dutkiewicz S, Follows MJ, d’Ovidio F. The dynamical landscape of marine phytoplankton diversity. J R Soc Interface. 2015;12:20150481.
Google Scholar 
Hellweger FL, van Sebille E, Calfee BC, Chandler JW, Zinser ER, Swan BK, et al. The role of ocean currents in the temperature selection of plankton: insights from an individual-based model. PLoS ONE. 2016;11:e0167010.
Google Scholar 
Giovannoni SJ, Tripp HJ, Givan S, Podar M, Vergin KL, Baptista D, et al. Genome streamlining in a cosmopolitan oceanic bacterium. Science. 2005;309:1242–5.CAS 

Google Scholar 
Brown SN, Giovannoni S, Cho JC. Polyphasic taxonomy of marine bacteria from the SAR11 group Ia: Pelagibacter ubiquis (strain HTCC1062) & Pelagibacter bermudensis (strain HTCC7211). Oregon State University; 2012.Auladell A, Barberán A, Logares R, Garcés E, Gasol JM, Ferrera I. Seasonal niche differentiation among closely related marine bacteria. ISME J. 2022;16:178–89.CAS 

Google Scholar 
Tsementzi D, Wu J, Deutsch S, Nath S, Rodriguez-R LM, Burns AS, et al. SAR11 bacteria linked to ocean anoxia and nitrogen loss. Nature. 2016;536:179–83.CAS 

Google Scholar 
Ruiz-Perez CA, Bertagnolli AD, Tsementzi D, Woyke T, Stewart FJ, Konstantinidis KT. Description of Candidatus Mesopelagibacter carboxydoxydans and Candidatus Anoxipelagibacter denitrificans: nitrate-reducing SAR11 genera that dominate mesopelagic and anoxic marine zones. Syst Appl Microbiol. 2021;44:126185.CAS 

Google Scholar 
Yeh YC, Fuhrman JA. Contrasting diversity patterns of prokaryotes and protists over time and depth at the San-Pedro Ocean Time series. ISME Commun. 2022;13:1–12.
Google Scholar 
McCarthy M, Spillane S, Walsh S, Kendon M. The meteorology of the exceptional winter of 2015/2016 across the UK and Ireland. Weather. 2016;71:305–13.
Google Scholar 
Met Office. UK Climate Projections: Headline Findings. 2021. More

Refinetti, R. The diversity of temporal niches in mammals. Biol. Rhythm Res. 39, 173–192 (2008).
Google Scholar 
Hut, R. A., Kronfeld-Schor, N., van der Vinne, V. & De la Iglesia, H. In search of a temporal niche: Environmental factors. Prog. Brain Res. 199, 281–304 (2012).PubMed 

Google Scholar 
Cox, D., Gardner, A. & Gaston, K. Diel niche variation in mammals associated with expanded trait space. Nat. Commun. 12, 1–10 (2021).
Google Scholar 
Grossnickle, D. M., Smith, S. M. & Wilson, G. P. Untangling the multiple ecological radiations of early mammals. Trends Ecol. Evol. 34, 936–949 (2019).PubMed 

Google Scholar 
Baker, J. & Venditti, C. Rapid change in mammalian eye shape is explained by activity pattern. Curr. Biol. 29, 1082–1088. e1083 (2019).PubMed 

Google Scholar 
Crompton, A., Taylor, C. R. & Jagger, J. A. Evolution of homeothermy in mammals. Nature 272, 333–336 (1978).ADS 
PubMed 

Google Scholar 
Maor, R., Dayan, T., Ferguson-Gow, H. & Jones, K. E. Temporal niche expansion in mammals from a nocturnal ancestor after dinosaur extinction. Nat. Ecol. Evol. 1, 1889–1895 (2017).PubMed 

Google Scholar 
Bennie, J. J., Duffy, J. P., Inger, R. & Gaston, K. J. Biogeography of time partitioning in mammals. Proc. Natl Acad. Sci. USA 111, 13727–13732 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Mccain, C. M. & King, S. R. Body size and activity times mediate mammalian responses to climate change. Glob. Change Biol. 20, 1760–1769 (2014).ADS 

Google Scholar 
Veldhuis, M. P. et al. Predation risk constrains herbivores’ adaptive capacity to warming. Nat. Ecol. Evol. 4, 1069–1074 (2020).PubMed 

Google Scholar 
Riede, S. J., van der Vinne, V. & Hut, R. A. The flexible clock: Predictive and reactive homeostasis, energy balance and the circadian regulation of sleep–wake timing. J. Exp. Biol. 220, 738–749 (2017).PubMed 

Google Scholar 
van der Vinne, V. et al. Maximising survival by shifting the daily timing of activity. Ecol. Lett. 22, 2097–2102 (2019).PubMed 
PubMed Central 

Google Scholar 
Harper, G. & Bunbury, N. Invasive rats on tropical islands: Their population biology and impacts on native species. Glob. Ecol. Conserv. 3, 607–627 (2015).
Google Scholar 
Sovie, A. R., Greene, D. U., Frock, C. F., Potash, A. D. & McCleery, R. A. Ephemeral temporal partitioning may facilitate coexistence in competing species. Anim. Behav. 150, 87–96 (2019).
Google Scholar 
Schoener, T. W. Resource partitioning in ecological communities. Science 185, 27–39 (1974).ADS 
PubMed 

Google Scholar 
Richards, S. A. Temporal partitioning and aggression among foragers: Modeling the effects of stochasticity and individual state. Behav. Ecol. 13, 427–438 (2002).
Google Scholar 
Kronfeld-Schor, N. & Dayan, T. Partitioning of time as an ecological resource. Annu. Rev. Ecol., Evol., Syst. 34, 153–181 (2003).
Google Scholar 
Sunarto, S., Kelly, M., Parakkasi, K. & Hutajulu, M. Cat coexistence in central Sumatra: Ecological characteristics, spatial and temporal overlap, and implications for management. J. Zool. 296, 104–115 (2015).
Google Scholar 
Lima, S. L. & Bednekoff, P. A. Temporal variation in danger drives antipredator behavior: The predation risk allocation hypothesis. Am. Naturalist 153, 649–659 (1999).
Google Scholar 
Beschta, R. L. & Ripple, W. J. Large predators and trophic cascades in terrestrial ecosystems of the western United States. Biol. Conserv. 142, 2401–2414 (2009).
Google Scholar 
Duffy, J. E. Biodiversity and ecosystem function: The consumer connection. Oikos 99, 201–219 (2002).
Google Scholar 
Sinclair, A., Mduma, S. & Brashares, J. S. Patterns of predation in a diverse predator–prey system. Nature 425, 288–290 (2003).ADS 
PubMed 

Google Scholar 
Cunningham, C. X., Scoleri, V., Johnson, C. N., Barmuta, L. A. & Jones, M. E. Temporal partitioning of activity: Rising and falling top‐predator abundance triggers community‐wide shifts in diel activity. Ecography 42, 2157–2168 (2019).
Google Scholar 
Hayward, M. W. & Slotow, R. Temporal partitioning of activity in large African carnivores: Tests of multiple hypotheses. South Afr. J. Wildl. Res. 39, 109–125 (2009).
Google Scholar 
Monterroso, P., Alves, P. C. & Ferreras, P. Catch me if you can: Diel activity patterns of mammalian prey and predators. Ethology 119, 1044–1056 (2013).
Google Scholar 
Schemske, D. W., Mittelbach, G. G., Cornell, H. V., Sobel, J. M. & Roy, K. Is there a latitudinal gradient in the importance of biotic interactions? Annu. Rev. Ecol. Evol. Syst. 40, 245–269 (2009).
Google Scholar 
Rovero, F. et al. A standardized assessment of forest mammal communities reveals consistent functional composition and vulnerability across the tropics. Ecography 43, 75–84 (2020).
Google Scholar 
Ahumada, J. A. et al. Community structure and diversity of tropical forest mammals: Data from a global camera trap network. Philos. Trans. R. Soc. B: Biol. Sci. 366, 2703–2711 (2011).
Google Scholar 
Zhang, J. et al. Trophic interactions among vertebrate guilds and plants shape global patterns in species diversity. Proc. R. Soc. B: Biol. Sci. 285, 20180949 (2018).
Google Scholar 
Beaudrot, L. et al. Local temperature and ecological similarity drive distributional dynamics of tropical mammals worldwide. Glob. Ecol. Biogeogr. 28, 976–991 (2019).
Google Scholar 
Janzen, D. H. Why mountain passes are higher in the tropics. Am. Naturalist 101, 233–249 (1967).
Google Scholar 
Khaliq, I., Hof, C., Prinzinger, R., Böhning-Gaese, K. & Pfenninger, M. Global variation in thermal tolerances and vulnerability of endotherms to climate change. Proc. R. Soc. B: Biol. Sci. 281, 20141097 (2014).
Google Scholar 
Willmer, P., Stone, G. & Johnston, I. Environmental Physiology of Animals (John Wiley & Sons, 2009).Cruz, P., Paviolo, A., Bó, R. F., Thompson, J. J. & Di Bitetti, M. S. Daily activity patterns and habitat use of the lowland tapir (Tapirus terrestris) in the Atlantic Forest. Mamm. Biol. 79, 376–383 (2014).
Google Scholar 
Taylor, W. & Skinner, J. Adaptations of the aardvark for survival in the Karoo: A review. Trans. R. Soc. South Afr. 59, 105–108 (2004).
Google Scholar 
Levy, O., Dayan, T., Porter, W. P. & Kronfeld‐Schor, N. Time and ecological resilience: Can diurnal animals compensate for climate change by shifting to nocturnal activity? Ecol. Monogr. 89, e01334 (2019).
Google Scholar 
Simpson, G. G. Splendid Isolation: The Curious History of South American Mammals Vol. 11 (Yale University Press, 1980).Gutiérrez-González, C. E. & López-González, C. A. Jaguar interactions with pumas and prey at the northern edge of jaguars’ range. PeerJ 5, e2886 (2017).PubMed 
PubMed Central 

Google Scholar 
Porfirio, G., Sarmento, P., Foster, V. & Fonseca, C. Activity patterns of jaguars and pumas and their relationship to those of their potential prey in the Brazilian Pantanal. Mammalia 81, 401–404 (2017).
Google Scholar 
Foster, V. C. et al. Jaguar and puma activity patterns and predator‐prey interactions in four Brazilian biomes. Biotropica 45, 373–379 (2013).
Google Scholar 
Ross, J., Hearn, A., Johnson, P. & Macdonald, D. Activity patterns and temporal avoidance by prey in response to S unda clouded leopard predation risk. J. Zool. 290, 96–106 (2013).
Google Scholar 
Lima, S. L. Nonlethal effects in the ecology of predator-prey interactions. Bioscience 48, 25–34 (1998).
Google Scholar 
Santos, F. et al. Prey availability and temporal partitioning modulate felid coexistence in Neotropical forests. PLoS One 14, e0213671 (2019).PubMed 
PubMed Central 

Google Scholar 
Herrera, H. et al. Time partitioning among jaguar Panthera onca, puma Puma concolor and ocelot Leopardus pardalis (Carnivora: Felidae) in Costa Rica’s dry and rainforests. Rev. de. Biol.ía Tropical 66, 1559–1568 (2018).
Google Scholar 
Pratas‐Santiago, L. P., Gonçalves, A. L. S., da Maia Soares, A. & Spironello, W. R. The moon cycle effect on the activity patterns of ocelots and their prey. J. Zool. 299, 275–283 (2016).
Google Scholar 
Gaynor, K. M., Hojnowski, C. E., Carter, N. H. & Brashares, J. S. The influence of human disturbance on wildlife nocturnality. Science 360, 1232–1235 (2018).ADS 
PubMed 

Google Scholar 
Espinosa, S. & Salvador, J. Hunters landscape accessibility and daily activity of ungulates in Yasuní Biosphere Reserve. Ecuad. Therya 8, 45–52 (2017).
Google Scholar 
Butynski, T. M. Ecological survey of the impenetrable (Bwindi) forest, Uganda, and recommendations for its conservation and management. https://doi.org/10.13140/RG.2.1.1719.0487 (1984).Rovero, F. & Ahumada, J. The Tropical Ecology, Assessment and Monitoring (TEAM) Network: An early warning system for tropical rain forests. Sci. Total Environ. 574, 914–923 (2017).ADS 
PubMed 

Google Scholar 
Barnosky, A. D., Koch, P. L., Feranec, R. S., Wing, S. L. & Shabel, A. B. Assessing the causes of late Pleistocene extinctions on the continents. Science 306, 70–75 (2004).ADS 
PubMed 

Google Scholar 
Gorczynski, D. et al. Tropical mammal functional diversity increases with productivity but decreases with anthropogenic disturbance. Proc. R. Soc. B 288, 20202098 (2021).PubMed 
PubMed Central 

Google Scholar 
Frey, S., Fisher, J. T., Burton, A. C. & Volpe, J. P. Investigating animal activity patterns and temporal niche partitioning using camera‐trap data: Challenges and opportunities. Remote Sens. Ecol. Conserv. 3, 123–132 (2017).
Google Scholar 
Bivand, R. et al. Maptools: Tools for Handling Spatial Objects. R package version 1.1-4. http://maptools.r-forge.r-project.org/reference/index.html (2021).Ensing, E. P. et al. GPS based daily activity patterns in European red deer and North American elk (Cervus elaphus): Indication for a weak circadian clock in ungulates. PLoS One 9, e106997 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Vazquez, C., Rowcliffe, J. M., Spoelstra, K. & Jansen, P. A. Comparing diel activity patterns of wildlife across latitudes and seasons: Time transformations using day length. Methods Ecol. Evol. 10, 2057–2066 (2019).
Google Scholar 
Rowcliffe, J. M., Kays, R., Kranstauber, B., Carbone, C. & Jansen, P. A. Quantifying levels of animal activity using camera trap data. Methods Ecol. Evol. 5, 1170–1179 (2014).
Google Scholar 
Rowcliffe, J. M. Activity: Animal Activity Statistics. R package version 1.3.2. https://cran.r-project.org/package=activity (2022).Faurby, S. et al. PHYLACINE 1.2: The phylogenetic atlas of mammal macroecology. Ecology 99, 2626 (2018).PubMed 

Google Scholar 
Wilman, H. et al. EltonTraits 1.0: Species-level foraging attributes of the world’s birds and mammals. Ecology 95, 2027–2027 (2014).
Google Scholar 
Elff, M., Heisig, J. P., Schaeffer, M. & Shikano, S. Multilevel analysis with few clusters: Improving likelihood-based methods to provide unbiased estimates and accurate inference. Br. J. Polit. Sci. 51, 412–426 (2020).Elff, M. Mclogit: mixed conditional logit models. R package version 0.5. 1. https://github.com/melff/mclogit/ (2018).Burnham, K & Anderson, D. Model Selection and Multi-model Inference 2nd edn, Vol. 63, 10 (Springer-Verlag 2004).Carbone, C., Teacher, A. & Rowcliffe, J. M. The costs of carnivory. PLoS Biol. 5, e22 (2007).PubMed 
PubMed Central 

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

Google Scholar 
Meredith, M. & Ridout, M. Overlap: Estimates of coefficient of overlapping for animal activity patterns. R package version 0.2. 4, https://cran.r-project.org/package=overlap (2014).Ridout, M. S. & Linkie, M. Estimating overlap of daily activity patterns from camera trap data. J. Agric., Biol., Environ. Stat. 14, 322–337 (2009).MathSciNet 
MATH 

Google Scholar 
RStudio Team. RStudio: Integrated Development for R (PBC, Boston, MA, 2020). More

Among the 325 municipalities in Greece during the period 2010–2021, WNV events, defined as the occurrence of at least one laboratory-confirmed human WNV case during a specific year, were reported in 154 (47%) municipalities, while the remaining 171 did not report any WNV case. WNV events were reported for a period ranging from one to eight years: 54 (35%) municipalities reported laboratory-confirmed WNV cases in only one year, 38 (25%) in two years, 30 (19%) in three years, 12 (8%) in four years, 10 (6%) in five years, 6 (4%) in six years, 1 (1%) in seven years, and 3 (2%) in eight years. This means that in 60% of the positive areas (82 municipalities out of 154), WNV appeared at most for two years, in 27% (42 out of 154) between three and four years, and in the remaining 13% (20 out of 154) for five years or more. Considering the total number of reported laboratory-confirmed human WNV cases across the twelve years (Fig. 1), in approximately 50% of the positive municipalities (78 out of 154), at most 4 cases were reported: 1, 2, 3, and 4 WNV cases were reported in 24, 32, 11, and 11 municipalities, respectively. Overall, 39 municipalities recorded a number of WNV cases ranging from 5 to 10 (third quartile), 34 a number ranging from 11 to 46, while the remaining 3 municipalities recorded a number of WNV cases equal to 56, 71 and 94.Figure 1(a) Map of Greece with total numbers of laboratory-confirmed human WNV cases throughout the 12-year period 2010–2021, with breakdowns by municipality. White denotes municipalities where no human cases were observed. (b) Map of Greece with total numbers of modelled human WNV cases throughout the 12-year period 2010–2021, with breakdowns by municipality. White denotes municipalities where no human cases were modelled.Full size imageModel evaluation and comparison with MIMESISWe investigated the ability of the MIMESIS-2 model to correctly identify the occurrence of WNV events, both in space and time, and its capacity to quantify the annual number of human WNV cases and the timing of the first WNV event in the year. The performance of many quantities of interest, such as the severity and timing of occurrence of human WNV cases, was also compared output from the original MIMESIS model26.Occurrence of WNV eventsStarting with the spatial analysis, we considered the fit of the model to replicate the observed 385 WNV events out of 3,900 (325*12) possible events across municipalities. MIMESIS-2 was able to correctly identify 356 of them, generated only one false alarm, and correctly modelled 3,514 true negatives.The performance of MIMESIS-2 was then evaluated according to four indices: the probability of detection (POD), false alarm rate (FAR), miss rate (MIS), and critical success index (CSI), described in the Methods section. For the POD, MIS and CSI, we considered the 154 municipalities with at least one reported and laboratory-confirmed human WNV case over the 12-year period, while for the FAR, we considered the 153 municipalities where at least one human WNV case was modelled over the same period. We split the (0.0–1.0) index interval into five equally sized bins to derive for each index, the fraction of municipalities falling into each bin. Both the POD and CSI were above 0.8 for 139 municipalities out of 154, while the MIS was below 0.2 for 142 municipalities (out of 154) and the FAR was always below 0.2, with one false alarm produced in a municipality where WNV events were observed in eight out of twelve years (Table 1).Table 1 Capacity of the MIMESIS-2 model to correctly model laboratory-confirmed human WNV cases.Full size tableWe also analysed how the model performed in different years by studying the multiannual evolution of the indices. Both the aggregated POD and CSI were equal to 0.92, with annual variations ranging from 0.72 (2021) to 1 (2011 and 2014). The aggregated MIS was 0.08, ranging from 0.0 (2011 and 2014) to 0.28 (2021). The FAR was virtually 0, being always equal to 0.0, with the only exception being 2017, when it was 0.1 (Table 2).Table 2 Capacity of the MIMESIS-2 model to correctly model laboratory-confirmed human WNV cases by year and in the whole observed time period.Full size tableMagnitude and timing of WNV events: performance and comparison with MIMESISTo evaluate the ability of MIMESIS-2 to capture the magnitude and timing of WNV events, we first considered the discrepancy between the overall number of observed and modelled WNV cases during the 12-year period for each municipality. Out of the 153 municipalities where at least one case was modelled across the 12 years, 76 (50%) had at most 4 modelled cases of WNV: 1, 2, 3, and 4 WNV cases were modelled in 22, 31, 13, and 10 municipalities, respectively. In 42 municipalities, the number of modelled cases ranged from 5 to 10 (which, as for the observed WNV cases, coincided with the third quartile), and in 32 municipalities, the number ranged from 11 to 47, while the remaining 3 municipalities had 55, 70, and 99 modelled cases (Fig. 1).The MIMESIS-2 model closely replicated the total number of laboratory-confirmed WNV cases during the 12-year period. When considering only the 154 municipalities that recorded at least one WNV event during the considered period (excluding the true negatives), for 140 of them, the modelled number of cases fell within a ± 10% error range of the observed value, whereas for 149 the modelled number of cases fell within the ± 25% error margin. Only two municipalities showed a percent error above 50%. These were particular instances where only one WNV case was reported throughout the considered period, while MIMESIS-2 fitted zero human cases. For the original MIMESIS model, 63 and 84 municipalities fell within the ± 10% and ± 25% error margins, respectively, while 31 municipalities—mainly those where few cases were observed— had a relative error ≥ 100% (Fig. 2).Figure 2(a) For MIMESIS-2, modelled (IHMOD) vs. observed (IHOBS) human WNV cases in each municipality in the period 2010–2021. The inner black line represents the main diagonal where ideally the points would lie in case of perfect fit, while the dashed green, black and red lines represent, respectively, the ± 10%, ± 25% and ± 50% error margin. (b) Same quantities for MIMESIS. (c). Breakdown of the week of first WNV incidence by year. Plotted are the modelled quantities (WYMOD) for each of the 325 municipalities on the y-axis and the observed quantities (WYOBS) on the x-axis. The continuous line represents the main diagonal where ideally the points would lie in case of perfect fit, while the dashed lines represent the ± 4-week error margins.Full size imageTo further evaluate the bias of the model across all municipalities and years, we explored the difference between the yearly modelled and observed human WNV cases both with MIMESIS-2 and the original MIMESIS (IHMOD-IHOBS) across municipalities. In MIMESIS-2, we excluded 3,514 true negative cases to avoid distorted conclusions. For the remaining 386 cases, the mean bias was -0.04 indicating a possibly unbiased model, with the standard deviation (SD) of the residuals equal to 0.66 (original MIMESIS: mean bias 0.33, SD 2.07, after removing 3,387 true negatives) (Supplementary Fig. 1).Across the 325 municipalities and the 12 years, 385 WNV events were observed, while on 3,515 occasions, no laboratory-confirmed human WNV cases were reported; on 162 occurrences, 1 case was reported, and on 67 and 39 occasions, 2 and 3 cases were reported, respectively. The maximum yearly number of human WNV cases observed in a single municipality was 38. Considering the modelled human WNV cases with MIMESIS-2, the distribution of the 356 hits ranged between 1 and 37 modelled cases, closely mimicking the distribution of the observed cases, since 1, 2 and 3 human WNV cases were modelled on 129, 72 and 37 occasions, respectively. For the 29 misses, the observed numbers of human cases were 1 (24 times), 2 (3 times), or 3 (2 times). The only false alarm was produced in the Pellas municipality, where WNV events were observed in 8 out of the 12 years.We evaluated the timing of the first occurrence of WNV in humans for any municipality and year. Ignoring the municipalities with zero cases, the observed and MIMESIS-2-modelled first WNV cases occurred between weeks 22 and 44 and weeks 24 and 36, respectively. Modelled values tended to be dispersed around the observed ones: excluding the 3514 true negatives, 290 (75.13%) of the remaining 386 cases fell into the ± 4-week error margins from the observed cases (Fig. 2). This translated into a much lower bias of the week of first appearance (WYMOD-WYOBS) with respect to MIMESIS (Supplementary Fig. 2).Case study: The Pellas municipalityIn addition to presenting the overall performance of the model throughout different years and Greek municipalities, we highlight here the capacity of the model to capture population-specific behaviour and epidemiological features, such as the force of infection, that is, the rate at which susceptible humans, birds, and mosquitoes become infected, by presenting a single municipality case study for the municipality of Pellas. The Pellas municipality had the highest number of observed WNV cases over the 12-year period with a total of 94 human WNV cases, 38 in 2010, 16 in 2018, and 13 in 2021, no cases from 2014 to 2017, and between 4 to 8 cases in the remaining years.We considered the impact arising from the changes in parameters defining the forces of infection. In addition to the introduction of bird (({psi }_{B})) and human (({psi }_{H})) host selections, changes included modifications for the mosquito-to-bird (({p}_{M})) and bird-to-mosquito (({p}_{B})) probabilities of transmission, whose values were made temperature-dependent following Vogels et al.21, and the replacement of the mosquito-to-bird (({varphi }_{B})) and mosquito-to-human (({varphi }_{H})) ratios with their dynamic counterparts, ({N}_{M}/{N}_{B}) and ({N}_{M}/{N}_{H}), respectively (Fig. 3). We used the May–October period for the 12 years that were considered, because this is the part of the year when Culex pipiens mosquitoes are reproductively active and the majority of human WNV cases are reported. In each year of the 12-year period, ({p}_{M}) started from 0.02, reached its peak—ranging from 0.16 to 0.25—in midsummer, and then decreased to the initial values (in the original model, ({p}_{M}=0.9)). Similarly, ({p}_{B}) started from 0.28, peaked in the same time interval—with maximal values ranging from 0.51 to 0.56—and then returned to the initial values (in the original model, ({p}_{B}=0.125)). Additionally, the dynamic specifications of ({varphi }_{B}) and ({varphi }_{H}) were shown to play an important role. Whereas in MIMESIS ({varphi }_{B}=30), in MIMESIS-2 the values started at approximately 8.6 and peaked in late summer when more human WNV cases are reported, reaching values of approximately 57, before decreasing to values ranging from 31.15 to 41.60 in late October. In MIMESIS, ({varphi }_{H}) was calibrated at the municipality level, and for Pellas municipality, it was 0.0001, whereas the dynamic counterpart in MIMESIS-2 showed a temporal evolution with a shape (but different scale) similar to that of ({varphi }_{B}), starting from values of approximately 1, peaking in late summer to values of approximately 7, and then decreasing to values of approximately 4 in late October.Figure 3The temporal evolution during May to October of the (a) mosquito-to-bird probability of transmission, ({p}_{M}), (b) bird-to-mosquito probability of transmission, ({p}_{B},) (c) mosquito-to-bird ratio, ({varphi }_{B},) and (d) mosquito-to-human ratio, ({varphi }_{H},) for both MIMESIS-2 across different years and MIMESIS for each of the 12 years from 2010 to 2021. The plots refer to the simulations for the municipality of Pella.Full size imageChanges in these parameters enter into the expression for the forces of infection. It is of major practical interest to investigate how the values for the forces of infection resulting from MIMESIS-2 may vary for different values of the relative abundance of the vectors with respect to the corresponding carrying capacity and the temperature in different months (Fig. 4). As expected, all forces of infection increased with both the temperature and the relative abundance of the infectious vertebrate hosts. It is worth noting the importance of day length, as this affects the fraction of nondiapausing mosquitoes, ({delta }_{M}), and causes the forces of infection, all other things being equal, to be potentially higher in June and July than in the other months. However, in these two months, the modelled forces of infection tend to be smaller than those in August due to the lower abundance of infectious hosts.Figure 4Contour plots of the forces of infection for May to September for different values of the relative abundance of infected hosts/vectors with respect to the carrying capacity and the temperature. All the other quantities were fixed to the amounts obtained in the simulations for Pellas municipality for 2021 at the end of the corresponding month. (a) Bird-to-mosquito force of infection (({lambda }_{BM})) as a function of the relative abundance of infected birds (({I}_{B})) with respect to the bird carrying capacity (({K}_{B})) and temperature. (b) Mosquito-to-bird force of infection (({lambda }_{MB})) as a function of the relative abundance of infected mosquitoes (({I}_{M})) with respect to the mosquito carrying capacity (({K}_{M})) and temperature. (c) Mosquito-to-human force of infection (({lambda }_{MH})) as a function of the relative abundance of infected mosquitoes (({I}_{M})) with respect to the mosquito carrying capacity (({K}_{M})) and temperature. The ranges for ({I}_{B}/{K}_{B}) and ({I}_{M}/{K}_{M}) were fixed, increasing the maximum modelled value by 20% for the considered period, while the range for the temperature was chosen considering that in the period of interest, the average daily temperature ranged from 16.6 to 27.1 degrees Celsius. Black crosses represent the modelled values for 2021.Full size imageThe bird-to-mosquito force of infection, ({uplambda }_{BM}), took values on the order of 10–4, with possible peaks of approximately 7 × 10–4 in the case of high temperature and high prevalence of birds in June and July, which were nevertheless not reached due to a low abundance of infected birds in that period. Considering the months of July and August 2021 for illustrative purposes, the resulting modelled values were 1.20 × 10–4 and 2.19 × 10–4, respectively, with the increase in August explained by a higher abundance of infected birds in that period. It is worth noting that if the infection across birds had a lead period of two weeks, the resulting ({uplambda }_{BM}) in July would become 3.91 × 10–4 (+ 226%), while an increase in the average temperature in August by 1 °C would result in ({uplambda }_{BM})= 2.36 × 10–4 (+ 8%). The mosquito-to-bird, ({uplambda }_{MB}), and mosquito-to-human, ({uplambda }_{MH}), forces of infection showed similar qualitative behaviours, albeit at different scales, and in this case, they were higher in August due to a higher prevalence of infected Culex mosquitoes in that month. More specifically, ({uplambda }_{MB}) equalled 1.06 × 10–3 and 1.22 × 10–3 at the end of July and August, respectively, and an expected two weeks for the infection of mosquitoes would result in ({uplambda }_{MB})=4.15 × 10–3 (+ 292%) at the end of July, while an increase in the average temperature in August by 1 °C would result in ({uplambda }_{MB})= 1.31 × 10–3 (+ 7%) at the end of August. Finally, ({uplambda }_{MH})=2.86 × 10–6 at the end of July, while ({uplambda }_{MH})=3.26 × 10–6 at the end of August, with the anticipation of the infection among mosquitoes by two weeks resulting in ({uplambda }_{MH})=1.12 × 10–5 (+ 290%) and an increase in the average August temperature by 1 °C leading to ({uplambda }_{MH})=3.51 × 10–6 (+ 8%). It is worth recalling that since we calibrated the model on the number of reported laboratory-confirmed human WNV cases, ({uplambda }_{MH}) represents the rate at which susceptible humans contract infection and become symptomatic leading to a recorded human WNV case.We explored changes in the populations of infectious hosts and the total population number for both mosquitoes and birds over 2010–2021 for the period spanning from May to October (Figs. 5 and 6). The population of infected mosquitoes (({I}_{M})) was initialised by calibration (see the Methods section). Each year, after a short period in which the population of infected mosquitoes slightly decreased due to a very small number of infectious birds (({I}_{B})) that prevented the infection from spreading, it started growing substantially during summer, reaching its peak in late summer, coinciding with the period when most human cases were recorded. The observed increase in ({I}_{M}) was combined with the growth ({I}_{B}) at approximately the same time (with a slightly anticipated peak), which had an amplification effect on the spread of the infection. Both ({I}_{M}) and ({I}_{B}) showed significant yearly variation, with higher modelled numbers in years where more human WNV cases were reported. The modelled total population of mosquitoes (({N}_{M})) did not show significant interannual variability, always peaking in late summer. Finally, the overall population of birds (({N}_{B})) did not show any variability in the first part of the year, when an increase due to immigration and offspring generation was observed, whereas it had a moderate interannual variability in the second half of the year. These differences may be due to heterogeneous numbers of observed infected, dead and immune birds.Figure 5The temporal evolution during May to October of (a) the number of infected mosquitoes modelled by MIMESIS-2 (({I}_{M})), (b) the total number of mosquitoes modelled by MIMESIS-2 (({N}_{M}))(,) (c) the number of infected birds modelled by MIMESIS-2 (({I}_{B}))(,) and (d) the total number of birds modelled by MIMESIS-2 (({N}_{B})) for each of the 12 years from 2010 to 2021. The plots refer to the simulations for the municipality of Pella.Full size imageFigure 6The temporal evolution during May to October of (a) the ratio between the number of WNV-infected mosquitoes modelled by MIMESIS-2 (({I}_{M,MIM-2})) and the ratio modelled by MIMESIS (({I}_{M,MIM})), (b) the ratio between the total number of mosquitoes modelled by MIMESIS-2 (({N}_{M,MIM-2})) and the ratio modelled by MIMESIS (({N}_{M,MIM}))(,) (c) the ratio between the number of infected birds modelled by MIMESIS-2 (({I}_{B,MIM-2})) and the ratio modelled by MIMESIS (({I}_{B,MIM}))(,) and (d) the ratio between the total number of birds modelled by MIMESIS-2 (({N}_{B,MIM-2})) and the ratio modelled by MIMESIS (({N}_{B,MIM})) for each of the years from 2010 to 2021. The plots refer to the simulations for the municipality of Pella.Full size imageComparison of these population dynamics with those of MIMESIS revealed interesting patterns (Fig. 6). Considering the relative number of mosquitoes in MIMESIS-2 with respect to MIMESIS, the populations in MIMESIS tended to grow faster due to a higher mosquito carrying capacity (({K}_{M})) in the original model (({K}_{M}) ≈ 8.3 × 105 in MIMESIS versus ({K}_{M}) ≈ 2.4 × 105 in MIMESIS-2), resulting in a decrease in the ratio between the amounts modelled by MIMESIS-2 and the ones modelled by MIMESIS. Significant interannual variability could be seen in the first part of the year for infectious mosquitoes, where different initial calibration values played an important role. For the populations of birds, until midsummer, the overall number modelled by MIMESIS-2 tended to be approximately 1/4 that of MIMESIS, while as of July, different patterns were observed due to the higher mortality of birds in the original MIMESIS model. In years with higher virus spread, higher mortality was reflected in a sharper decrease in bird populations; therefore, the ratio between the population modelled by MIMESIS-2 and that modelled by MIMESIS increased up to approximately 0.6 (2010). More

Colpaert, J. V. In Stress in Yeasts and Filamentous Fungi Vol. 27 (eds Avery, S. V. et al.) 157–173 (Elsevier, Academic Press, 2008).Chapter 

Google Scholar 
Dalvi, A. A. & Bhalerao, S. A. Response of plants towards heavy metal toxicity: An overview of avoidance, tolerance and uptake mechanism. Ann. Plant Sci. 2, 362–368 (2013).
Google Scholar 
Langer, I. et al. Ectomycorrhizal impact on Zn accumulation of Populus tremula L. grown in metalliferous soil with increasing levels of Zn concentration. Plant Soil 355, 283–297 (2012).Article 

Google Scholar 
Agrahar-Murugkar, D. & Subbulakshmi, G. Nutritional value of edible wild mushrooms collected from the Khasi hills of Meghalaya. Food Chem. 89, 599–603 (2005).Article 

Google Scholar 
Das, N., Charumathi, D. & Vimala, R. Effect of pretreatment on Cd2+ biosorption by mycelial biomass of Pleurotus florida. Afr. J. Biotechnol. 6, 2555–2558 (2007).Article 

Google Scholar 
Dulay, R. M. R. et al. Effects and myco-accumulation of lead (Pb) in five Pleurotus mushrooms. Int. J. Biol. Pharm. Allies Sci. 4, 1664–1677 (2015).
Google Scholar 
Prasad, A. A., Varatharaju, G., Anushri, C. & Dhivyasree, S. Biosorption of lead by Pleurotus florida and Trichoderma viride. Br. Biotechnol. J. 3, 66–78 (2013).Article 

Google Scholar 
Tay, C. C. et al. Biosorption of cadmium ions using Pleurotus ostreatus: Growth kinetics, isotherm study and biosorption mechanism. Korean J. Chem. Eng. 28, 825–830 (2011).Article 

Google Scholar 
Kulshreshtha, S., Mathur, N. & Bhatnagar, P. Mushroom as a product and their role in mycoremediation. AMB Express 4, 1–7 (2014).Article 

Google Scholar 
Khan, I. et al. Mycoremediation of heavy metal (Cd and Cr)-polluted soil through indigenous metallotolerant fungal isolates. Environ. Monit. Assess. 191, 1–11 (2019).
Google Scholar 
Mosa, K. A., Saadoun, I., Kumar, K., Helmy, M. & Dhankher, O. P. Potential biotechnological strategies for the cleanup of heavy metals and metalloids. Front. Plant Sci. 7, 303 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Velásquez, L. & Dussan, J. Biosorption and bioaccumulation of heavy metals on dead and living biomass of Bacillus sphaericus. J. Hazard. Mater. 167, 713–716 (2009).Article 
PubMed 

Google Scholar 
Fawzy, E. M., Abdel-Motaal, F. F. & El-zayat, S. A. Biosorption of heavy metals onto different eco-friendly substrates. J. Toxicol. Environ. Health Sci. 9, 35–44 (2017).
Google Scholar 
Kumar, V. & Dwivedi, S. K. Mycoremediation of heavy metals: Processes, mechanisms, and affecting factors. Environ. Sci. Pollut. Res. 28, 10375–10412 (2021).Article 

Google Scholar 
Barros, L., Baptista, P., Estevinho, L. M. & Ferreira, I. C. Bioactive properties of the medicinal mushroom Leucopaxillus giganteus mycelium obtained in the presence of different nitrogen sources. Food Chem. 105, 179–186 (2007).Article 

Google Scholar 
Kim, H. G. et al. Phellinus linteus inhibits inflammatory mediators by suppressing redox-based NF-κB and MAPKs activation in lipopolysaccharide-induced RAW 264.7 macrophage. J. Ethnopharmacol. 114, 307–315 (2007).Article 
PubMed 

Google Scholar 
Sarikurkcu, C., Tepe, B. & Yamac, M. Evaluation of the antioxidant activity of four edible mushrooms from the Central Anatolia, Eskisehir–Turkey: Lactarius deterrimus, Suillus collitinus, Boletus edulis, Xerocomus chrysenteron. Bioresour. Technol. 99, 6651–6655 (2008).Article 
PubMed 

Google Scholar 
Giannopolitis, C.N. & Ries, S.K. Superoxide Dismutases II. Purification and quantitative relationship with water-soluble protein in seedlings.Plant Physiol. 59, 315–318 (1977).Article 
PubMed 
PubMed Central 

Google Scholar 
Dhindsa, R.S., Plumb-Dhindsa, P. & Thorpe, T.A. Leaf senescent: correlated with increased levels of membrane permeability and lipid peroxidation, and decreased levels of superoxide dismutase and catalase. J. Exp. Bot. 32, 93–101 (1981).Article 

Google Scholar 
Marx, D. H. Ectomycorrhizae as biological deterrents to pathogenic root infections. Annu. Rev. Phytopathol. 10, 429–454 (1972).Article 
PubMed 

Google Scholar 
Calvillo-Medina, R. P. Determination of fungal tolerance index to heavy metals and heavy metal resistance tests. Bio-Protoc. 11, e4218 (2021).PubMed 
PubMed Central 

Google Scholar 
Ouzounidou, G., Eleftheriou, E. & Karataglis, S. Ecophysical and ultrastructural effects of copper in Thlaspi ochroleucum (Cruciferae). Can. J. Bot. 70, 947–957 (1992).Article 

Google Scholar 
Carrillo-González, R., González-Chávez, M. & del Carmen, A. Tolerance to and accumulation of cadmium by the mycelium of the fungi Scleroderma citrinum and Pisolithus tinctorius. Biol. Trace Elem. Res. 146, 388–395 (2012).Article 
PubMed 

Google Scholar 
Mkhize, S. S., Simelane, M. B. C., Gasa, N. L. & Pooe, O. J. Evaluating the antioxidant and heavy metal content of Pleurotus ostreatus mushrooms cultivated using sugar cane agro-waste. Pharmacogn. J. 13, 1–9 (2021).Article 

Google Scholar 
Kumar, V. In Microbial Cell Factories (eds Sharma, D. & Saharan, B. S.) 149–174 (CRC Press, 2018).Chapter 

Google Scholar 
El-Sayed, M. T., Ezzat, S. M., Taha, A. S. & Ismaiel, A. A. Iron stress response and bioaccumulation potential of three fungal strains isolated from sewage-irrigated soil. J. Appl. Microbiol. 132, 1936–1953 (2022).Article 
PubMed 

Google Scholar 
Li, X. et al. Mechanisms of Cd and Cr removal and tolerance by macrofungus Pleurotus ostreatus HAU-2. J. Hazard. Mater. 330, 1–8 (2017).Article 
ADS 
PubMed 

Google Scholar 
Chuang, H.-W., Wang, I.-W., Lin, S.-Y. & Chang, Y.-L. Transcriptome analysis of cadmium response in Ganoderma lucidum. FEMS Microbiol. Lett. 293, 205–213 (2009).Article 
PubMed 

Google Scholar 
Jones, D. & Muehlchen, A. Effects of the potentially toxic metals, aluminium, zinc and copper on ectomycorrhizal fungi. J. Environ. Sci. Health A Environ. Sci. Eng. Technol. 29, 949–966 (1994).
Google Scholar 
Tam, P. C. Heavy metal tolerance by ectomycorrhizal fungi and metal amelioration by Pisolithus tinctorius. Mycorrhiza 5, 181–187 (1995).Article 

Google Scholar 
Purkayastha, R., Mitra, A. K. & Bhattacharyya, B. Uptake and toxicological effects of some heavy metals on Pleurotus sajor-caju (Fr.) Singer. Ecotoxicol. Environ. Saf. 27, 7–13 (1994).Article 
PubMed 

Google Scholar 
Akkin, N. Mycoremediation by oyster mushroom. Acta Sci. Agric. 5, 47–48 (2021).Article 

Google Scholar 
Ogbo, E. & Okhuoya, J. Bio-absorption of some heavy metals by Pleurotus tuber-regium Fr. Singer (an edible mushroom) from crude oil polluted soils amended with fertilizers and cellulosic wastes. Int. J. Soil Sci. 6, 34–48 (2011).Article 

Google Scholar 
Brunnert, H. & Zadražil, F. The translocation of mercury and cadmium into the fruiting bodies of six higher fungi. Eur. J. Appl. Microbiol. Biotechnol. 17, 358–364 (1983).Article 

Google Scholar 
Singh, V., Singh, M. P. & Mishra, V. Bioremediation of toxic metal ions from coal washery effluent. Desalin. Water Treat. 197, 300–318 (2020).Article 

Google Scholar 
Lazarova, N., Krumova, E., Stefanova, T., Georgieva, N. & Angelova, M. The oxidative stress response of the filamentous yeast Trichosporon cutaneum R57 to copper, cadmium and chromium exposure. Biotechnol. Biotechnol. Equip. 28, 855–862 (2014).Article 
PubMed 
PubMed Central 

Google Scholar  More

Study settingUganda covers 241,037 km2 (4° N, 2° S, 29° E, 35° E) and averages 1100 m above sea level25. Uganda has a total population of 34.6 million people, with 7.4 and 27.2 million living in urban and rural areas, respectively26. It borders Lake Victoria and has an equatorial climate. The study area is mainly plateau, with a few mountains. About 20 percent of Uganda’s area is covered by swamps and water bodies, including the four Great Lakes of East Africa (Lake Edward, Lake Victoria, Lake Albert, and Lake Kyoga). The country has ten national parks housing a high diversity of wildlife and endangered species26. By 2009, agriculture ranked as the second leading contributor to the country’s Gross Domestic Product27.Surveillance data of Bacillus anthracis infectionSurveillance data of livestock and human cases from 2018 were provided by the Ministry of Health (MOH) in Uganda through the Field Epidemiology and Laboratory Training Program and the Uganda Ministry of Agriculture, Animal Industry and Fisheries (MAAIF). Geographical coordinates of anthrax cases among wildlife (from 2004 to 2010) in Queen Elizabeth National Park were obtained from a recently published study4. These cases were mapped to show the distribution of anthrax across Uganda (Fig. 1). It is only recently that Uganda has mandated systematic anthrax surveillance across the country following the outbreak that started in 2018. GPS coordinates gathered by field personnel during outbreak responses were used to map outbreak locations. The human anthrax cases were defined based on the CDC’s clinical criteria (signs and symptoms), presumptive laboratory diagnosis (Gram staining), and confirmatory laboratory diagnosis (bacterial culture, immunohistochemistry, ELISA, and PCR)28. The animal anthrax cases were also defined based on the clinical presentation, presumptive laboratory diagnosis, and confirmatory laboratory diagnosis. All cases were classified as either ‘probable’ or ‘confirmed,’ with probable defined as cases that met both the clinical and presumptive laboratory diagnostic criteria and confirmed defined as cases that met the clinical and confirmatory laboratory diagnostic criteria. A total of 497 livestock (n = 171), humans (n = 32), and wildlife cases (n = 294), both confirmed (n = 32) and probable (n = 465), occurring from 2004 to 2018 were compiled. All methods were performed in accordance with the relevant guidelines and regulations.Figure 1Distribution of anthrax presence and pseudo-absence locations across Uganda. The navy-blue circles show wildlife cases (n = 294) used for model training, the blue triangles show livestock cases (n = 171), and the red diamonds represent human cases (n = 32). The blue polygons show the locations of the 50 km, 75 km, and 100 km buffers which were constructed around the wildlife cases with a distance of 10 km between the buffers and the presence locations. The pink dots show the pseudo-absence points selected within the 50 km buffer, the orange dots show the pseudo-absence points selected within the 75 km buffer, and the white dots show the pseudo-absence points selected within the 100 km buffer. Prediction maps were developed using the Quantum Geographic Information System software (QGIS). URL: https://qgis.osgeo.org (2020).Full size imageEnvironmental variable processingCorrelative studies of environmental risk factors for anthrax outbreaks suggest that temperature6,29,30,31,32,33,34,35,36,37, precipitation6,29,30,31,32,33,34,35,36,37,38,39, elevation6,29,31,32,34,35,37,38,39, soil (type, calcium concentration, pH, carbon content, and moisture)6,30,31,33,34,36,37,39,40,41,42,43, vegetation6,29,31,34,36,37,38,39,40,42,43, and hydrology37,44 are some of the major drivers of B. anthracis suitability. A total of 26 environmental predictors (Fig. 2) were selected for this study based on these known variables. These comprised 19 bioclimatic variables (the mean for the years 1970–2000) collected from the WorldClim database version 2 (https://www.worldclim.org/data/worldclim21.html) at a resolution of 30 s (~ 1 km2)45. Four soil variables, including exchangeable calcium at a depth of 0–20 cm, soil water availability, soil pH (10×) in H2O at a depth of 0 cm, and soil organic carbon at a depth of 0–5 cm, were retrieved from the International Soil Reference and Information Centre (ISRIC) data hub at a resolution of 250 m (https://data.isric.org/geonetwork/srv/eng/catalog.search#/home). Distance to permanent water bodies was derived from a global hydrology map provided by ArcGIS online version 10.6.146. Elevation data of 1 km2 in resolution was obtained from the Global Multi-resolution Terrain Elevation Data (GMTED2010) dataset available from the United States Geological Service. Finally, the monthly Enhanced Vegetation Index (EVI) data for the years 2004, 2005, and 2010 (36 tiles in total) were obtained from The Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MYD13A3 v.6) at a spatial resolution of 1 km (https://lpdaac.usgs.gov/products/myd13a3v006/). The single variable, mean EVI, was calculated in QGIS by averaging all 36 tiles. The EVI minimizes variations in the canopy background and maintains precision over conditions with dense vegetation.Figure 2Results of correlation between covariates using Pearson’s correlation test. Correlation between covariates was shown by red numbers (negative correlation) and blue numbers (positive correlation). BIO1 = Annual Mean Temperature, BIO2 = Mean Diurnal Range, BIO3 = Isothermality, BIO4 = Temperature Seasonality, BIO5 = Max Temperature of Warmest Month, BIO6 = Min Temperature of Coldest Month, BIO7 = Temperature Annual Range, BIO8 = Mean Temperature of Wettest Quarter, BIO9 = Mean Temperature of Driest Quarter, BIO10 = Mean Temperature of Warmest Quarter, BIO11 = Mean Temperature of Coldest Quarter, BIO12 = Annual Precipitation, BIO13 = Precipitation of Wettest Month, BIO14 = Precipitation of Driest Month, BIO15 = Precipitation Seasonality, BIO16 = Precipitation of Wettest Quarter, BIO17 = Precipitation of Driest Quarter, BIO18 = Precipitation of Warmest Quarter, BIO19 = Precipitation of Coldest Quarter.Full size imageAll environmental variables were resampled using the R package ‘Resample’47 to a resolution of 1 km and clipped to the extent of Uganda. Since data sampling for wildlife data (used in model training) was not done systematically across the study area, target backgrounds buffers (polygon buffers created at certain radii from the training points and used for the random selection of pseudo-absences) were created for model calibration to reduce sampling bias. As there was no information on the sampling radius, a sensitivity analysis was conducted by creating target backgrounds using circular buffers of radius 50 km, 75 km, and 100 km around the presence points used for model training, leaving 10 km between the presence points and the various buffers (Fig. 1). Quantum Global Information System (QGIS) version 3.16 (https://qgis.org) software was then used to add 294 random pseudo-absence points within each buffer polygon giving a ratio of 1:1 for the training presences to pseudo-absences. A recent study explored how four approaches of pseudo-absence creation affect the performance of models across different species and three model types by building both terrestrial and marine models using boosted regression trees, generalised additive mixed models, and generalised linear mixed models48. They then tested four methods for generating pseudoabsences across all the different model types: (1) correlated random walks (RW); (2) reverse correlated RW; (3) sampling pseudoabsences within a buffer area surrounding the presence points; (4) background sampling48. The findings of the study suggested that the separation or distance in the environmental space between the presence locations and the pseudoabsences was the most significant driver of the model predictive ability and explanatory power, and thus finding was consistent across the three model types (boosted regression trees, generalised additive mixed models, and generalised linear mixed models) and both the terrestrial and marine habitats48.The values of the environmental variables were then extracted for each presence and pseudo-absence location using the raster package in R. With these, we did an initial data exploration to check for outliers within the covariates, collinearity, and to explore the relationships between the covariates and the response variables (presence or absence of anthrax). Cleveland dot plots were used to check for possible outliers. Following the outlier checks, variance inflation factors (VIF), pairwise plots, and Pearson correlation coefficients with correction for multiple comparisons were used to measure the statistically significant correlation between the covariates (Fig. 2). For variables that were highly correlated (correlation greater than 0.6) or those with high variance inflation (VIF  > 3), only one was used in the modelling process. Five variables were selected following this analysis: Temperature seasonality (BIO4), elevation, distance to water, soil calcium, and soil water (Table 1).Table 1 Summary of the environmental variables used.Full size tableModelling anthrax suitability across UgandaQGIS v.3.16 (https://qgis.org) and the R statistical package version 4.1.049 were used to conduct data visualization, cleaning, and model analysis (R code used available in a Github repository: https://github.com/valentinandolo/Uganda-Spatial_Model/tree/master). The wildlife cases (n = 294) were used for model training and testing. The remaining human (n = 32) and livestock (n = 171) cases were used for model evaluation. Since the wildlife case locations were recorded from 2004 to 2010, while the livestock and human cases were recorded in 2018, the latter locations were both spatially and temporally distinct from the wildlife cases, making an excellent basis for block cross-validation of the final model performance50. Random partitioning of the data into training and testing sets can inflate the performance of a model and underestimate the error in the spatial prediction evaluation50. Block cross-validation uses spatial blocks that separate the testing and training datasets; thus, the method has been suggested to be a good technique for error estimation and a robust approach for measuring a model’s predictive performance50.The INLA package in R was applied to model the suitability of B. anthracis across Uganda. INLA calculates the spatial interaction effects using a Stochastic Partial Differential Equation (SPDE) method, which estimates a continuous Gaussian Markov Random Field (GMRF) where the correlation between two locations in space is specified by the Matérn correlation which is explained in more detail elsewhere51. The initial step in fitting an SPDE model is the creation of a Constrained Refined Delaunay Triangulation or a mesh to illustrate the spatial process51. R-INLA uses a function called ‘inla.mesh.2d()’ that applies a variety of arguments to build the mesh, these include: loc, loc.domain, boundary, max.edge, and cutoff51. The loc argument contains the point locations which provide information about the area of study and are used to create the triangulation nodes. Alternatively, a polygon of the study area can be used to identify the extent of the domain via loc.domain. We applied the point locations using the loc argument. We then specified the boundary of the mesh as a convex hull. We used the max.edge argument to specify the maximum edge length for the inner mesh domain/triangles and the outer triangles. We did this by first studying the distribution of distances between the point locations for the training presences and pseudo-absences. Most points were within a distance of about 90–100 km away from each other, thus, a possible guess for the range at which spatial autocorrelation persists was 100 km. We used a distance of 20 km as a range guess to create a finer mesh which has been shown to produce more precise models. We specified the maximum edge length for the inner triangles as 20 km divided by 5 (4 km) and the maximum edge length for the outer triangles as 20 km. Finally, the cutoff argument sets the minimum distance allowed between point locations. We divided the maximum edge length for the inner triangles by 5 (4 km divided by 5 = 0.8 km), meaning that points that were closer in distance than 0.8 km were replaced by one vertex to avoid the occurrence of small triangles.The spatial effect, which is a numeric vector, then links each observation within the data to a spatial location, thus, accounting for region-specific variation that cannot be accounted for by the covariates. Following the recommendation by Lindgren and Rue52, multivariate Gaussian distributions with means of zero and a spatially defined covariance matrix were used to model the spatial effect. Several versions of Bayesian hierarchical additive models were created by estimating a Bernoulli generalized additive model (GAM) with and without spatially correlated random effects. The Bernoulli GAM is defined as shown in Eqs. (1) and (2)$${C}_{i} sim Bernoulli left({p}_{i}right),$$
(1)
$$logit left({p}_{i}right)= alpha +sum_{j=1}^{m}{beta }_{j}left({X}_{j,i}right)+ sum_{k=1}^{l}{f}_{k}left({X}_{k,i}right)+ {mu }_{i},$$
(2)
where ({C}_{i}) denotes the observed value, such that: B. anthracis presence or absence at a given location i (i = 1, …, n; n = 588) is given as ({C}_{i}), where ({C}_{i}) =1 if B. anthracis was present, and ({C}_{i})=0 if absent. Logit is the link function for binomial family, ({p}_{i}) is the expected value of the response variable (the probability of B. anthracis suitability) at location i, α is the intercept, ({X}_{j,i}) and ({X}_{k,i}) are the j th and the k th covariates at a location i, ({beta }_{j}) are the beta coefficients, ({f}_{k}) are the smooth functions (cubic regression splines) for k th covariates, and ({mu }_{i}) is the spatial random effect at the location i53. The number of variables in linear term (m) and the non-liner term (l) are different because the variables employed in each term are different. We estimated both linear and non-linear effects for the covariates. Our overall database had 294 B. anthracis pseudo-absences generated randomly across the target background buffers to match the number of species presences recorded. Because we had no prior information, a Gaussian prior distribution with a mean of zero (default no effect prior unless data is informative) was applied for all the model parameters. The posterior mean, standard deviation, and 95% credible intervals were estimated for all the parameters.Model selectionSeveral different candidate models were examined. First, a baseline model was built using only the intercept. A second baseline model was then built, which included the intercept and spatial random effects only. Covariates were added to the second baseline model (intercept plus spatial effects model) without any smoothing function (i.e., only linear effects). The contribution of the spatial random effect was then re-examined by taking it out from the model. Smoothing functions were then added to all covariates, and the same procedure was repeated. Model selection was done using this forward stepwise approach. The final model was run using the three different target background buffers to identify the buffer distance with the lowest Deviance Information Criterion (DIC). The various options were assessed using the DIC54, Watanabe-Akaike information criterion (WAIC), and the Conditional Predictive Ordinate (CPO). The DIC and WAIC were chosen because they are commonly used to assess model performance by measuring the compromise between the goodness of fit and complexity. For CPO, the logarithmic score (− mean(log (CPO)) (LCPO) was calculated and used55. CPO can also be used to conduct internal cross-validation of models using a leave-one-out approach to evaluate the predictive performance of the model. Lower LCPO, DIC, and WAIC estimates suggest superior model performance. Thus, the favoured model had the lowest values across the 3 metrics.Model validation and evaluationModel validation for the favoured model was conducted using an independent evaluation dataset comprising of livestock and human outbreaks occurring at different spatial locations and 8 years after the training data. The omission rate, which indicates the proportion of positive test locations that end up in pixels predicted to be unsuitable for B. anthracis56, was used to validate the model. A low omission rate indicates good model performance. The threshold for suitability was the probability threshold that maximized the sensitivity and specificity of the model. The sensitivity was then derived by calculating the proportion of positive test locations that end up in pixels predicted to be suitable for B. anthracis56. A high sensitivity indicates good model performance.Model predictionThe favoured model selected using the criteria described above was used to generate countrywide prediction maps showing the posterior mean values, standard deviation, and the 95% credible intervals of the probability of B. anthracis suitability. The raster package in R was used to make the prediction maps. Bayesian kriging was done by treating all model parameters as random variables to include uncertainty in the prediction57. This kriging is built into the INLA framework via the SPDE, which allows a Delaunay triangulation to be constructed around the presence and absence locations within the sampling frame52. INLA then conducts model inference and prediction at the same time by considering the prediction points as locations missing the response variable (set to NA)51. Following successful model prediction, additional linear interpolation functions then generate the output for the whole study area scaled from 0 to 1.Ethical approvalEthical approval for this study was provided by the Human Biology Research Ethics Committee, University of Cambridge, UK (Ref: HBREC.2019.02) the School of Veterinary Medicine and Animal Resources Institutional Animal Care and Use Committee, Makerere University, Uganda (Ref: SVAREC/21/2019); and the Uganda National Council for Science and Technology. Informed consent was obtained from all subjects and/or their legal guardian(s). All methods were performed in accordance with the relevant guidelines and regulations. More

Chironomids have the ability to survive and reproduce in polluted environments, and thus they are included in many ecological studies where approaches may be taxonomic or functional16. The diversity of most macroinvertebrates is controlled by the oxygen level of water, but chironomids may survive in hypoxic conditions where the oxygen concentration may be less than 3 mg l−117. The current study demonstrates that changing seasons, as well as anthropogenic activities, have a significant impact on the levels of DO in aquatic bodies. As observed from the result, DO highly influences the tube length of the chironomid larvae. Since KWC is a wastewater canal, the average oxygen level is lower (5.24 ± 1.14 mg l−1) than KFP (6.63 ± 1.28 mg l−1) which is a normal fish culturing pond. It has also been observed that the average tube length of the chironomid larvae of KWC (8.66 ± 0.88 mm) is higher than KFP (7.68 ± 0.62 mm), which indicates that a low concentration of DO promotes the building of longer tubes in natural conditions. Similar observations were also observed in laboratory conditions. When the oxygen level (7.03 ± 0.41 mg l−1) in the experiment was kept in the normal range, there was negligible variation in tube length (7.61 ± 0.31 mm). But when the concentration of oxygen is gradually reduced by dilution, the tube length starts to increase accordingly, which is explained graphically in Fig. 4. The regression model of both the experimental conditions also supports the hypothesis that the tube length has an inverse relationship with DO. The scatter plot and simple linear regression confirmed the inverse relationship between DO and tube length (Figs. 1 and 2).Chironomid larvae are able to grow in the polluted water of a wastewater pond as dominant macroinvertebrates18. It is observed that those larvae living in the sand tubes are more susceptible to chemical pollutants than the larvae living in silt tubes7. Sand particles are bigger than silt and are not suitable for the survival of larvae19. Chironomus riparius larvae make their tubes from different external particles and their own proteins20. Midge larvae are the inhabitants of sediments, and at the same time, sediment is the depository of different inorganic, organic, and heavy metals. In such cases, the tube of chironomid larvae may act as a defensive structure, which protects them from the adverse effects of undesirable pollutants and may increase their tolerance against such chemicals21,22,23.Larvae can thrive in benthic sediments with high decaying organic content and very low DO concentrations in water bodies24. In poor DO concentration, larvae can survive due to the presence of haemoglobin in their body tissue fluid, which plays an important physiological role in increasing respiratory efficiency, as was observed in Chironomus plumosus. Longer tube length may help larvae generate better respiratory currents so that they can cope with a low DO environment.Tube length is crucial for living in water because primarily tubes protect them from outer environmental factors like predators, and pollution. It was observed during this study that when the DO of water is low, larvae make elongated tubes to reach the upper layer of water, where the DO level is comparatively high. To get their required amount of oxygen, the larvae increase the tube length towards the water surface and increase the DO in tube water by undulating the body and other structures, creating a current inside the tube25,26. On contrary, when the DO level of the surrounding water of chironomid is sufficient, they can manage their normal physiological activities with the available oxygen. They need not to elongate their tube length. That’s why their tube length is inversely related to the DO of their surrounding medium.If tube length does not increase in size in hypoxic water, larvae will not be able to meet their oxygen demand. If the DO of water decreases, tube length will increase and vice versa. Behavioural and physiological adaptations of chironomids larvae make them successful to live in a hypoxic environment. Thus, in hypoxic conditions, larvae with longer tubes are able to gather more oxygen from the upper layer of water and get more space to create a current of water to increase the amount of O2 inside the tube by undulating the preanal papillae, anal gill, ventral gills. This would explain why the tube length of chironomids depends on the DO of water. Hence by measuring the tube length with a standard measuring scale, one may get an idea about the quality of water, especially DO, before doing any chemical analysis. The work seems to be unique and novel for its own kind. More

Chen, I. C., Hill, J. K., Ohlemüller, R., Roy, D. B. & Thomas, C. D. Rapid range shifts of species associated with high levels of climate warming. Science 333(6045), 1024–1026. https://doi.org/10.1126/SCIENCE.1206432 (2011).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Newton, I. Migration within the annual cycle: Species, sex and age differences. J. Ornithol. 152, 169–185. https://doi.org/10.1007/S10336-011-0689-Y/TABLES/1 (2011).Article 

Google Scholar 
Dodson, S., Abrahms, B., Bograd, S. J., Fiechter, J. & Hazen, E. L. Disentangling the biotic and abiotic drivers of emergent migratory behavior using individual-based models. Ecol. Model. 432, 109225. https://doi.org/10.1016/J.ECOLMODEL.2020.109225 (2020).Article 

Google Scholar 
Lehikoinen, A. et al. Sex-specific timing of autumn migration in birds: the role of sexual size dimorphism, migration distance and differences in breeding investment. Ornis Fennica 94, 53–65 (2017).
Google Scholar 
Stewart, B. S. Ontogeny of differential migration and sexual segregation in northern elephant seals. J. Mammol. 78(4), 1101–1116 (1997).Somveille, M., Rodrigues, A. S. L. & Manica, A. Why do birds migrate? A macroecological perspective. Glob. Ecol. Biogeogr. 24(6), 664–674. https://doi.org/10.1111/geb.12298 (2015).Article 

Google Scholar 
Corkeron, P. J. & Connor, R. C. Why do baleen whales migrate?. Mar. Mamm. Sci. 15(4), 1228–1245. https://doi.org/10.1111/J.1748-7692.1999.TB00887.X (1999).Article 

Google Scholar 
Mourier, J., Mills, S. C. & Planes, S. Population structure, spatial distribution and life-history traits of blacktip reef sharks Carcharhinus melanopterus. J. Fish Biol. 82(3), 979–993. https://doi.org/10.1111/JFB.12039 (2013).Article 
CAS 
PubMed 

Google Scholar 
Avgar, T., Mosser, A., Brown, G. S. & Fryxell, J. M. Environmental and individual drivers of animal movement patterns across a wide geographical gradient. J. Anim. Ecol. 82, 96–106. https://doi.org/10.1111/j.1365-2656.2012.02035.x (2013).Article 
PubMed 

Google Scholar 
Crawshaw, L. I. Physiological and behavioral reactions of fishes to temperature change. J. Fish. Res. Board Can. 34(5), 730–734. https://doi.org/10.1139/f77-113 (1977).Article 

Google Scholar 
Heithaus, M., Dill, L., Marshall, G. J. & Buhleier, B. Habitat use and foraging behavior of tiger sharks (Galeocerdo cuvier) in a seagrass ecosystem. Mar. Biol. 140, 337–348. https://doi.org/10.1007/s00227-001-0711-7 (2002).Article 

Google Scholar 
Magnuson, J. J., Crowder, L. B. & Medvick, P. A. Temperature as an ecological resource. Integr. Comp. Biol. 19(1), 331–343. https://doi.org/10.1093/icb/19.1.331 (1979).Article 

Google Scholar 
Matern, S. A., Cech, J. J. & Hopkins, T. E. Diel movements of bat rays, Myliobatis californica, in Tomales Bay, California: Evidence for behavioral thermoregulation?. Environ. Biol. Fishes 58(2), 173–182. https://doi.org/10.1023/A:1007625212099 (2000).Article 

Google Scholar 
Speed, C. W., Meekan, M. G., Field, I. C., McMahon, C. R. & Bradshaw, C. J. A. Heat-seeking sharks: Support for behavioural thermoregulation in reef sharks. Mar. Ecol. Prog. Ser. 463, 231–244. https://doi.org/10.3354/meps09864 (2012).Article 
ADS 

Google Scholar 
Dewar, H., Domeier, M. & Nasby-Lucas, N. Insights into young of the year white shark, Carcharodon carcharias, behavior in the Southern California Bight. Environ. Biol. Fishes https://doi.org/10.1023/B:EBFI.0000029343.54027.6a.pdf (2004).Article 

Google Scholar 
Hertz, P. E., Huey, R. & Stevenson, R. D. Evaluating temperature regulation by field-active ectotherms. Am. Nat. 142, 796–818 (1993).Article 
CAS 
PubMed 

Google Scholar 
Heupel, M. R., Simpfendorfer, C. A. & Hueter, R. E. Estimation of shark home ranges using passive monitoring techniques. Environ. Biol. Fishes 71(2), 135–142. https://doi.org/10.1023/b:ebfi.0000045710.18997.f7 (2004).Article 

Google Scholar 
Topping, D. T., Lowe, C. G. & Caselle, J. E. Site fidelity and seasonal movement patterns of adult California sheephead Semicossyphus pulcher (Labridae): An acoustic monitoring study. Mar. Ecol. Progr. Ser. 326, 257–267 (2006).Weng, K. C. et al. Movements, behavior and habitat preferences of juvenile white sharks Carcharodon carcharias in the eastern Pacific. Mar. Ecol. Prog. Ser. 338, 211–224. https://doi.org/10.3354/meps338211 (2007).Article 
ADS 

Google Scholar 
Lyons, K. et al. The degree and result of gillnet fishery interactions with juvenile white sharks in southern California assessed by fishery-independent and -dependent methods. Fish. Res. 147, 370–380. https://doi.org/10.1016/J.FISHRES.2013.07.009 (2013).Article 
ADS 

Google Scholar 
Papastamatiou, Y. P. et al. Drivers of daily routines in an ectothermic marine predator: Hunt warm, rest warmer?. PLoS ONE. https://doi.org/10.1371/journal.pone.0127807 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Adolph, S. C. Influence of behavioral thermoregulation on microhabitat use by two sceloporus lizards. Ecology 71(1), 315–327. https://doi.org/10.2307/1940271 (1990).Article 

Google Scholar 
Heithaus, M. R. The biology of tiger sharks, Galeocerdo cuvier, in Shark Bay, Western Australia: sex ratio, size distribution, diet, and seasonal changes in catch rates. Environ. Biol. Fishes 61, 25–36 (2001).Article 

Google Scholar 
Vaudo, J. J. & Lowe, C. G. Movement patterns of the round stingray Urobatis halleri(Cooper) near a thermal outfall. J. Fish Biol. 68(6), 1756–1766. https://doi.org/10.1111/j.0022-1112.2006.01054.x (2006).Article 

Google Scholar 
Vaudo, J. J. & Heithaus, M. R. Microhabitat selection by marine mesoconsumers in a thermally heterogeneous habitat: Behavioral thermoregulation or avoiding predation risk?. PLoS ONE. 8(4), e61907. https://doi.org/10.1371/journal.pone.0061907 (2013).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Weng, K. C. et al. Migration and habitat of white sharks (Carcharodon carcharias) in the eastern Pacific Ocean. Mar. Biol. 152(4), 877–894. https://doi.org/10.1007/s00227-007-0739-4 (2007).Article 

Google Scholar 
White, C. F. et al. Quantifying habitat selection and variability in habitat suitability for juvenile white sharks. PLoS ONE 14(5), e0214642. https://doi.org/10.1371/journal.pone.0214642 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Curtis, T. H. et al. First insights into the movements of young-of-the-year white sharks (Carcharodon carcharias) in the western North Atlantic Ocean. Sci. Rep. 8(1), 1–8. https://doi.org/10.1038/s41598-018-29180-5 (2018).Article 
ADS 
CAS 

Google Scholar 
Bruce, B. D., Harasti, D., Lee, K., Gallen, C. & Bradford, R. Broad-scale movements of juvenile white sharks Carcharodon carcharias in eastern Australia from acoustic and satellite telemetry. Mar. Ecol. Prog. Ser. 619, 1–15. https://doi.org/10.3354/MEPS12969 (2019).Article 
ADS 

Google Scholar 
Carey, F. G. et al. Temperature and activities of a white shark Carcharodon carcharias. Copeia 2, 254–260. https://doi.org/10.2307/1444603 (1982).Article 

Google Scholar 
Klimley, A. P., Beavers, S. C., Curtis, T. H. & Jorgensen, S. J. Movements and swimming behavior of three species of sharks in La Jolla Canyon, California. Environ. Biol. Fish. 63, 117–135. https://doi.org/10.1023/A:1014200301213.pdf (2002).Article 

Google Scholar 
Towner, A. V., Underhill, L. G., Jewell, O. J. D. & Smale, M. J. Environmental Influences on the abundance and sexual composition of white sharks Carcharodon carcharias in Gansbaai, South Africa. PLoS ONE. 8(8), e71197. https://doi.org/10.1371/journal.pone.0071197 (2013).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Anderson, J. M. et al. High-resolution acoustic telemetry reveals swim speeds and inferred field metabolic rates in juvenile white sharks (Carcharodon carcharias). PLoS ONE 17(6), e0268914. https://doi.org/10.1371/JOURNAL.PONE.0268914 (2022).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Anderson, J. M. et al. Interannual nearshore habitat use of young of the year white sharks off Southern California. Front. Mar. Sci. 8, 238. https://doi.org/10.3389/fmars.2021.645142 (2021).Article 

Google Scholar 
Domeier, M. L. & Nasby-Lucas, N. Two-year migration of adult female white sharks (Carcharodon carcharias) reveals widely separated nursery areas and conservation concerns. Anim. Biotelemet. 1(1), 1–10. https://doi.org/10.1186/2050-3385-1-2/FIGURES/3 (2013).Article 

Google Scholar 
Oñate-González, E. C. et al. Importance of Bahia Sebastian Vizcaino as a nursery area for white sharks (Carcharodon carcharias) in the Northeastern Pacific: A fishery dependent analysis. Fish. Res. 188, 125–137. https://doi.org/10.1016/J.FISHRES.2016.12.014 (2017).Article 

Google Scholar 
Lowe, C. G. et al. Historic fishery interactions with white sharks in the Southern California Bight. Glob. Perspect. Biol. Life Hist. White Shark 14, 169–190 (2012).
Google Scholar 
Anderson, J. M. et al. Non-random Co-occurrence of Juvenile White Sharks (Carcharodon carcharias) at Seasonal Aggregation Sites in Southern California. Front. Mar. Sci. 8, 1–14. https://doi.org/10.3389/fmars.2021.688505 (2021).Article 
ADS 
CAS 

Google Scholar 
Benson, J. F. et al. Juvenile survival, competing risks, and spatial variation in mortality risk of a marine apex predator. J. Appl. Ecol. 55, 2888–2897. https://doi.org/10.1111/1365-2664.13158 (2018).Article 

Google Scholar 
RStudio Team. RStudio: Integrated Development for R. (RStudio, PBC, 2020) http://www.rstudio.com/.Derrick, T., & Thomas, J. Time Series Analysis: The Cross-Correlation Function. Innovative Analyses of Human Movement, Chapter 7. https://lib.dr.iastate.edu/kin_pubs/46 (2004).Killick, R., Fearnhead, P. & Eckley, I. A. Optimal detection of changepoints with a linear computational cost. J. Am. Stat. Assoc. 107, 1590–1598. https://doi.org/10.1080/01621459.2012.737745 (2012).Article 
MathSciNet 
CAS 
MATH 

Google Scholar 
Bakun, A. Coastal Upwelling Indices, West Coast of North America. US Department of Commerce. NOAA Technical Report, NMFS SSRF-671 (1973).Di Lorenzo, E. Seasonal dynamics of the surface circulation in the Southern California Current System. Deep-Sea Res. Part II 50(14–16), 2371–2388. https://doi.org/10.1016/S0967-0645(03)00125-5 (2003).Article 
ADS 

Google Scholar 
Lynn, R. J. & Simpson, J. J. The California Current System: The seasonal variability of its physical characteristics. J. Geophys. Res. 92(C12), 12947. https://doi.org/10.1029/jc092ic12p12947 (1987).Article 
ADS 

Google Scholar 
Sinnett, G. & Feddersen, F. The surf zone heat budget: The effect of wave heating. Geophys. Res. Lett. 41(20), 7217–7226. https://doi.org/10.1002/2014GL061398 (2014).Article 
ADS 

Google Scholar 
Wei, X., Li, K.-Y., Kilpatrick, T., Wang, M. & Xie, S.-P. Large-scale conditions for the record-setting Southern California marine heatwave of August 2018. Geophys. Res. Lett. 48(7), e2020GL091803 (2021).Article 
ADS 

Google Scholar 
Freedman, R. M., Brown, J. A., Caldow, C. & Caselle, J. E. Marine protected areas do not prevent marine heatwave-induced fish community structure changes in a temperate transition zone. Sci. Rep. 10(1), 1–8. https://doi.org/10.1038/s41598-020-77885-3 (2020).Article 
CAS 

Google Scholar 
Heupel, M. R., Simpfendorfer, C. A. & Hueter, R. E. Running before the storm: blacktip sharks respond to falling barometric pressure associated with Tropical Storm Gabrielle. J. Fish Biol. 63(5), 1357–1363. https://doi.org/10.1046/J.1095-8649.2003.00250.X (2003).Article 

Google Scholar 
Guttridge, T. L. et al. Deep danger: Intra-specific predation risk influences habitat use and aggregation formation of juvenile lemon sharks Negaprion brevirostris. Mar. Ecol. Progr. Ser. 445, 279–291 (2012).Article 
ADS 

Google Scholar 
Grainger, R. et al. Diet composition and nutritional niche breadth variability in juvenile white sharks (Carcharodon carcharias). Front. Mar. Sci. 7, 422 (2020).Article 

Google Scholar 
Hussey, N. E., Christiansen, H. M. & Dudley, S. F. J. Size-based analysis of diet and trophic position of the white shark, carcharodon carcharias, in South African waters. Glob. Perspect. Biol. Life Hist. White Shark 3, 27–49. https://doi.org/10.1201/b11532-5 (2012).Article 

Google Scholar 
Kim, S. L., Tinker, M. T., Estes, J. A. & Koch, P. L. Ontogenetic and among-individual variation in foraging strategies of northeast Pacific white sharks based on stable isotope analysis. PLoS ONE 7(9), e45068. https://doi.org/10.1371/JOURNAL.PONE.0045068 (2012).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Tinker, M. T. et al. Dramatic increase in sea otter mortality from white sharks in California. Mar. Mamm. Sci. 32(1), 309–326. https://doi.org/10.1111/mms.12261 (2015).Article 

Google Scholar  More

Sample collectionA survey was conducted among the residents of nine counties in Kenya (Mombasa, Kisumu, Machakos, Nairobi, Makueni, Bomet, Kericho, Kiambu, and Narok) and GPS location coordinates were recorded and later used to build the predictive model (“Infestation dynamics of bedbugs in residential communities” section). These counties represent diversity in cultural practices, livelihood strategies (such as fishing, tourism, farming), and infrastructure development. Also, they comprise different altitudes above sea level, temperatures, and differing in average annual rainfall.Samples identification using morphological identification keysIn each county where the survey was conducted, bedbug samples was taken and preserved in ethanol 70% for morphological identification. Cimex belonging to Cimicidae family is the common genus adapted to human environment and reported throughout the world and comprising species such as Cimex lectularius and C. hemipterus that are hematophagous mainly feeding on human blood5. The key morphological features used in identifying bedbugs include: (1) the head has a labrum that appears as a free sclerite at the extreme anterior margin, ecdysial lines form a broad V, eyes project from the sides composed of several facets and the antennae are 4-segmented, (2) thorax is subdivided into prothorax, mesothorax and metathorax, (3) legs have all other normal parts except pulvilli and arolia, tarsus is 3-segmented with 2 simple claws, (4) the abdomen has 11 more-or-less segmented recognizable segments, 7 pairs of spiracles borne on the second to eighth segments, hosts the genital structures, paramere in males and mesospermalege in females45. Bedbug specimen morphological features were examined using Leica EZ24 HD dissecting microscope (Leica Microsystems, UK) and photos documented using the associated software.Survey for household’s knowledge and perceptions on bedbugsThis study was a community-based cross-sectional survey conducted from November–December 2020 with respect of the rules/guidlines introduced by the Ministry of Health to contain the COVID-19 pandemic in Kenya (wearing mask, social distance, washing hand, etc.). It was based on a stratified, systematic random sampling where 100 respondents were selected from each county.A total number of 900 respondents were randomly selected and the household head or the representative showing willingness and consent was interviewed face-to-face. The interview was conducted using a semi-structured questionnaire prepared in the English language (Appendix A). The questionnaire was translated into the local native language (Kiswahili) to avoid biasness and improve the understanding between the enumerator and the respondent. Prior to the commencement of the survey and authentic data collection, a pre-testing exercise was performed by training enumerators on a similar socio-demographic pattern. This was useful for improving the quality of data, ensuring validity, familiarizing the enumerators with the questionnaire, and data handling.The information collected using the semi-structured questionnaire included residents’ socio-economic profiles, knowledge, and perceptions on the pest, bedbug incidence, and management practices. The socio-economic profile factors addressed in the survey comprised gender, age, education, access to basic social amenities, and household size. The study also prioritized the financial consequences, the severity of the bites, perceptions of respondents on the pest, and management practices for its control.Survey data were checked for errors, completeness, summarized, and entered in Microsoft-Excel. It was then cleaned and transferred to Statistical Package for Social Science (SPSS) version 25 software (IBM Corp., Armonk, NY) for purposes of descriptive statistics (means and percentages).In contrast, in instances where more than one reason was given for a single question, percentages were calculated based on each group of similar responses. Chi-square was performed to determine the differences regarding socio-demographic characteristics, knowledge, and perceptions on bedbugs and control practices. Additionally, data were disaggregated by gender and age categories to understand the existing differences among the various respondent categories. Besides, F-test statistics was performed on the ages of respondents to determine the mean, standard deviation and statistical significance. The level of significance was considered when the p-value was below 5%.Infestation dynamics model of bedbugModel simulation assumptionsHouses infestation dynamics was studied following Susceptible-Infested-Treatment (SIT) model46. Therefore, houses in the community are classified into three groups: susceptible, infested or treated. Within a house, bedbug population dynamics was ignored, while it was considered from one house to another where infested houses have some potential to spread the infestation to other houses in the community. A population of bedbugs in an infested house has some probability per unit of time of becoming extinct either naturally or after treatment. In the infestation dynamics, the rate of house infestation depends on the number of infested houses, the movement of people from one house to another and the proportion of treated houses in the community. We assume that infested houses (I) spread the infestation at the rate β and only a fraction S/N of the houses is susceptible (S) to infestation. Infested houses become extinct at a certain rate known as rate γ. Infested houses are treated at the rate τ and the protection conferred is lost at the rate α. Ordinary differential equation developed to study SIT model were used in this study46. All the models used have the generic formulations displayed below:$$frac{dS}{dt}=frac{beta }{N}SI+gamma I+alpha T$$
(1)
$$frac{dI}{dt}=frac{beta }{N}SI-(gamma +tau )I$$
(2)
$$frac{dT}{dt}=tau I-alpha T$$
(3)
where β  > 0, τ  > 0, α ≥ 0 and γ  > 0. The total population size is N = S(t) + I(t) + t(t). The initial conditions satisfy at S(0)  > 0, I(0)  > 0, T(0) ≥ 0 and S(0) + I(0) = N, where N is the constant total population size, dN/dt = 0.Infestation dynamics models implementationThe method used to implement the infestation dynamics model of the pest is based on the system thinking approach with its archetypes [Causal Loop Diagram (CLD), Reinforcing (R) and Balancing (B)] by a mental and holistic conceptual framework. This is important for mapping how the variables, issues, and processes influence each other in the complex interactions of bedbugs within and between houses and their impacts. Despite these archetypes being qualitative, they are necessary for elucidating and disclosing the basic feedback configurations that occur in houses and their environs when infested with pests like bedbugs. A dynamic model was generated by converting the causal loop diagram (CLD) obtained using stocks, flows, auxiliary links, and clouds. Consequently, these in turn were translated into coupled differential equations for simulations.The SIT model was translated into causal loop diagram where arrows show the cause-effect relations where positive sign indicates direct proportionality of cause and effect while negative sign shows inverse proportionality relations, and two different scenarios have been assessed: (1) homogeneous houses where there is a single community of houses of the same quality, and (2) heterogeneous houses where there is a community of good and bad houses. Ancient houses presenting slits/fissures with less cleanliness and filled with old or secondhand furniture at low grade are considered bad houses as they may sustain high level of bedbug infestation; and new houses don’t provide well enough conditions for bedbug population to survive, and they are called in the model good houses47. Bad houses are considered to act as sources while good houses act as sinks, but all together are randomly distributed where each house has the same probability to contact good or bad houses.In the scenarios of homogeneous houses, the causal loop diagram (Fig. 7) has two feedback loops: (a) one positive, as the number of infested houses increases, the probability to get susceptible houses infested also increases resulting in infested houses increase; (b) one negative, as the infested houses increases, the treated houses increase resulting in susceptible houses decrease. The causal loop diagram is displayed in Fig. 7A while Fig. 7B showed the stocked and flows diagram and axillary variables obtained from causal loop diagram.Figure 7Susceptible-Infested-Treatment (SIT) model translated into causal loop diagram (A) and stock and flow diagram (B) for homogeneous houses and causal loop diagram (C) and stock and flow diagram (D) for heterogeneous houses in the community.Full size imageSusceptible, infested, and treated houses are stocks in the system, representing the number of houses susceptible, infested, and treated, respectively at a given point of time. The rates represent in and out-flows of the diagram. Auxiliary and constants that drive the behavior of the system were connected using information arrows within them and flows and stocks to represent the relations among variables in terms of equations.In the scenarios of heterogeneous houses, the causal loop diagram (Fig. 7C) comes with the two previous feedback loops but for each category of house. In addition, there is a fifth feedback loop that connect bad house to good house and vice versa.Therefore, as the infested bad houses increase, the probability to infest good houses increases. The more they are exposed the more they get infested. In turn, as the infested good houses increase, the chance to infest susceptible bad houses increases and the more they are exposed, the more they get infested, resulting in the increase of infested bad houses. The stocks and flows diagram of each of the two categories of houses occurred with interconnexion relationships between the two categories (Fig. 7D).Models’ simulationsThe survey data (“Bedbug Genus identification” section) on prevalence, knowledge, perceptions and self-reported; in addition, the respondents’ reported control mechanisms and their average time of effectiveness (Appendix B, Table S1) were used for model simulations. The different control methods reported were reclassified in three control approaches: chemical control, other control methods (including exposure to direct sunlight, use of hot water, painting, application of diesel, paraffin and wood ash, use of Aloe Vera extract and Herbs), and combination of chemical and other control methods. All the models commodities and units were checked before performing the simulations. Simulation and implementation of the models were done using Vensim PLP 8.1 platform (Ventana systems, Harvard, USA). It consists of a graphical environment that usually permits drawing of Causal Loop Diagram (CLD), stocks, flow diagrams and to carry out simulations. After we simulated the infestation dynamics under the two scenarios, we explored the effect of the different control methods.Spatial distribution analysis of bedbugs using MaxEnt modelEnvironmental data for MaxEntThe environmental variables used as the other maxent input were obtained by deriving bioclimatic, land cover, and elevation data. Bioclimatic variables and elevation (Digital Elevation Model; DEM) data were obtained from the Global Climate Data official website, Worldclim (http://www.worldclim.org/bioclim.htm)48 including 19 bioclimatic variables (Appendix B, Table S2). The land cover data were downloaded from the Global Land Cover Facility (GLCF).In order to reduce collinearity between predictors, a collinearity test was performed on all the variables by filtering them according to the following steps36: firstly, the MaxEnt model was run using the distribution data of bedbugs and 19 bioclimatic variables to obtain the percent contribution of each variable to the preliminary prediction results. Secondly, following the generation of the percentage contribution of all the variables, we then imported all distribution points in Arc-GIS and extracted the attribute values of the 19 variables. Furthermore, the “virtual species” package49 in R-software (R Foundation for Statistical Computing, Vienna, Australia) was used to explore the extracted variables’ clusters spatial correlation using Pearson’s correlation coefficient and the cluster tree (Fig. 8). Thus, the final number of predictor variables after screening was 5 establishing the potential geographical distribution of bedbug, which includes Temperature Seasonality (bio4), Precipitation of Driest Month (bio14), Temperature Annual Range (bio7), Precipitation of Driest Quarter (bio17) and Precipitation of Warmest Quarter (bio18) (Appendix B, Table S2). The land cover was considered because studies have shown its importance on insect spatial distribution50,51,52 and it was setled as a categorical variable53. Elevation was selected as variable because it greatly influences species’ occurrence and dispersal by affecting the temperature, precipitation, vegetation, and sun characteristics (direction, intensity, etc.) on the earth’s surface54,55,56. The study variables had different resolutions and were therefore, resampled to 1 km. The variables were clipped to Kenya and Africa boundaries and converted to ASCII (Stands for “American Standard Code for Information Interchange”) format using the ‘raster’ package49 in R statistical software (R Foundation for Statistical Computing, Vienna, Australia).Figure 8Key model predictor variables.Full size imageDistribution modelling in Kenya and AfricaIn our study, we used the maximum entropy distribution modelling method. This is because it has been recommended to have the ability to perform best and remain effective despite the use of small sample size relative to the other modelling methods57.Our selected bioclimatic variables (5) and occurrence/prevalence data for bedbugs were then imported into MaxEnt model and the options of ‘Create response curves’ and ‘Do jackknife’ were selected to measure variable importance’ options. The model output file was selected as ‘Logistic’, the commonly used approach is the random portioning of distribution datasets into ‘training’, and ‘test’ sets57,58. MaxEnt model was run with a total number of 5000 iterations and five replicates for better convergence of the model and rescaled within the range of 0–1000 suitability scores using ‘raster’ package49 in R statistical software (R Foundation for Statistical Computing, Vienna, Australia).The modelling performance/MaxEnt accuracy was evaluated by choosing the area under the receiver operating characteristics (ROC) curve (AUC) as the estimation index. This was important for the calibration and validation of the robustness of MaxEnt model evaluation. Furthermore, the area under the ROC curve (AUC) was necessary as an additional precision analysis59. The range of AUC values greater than 0.7 was considered a fair model performance, while those greater than 0.9 indicated that the model was considered an excellent model performance. Therefore, by considering the AUC values, the excellently performing model was selected to analyze the suitability of bedbugs in Kenya and Africa59,60,61,62.The ASCII format output was then imported into QGIS 3.10.2 (using the QGIS 3.10.2 software, https://qgis.org/downloads/), following its conversion into a raster format file using R software. This was useful for the classification and visualization of the distribution area63,64. The potential suitable distribution of bedbugs was extracted using the Kenyan and African maps. At the same time, Jenks’ natural breaks were also used to reclassify and classify the suitability into five categories, namely: unsuitable (P  More