Philippa Kaur

The haplotype networks, PCoA analysis and STRUCTURE analysis based on the six nuclear loci confirm that S. longifolium is a hybrid between S. emersum and S. gramineum, providing molecular support for previous morphological analyses5. Furthermore, all individuals with intermediate admixture coefficient (Fig. 2b) and private haplotypes only present in one out of six nuclear loci (Fig. 1) suggest that S. longifolium is most likely a F1 hybrid. We thus hypothesized that S. emersum and S. gramineum could likely maintain their species boundary through the post-zygote reproductive isolation mechanism of F1 generation sterility. This hypothesis is possible based on the observations from hybrids in European Russia. The pollen viability was checked in S. longifolium samples from Vysokovskoe Lake and Sabro Lake, and the vast majority of checked pollens were sterile5. In addition, flowering plants of S. longifolium often do not form seeds, or the seeds are puny and significantly inferior to normal seeds in size5. However, the hypothesis is only based on our limited sampling, which is contrary to the conclusion inferred from morphological characteristics that it is fertile and may backcross with parental species1. Further studies with extensive sampling are necessary to test our hypothesis.The chloroplast DNA fragment trnH-psbA was used to infer the direction of hybridization between S. emersum and S. gramineum because chloroplast DNA is maternal inheritance in Sparganium3,4. The hybrid S. longifolium shared haplotypes with S. emersum and S. gramineum simultaneously (Fig. 1). This finding clearly indicates that bidirectional hybridization exists between S. emersum and S. gramineum. At the same time, the different frequency of these two haplotypes in the hybrid (H1, 19.1% vs. H2, 80.9%) means that the direction of hybridization is asymmetric. A variety of factors can lead to asymmetry in natural hybridization, such as flowering time, preference of pollinators, quality and quantity of pollen, cross incompatibility and the abundance of parent species7,8. Rare species usually act as maternal species relative to abundant species9,10. S. gramineum is confined to oligotrophic lakes and its abundance is obviously lower than that of S. emersum1,11. The relatively scarcity combined with the ecology of S. gramineum make it more often act as maternal species when hybridizing with S. emersum.As described by5, the morphological diversification of S. longifolium was also observed in this study. For example, individuals of S. longifolium with emergent and floating-leaved life forms occur concurrently in Zaozer’ye Lake (Supplementary Fig. S2). However, all individuals had the same haplotype H2 as S. gramineum (Fig. 1), suggesting that the direction of hybridization do not determine life form of S. longifolium. In addition, all individuals of S. longifolium sampled here are likely F1 hybrid. Their variable phenotypes could not be associated with traits segregation due to F2 generation or backcross. Detailed ecological investigation combining with research at the genomic level are essential to find out the potential factors leading to morphological diversification of S. longifolium.Here, using sequences of six nuclear loci and one chloroplast DNA fragment, we confirmed that S. longifolium is the hybrid between S. emersum and S. gramineum. The natural hybridization between S. emersum and S. gramineum is bidirectional but the latter mainly acts as maternal species. We also found that all samples of S. longifolium were F1 generations, indicating that S. emersum and S. gramineum could maintain their species boundary through the post-zygote reproductive isolation mechanism of F1 generation sterility. More

Miyata, T., Miyazawa, S. & Yasunaga, T. Two types of amino acid substitutions in protein evolution. J. Mol. Evol. 12, 219–236 (1979).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Li, W.-H., Wu, C.-I. & Luo, C.-C. A new method for estimating synonymous and nonsynonymous rates of nucleotide substitution considering the relative likelihood of nucleotide and codon changes. Mol. Biol. Evol. 2, 150–174 (1985).PubMed 

Google Scholar 
Bielawski, J. P. & Yang, Z. Positive and negative selection in the DAZ gene family. Mol. Biol. Evol. 18, 523–529 (2001).CAS 
PubMed 
Article 

Google Scholar 
Ohta, T. Slightly deleterious mutant substitutions in evolution. Nature 246, 96–98 (1973).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Ohta, T. The nearly neutral theory of molecular evolution. Annu. Rev. Ecol. Evol. Syst. 23, 263–286 (1992).Article 

Google Scholar 
Johnson, K. P. & Seger, J. Elevated rates of nonsynonymous substitution in island birds. Mol. Biol. Evol. 18, 874–881 (2001).CAS 
PubMed 
Article 

Google Scholar 
Woolfit, M. & Bromham, L. Population size and molecular evolution on islands. Proc. Biol. Sci. 272, 2277–2282 (2005).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ross, L., Hardy, N. B., Okusu, A. & Normark, B. B. Large population size predicts the distribution of asexuality in scale insects. Evolution 67, 196–206 (2013).PubMed 
Article 

Google Scholar 
Weber, C. C., Nabholz, B., Romiguier, J. & Ellegren, H. Kr/Kc but not dN/dS correlates positively with body mass in birds, raising implications for inferring lineage-specific selection. Genome Biol. 15, 542 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Brandt, A. et al. Effective purifying selection in ancient asexual oribatid mites. Nat. Commun. 8, 873 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Figuet, E. et al. Life history traits, protein evolution, and the nearly neutral theory in amniotes. Mol. Biol. Evol. 33(6), 1517–1527 (2016).CAS 
PubMed 
Article 

Google Scholar 
Saclier, N. et al. Life history traits impact the nuclear rate of substitution but not the mitochondrial rate in isopods. Mol. Biol. Evol. 35, 2900–2912 (2018).CAS 
PubMed 
Article 

Google Scholar 
Hebert, P. D. The Daphnia of North America: An Illustrated Fauna (on CD-ROM) (CyberNatural Software, Guelph, 1995).
Google Scholar 
Colbourne, J. K. et al. Phylogenetics and evolution of a circumarctic species complex (Cladocera: Daphnia pulex). Biol. J. Linn. Soc. 65, 347–365 (1998).
Google Scholar 
Crease, T. J., Omilian, A. R., Costanzo, K. S. & Taylor, D. J. Transcontinental phylogeography of the Daphnia pulex species complex. PLoS ONE 7, e46620 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Mergeay, J., Verschuren, D. & De Meester, L. Cryptic invasion and dispersal of an American Daphnia in East Africa. Limnol. Oceanogr. 50, 1278–1283 (2005).ADS 
CAS 
Article 

Google Scholar 
Ma, X. et al. Lineage diversity and reproductive modes of the Daphnia pulex group in Chinese lakes and reservoirs. Mol. Phylogenet. Evol. 130, 424–433 (2019).PubMed 
Article 

Google Scholar 
So, M. et al. Invasion and molecular evolution of Daphnia pulex in Japan. Limnol. Oceanogr. 60, 1129–1138 (2015).ADS 
Article 

Google Scholar 
Duggan, I. C. et al. Identifying invertebrate invasions using morphological and molecular analyses: North American Daphnia ‘pulex’ in New Zealand fresh waters. Aquat. Invasions 7, 585–590 (2012).Article 

Google Scholar 
Ye, Z. et al. The rapid, mass invasion of New Zealand by North American Daphnia “pulex”. Limnol. Oceanogr. 66, 2673–2683 (2021).ADS 
Article 

Google Scholar 
Paland, S., Colbourne, J. K. & Lynch, M. Evolutionary history of contagious asexuality in Daphnia pulex. Evolution 59, 800–813 (2005).CAS 
PubMed 
Article 

Google Scholar 
Muller, H. J. The relation of recombination to mutational advance. Mutat. Res. 106, 2–9 (1964).CAS 
PubMed 
Article 

Google Scholar 
Felsenstein, J. The evolutionary advantage of recombination. Genetics 78, 737–756 (1974).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Paland, S. & Lynch, M. Transitions to asexuality result in excess amino acid substitutions. Science 311, 990–992 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Johnson, S. G. & Howard, R. S. Contrasting patterns of synonymous and nonsynonymous sequence evolution in asexual and sexual freshwater snail lineages. Evolution 61, 2728–2735 (2007).PubMed 
Article 

Google Scholar 
Neiman, M. et al. Accelerated mutation accumulation in asexual lineages of a freshwater snail. Mol. Biol. Evol. 27, 954–963 (2010).CAS 
PubMed 
Article 

Google Scholar 
Henry, L., Schwander, T. & Crespi, B. J. Deleterious mutation accumulation in asexual Timema stick insects. Mol. Biol. Evol. 29, 401–408 (2012).CAS 
PubMed 
Article 

Google Scholar 
Tucker, A. E. et al. Population-genomic insights into the evolutionary origin and fate of obligately asexual Daphnia pulex. Proc. Natl. Acad. Sci. 110, 15740–15745 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Colbourne, J. K. et al. The ecoresponsive genome of Daphnia pulex. Science 331, 555–561 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ye, Z. et al. A new reference genome assembly for the microcrustacean Daphnia pulex. G3 (Bethesda) 7, 1405–1416 (2017).CAS 
Article 

Google Scholar 
Keith, N. et al. High mutational rates of large-scale duplication and deletion in Daphnia pulex. Genome Res. 26, 60–69 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hall, D. J. An experimental approach to the dynamics of a natural population of Daphnia galeata mendotae. Ecology 45, 94–112 (1964).Article 

Google Scholar 
McCauley, E., Murdoch, W. W. & Nisbet, R. M. Growth, reproduction, and mortality of Daphnia pulex Leydig: Life at low food. Ecology 4, 505–514 (1990).
Google Scholar 
Xu, S. et al. High mutation rates in the mitochondrial genomes of Daphnia pulex. Mol. Biol. Evol. 29, 763–769 (2012).CAS 
PubMed 
Article 

Google Scholar 
Zheng, Y., Peng, R., Kuro-o, M. & Zeng, X. Exploring patterns and extent of bias in estimating divergence time from mitochondrial DNA sequence data in a particular lineage: A case study of salamanders (Order Caudata). Mol. Biol. Evol. 28, 2521–2535 (2011).CAS 
PubMed 
Article 

Google Scholar 
Zaret, T. M. Predation and Freshwater Communities (Yale University Press, New Haven, 1980).
Google Scholar 
Lynch, M. Predation, competition, and zooplankton community structure: An experimental study. Limnol. Oceanogr. 24, 253–272 (1979).ADS 
Article 

Google Scholar 
Mills, E. L. & Forney, J. L. Impact on Daphnia pulex of predation by young yellow perch in Oneida Lake, New York. Trans. Am. Fish. Soc. 112(2A), 154–161 (1983).Article 

Google Scholar 
Craddock, D. R. Effects of increased water temperature on Daphnia pulex. Fish. Bull. 74, 403–408 (1976).
Google Scholar 
Maruoka, N. & Urabe, J. Inter and intraspecific competitive abilities and the distribution ranges of two Daphnia species in Eurasian continental islands. Popul. Ecol. 62, 353–363 (2020).Article 

Google Scholar 
Dodson, S. I. & Hanazato, T. Commentary on effects of anthropogenic and natural organic chemicals on development, swimming behavior, and reproduction of Daphnia, a key member of aquatic ecosystems. Environ. Health Perspect. 103(Suppl 4), 7–11 (1995).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Claska, M. E. & Gilbert, J. J. The effect of temperature on the response of Daphnia to toxic cyanobacteria. Freshw. Biol. 39, 221–232 (1998).Article 

Google Scholar 
Bast, J. et al. Consequences of asexuality in natural populations: Insights from stick insects. Mol. Biol. Evol. 35, 1668–1677 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hartfield, M. Evolutionary genetic consequences of facultative sex and outcrossing. J Evol Biol 29, 5–22 (2016).CAS 
PubMed 
Article 

Google Scholar 
Hörandl, E. et al. Genome evolution of asexual organisms and the paradox of sex in eukaryotes. In Evolutionary Biology—A Transdisciplinary Approach (ed. Pontarotti, P.) (Springer, Cham, 2020). https://doi.org/10.1007/978-3-030-57246-4_7.Chapter 

Google Scholar 
Lynch, M., Bürger, R., Butcher, D. & Gabriel, W. The mutational meltdown in asexual populations. J. Hered. 84, 339–344 (1993).CAS 
PubMed 
Article 

Google Scholar 
Gordo, I. & Charlesworth, B. The degeneration of asexual haploid populations and the speed of Muller’s ratchet. Genetics 154, 1379–1387 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Downing, J. A. et al. The global abundance and size distribution of lakes, ponds, and impoundments. Limnol. Oceanogr. 51, 2388–2397 (2006).ADS 
Article 

Google Scholar 
McDonald, C. P., Rover, J. A., Stets, E. G. & Striegl, R. G. The regional abundance and size distribution of lakes and reservoirs in the United States and implications for estimates of global lake extent. Limnol. Oceanogr. 57, 597–606 (2012).ADS 
Article 

Google Scholar 
De Meester, L., Góme, A., Okamura, B. & Schwenk, K. The monopolization hypothesis and the dispersal-gene flow paradox in aquatic organisms. Acta Oecol. 23, 121–135 (2002).ADS 
Article 

Google Scholar 
Fukami, T., Bezemer, T. M., Mortimer, S. R. & Van Der Putten, W. H. Species divergence and trait convergence in experimental plant community assembly. Ecol. Lett. 8, 1283–1290 (2005).Article 

Google Scholar 
Makino, T. & Kawata, M. Invasive invertebrates associated with highly duplicated gene content. Mol. Ecol. 28, 1652–1663 (2019).CAS 
PubMed 
Article 

Google Scholar 
Kondrashov, F. A. Gene duplication as a mechanism of genomic adaptation to a changing environment. Proc. R. Soc. Lond. B Biol. Sci. 279, 5048–5057 (2012).
Google Scholar 
Barrick, J. E. & Lenski, R. E. Genome dynamics during experimental evolution. Nat. Rev. Genet. 14, 827–839 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Rocha, E. P. C. Neutral theory, microbial practice: Challenges in bacterial population genetics. Mol. Biol. Evol. 35, 1338–1347 (2018).CAS 
PubMed 
Article 

Google Scholar 
Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol. Biol. Evol. 30, 772–780 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Tanabe, A. S. Kakusan4 and Aminosan: Two programs for comparing nonpartitioned, proportional and separate models for combined molecular phylogenetic analyses of multilocus sequence data. Mol. Ecol. Resour. 11, 914–921 (2011).PubMed 
Article 

Google Scholar 
Stamatakis, A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Tian, X., Ohtsuki, H. & Urabe, J. Evolution of asexual Daphnia pulex in Japan: Variations and covariations of the digestive, morphological and life history traits. BMC Evol. Biol. 19, 122 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Chen, Y. et al. SOAPnuke: A MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience 7, 1–6 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2019).Article 
CAS 

Google Scholar 
Lee, T. H. et al. SNPhylo: A pipeline to construct a phylogenetic tree from huge SNP data. BMC Genomics 15, 162 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
R Core Team, R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, Austria, 2019). https://www.R-project.org/ More

Meropol, S. B., Haupt, A. A. & Debanne, S. M. Incidence and outcomes of infections caused by multidrug-resistant Enterobacteriaceae in Children, 2007–2015. J. Pediatr. Infect. Dis. Soc. 7, 36–45 (2018).Article 

Google Scholar 
Moxon, C. A. & Paulus, S. Beta-lactamases in Enterobacteriaceae infections in children. J. Infect. 72, S41–S49 (2016).PubMed 
Article 

Google Scholar 
Morrissey, I. et al. A review of ten years of the study for monitoring antimicrobial resistance trends (SMART) from 2002 to 2011. Pharmaceuticals 6, 1335–1346 (2013).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Junnila, J. et al. Changing epidemiology of methicillin-resistant Staphylococcus aureus in a low endemicity area—new challenges for MRSA control. Eur. J. Clin. Microbiol. Infect. Dis. 39, 2299–2307 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Milstone, A. M. et al. Methicillin-resistant Staphylococcus aureus colonization and risk of subsequent infection in critically ill children: Importance of preventing nosocomial methicillin-resistant Staphylococcus aureus transmission. Clin. Infect. Dis. 53, 853–859 (2011).PubMed 
PubMed Central 
Article 

Google Scholar 
Lakhundi, S. & Zhang, K. Methicillin-resistant Staphylococcus aureus: Molecular characterization, evolution, and epidemiology. Clin. Microbiol. Rev. 31, e00020-18 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Schlesinger, Y. et al. Methicillin-resistant nasal colonization in children in Jerusalem: Community vs. chronic care institutions. Isr. Med. Assoc. J. 5, 847–851 (2003).PubMed 

Google Scholar 
Liang, B. et al. Active surveillance, drug resistance, and genotypic profiling of Staphylococcus aureus among school-age children in China. Front. Med. 8, 701494 (2021).Article 

Google Scholar 
Del Rosal, T. et al. Staphylococcus aureus nasal colonization in Spanish children. The COSACO Nationwide Surveillance Study. Infect. Drug Resist. 13, 4643–4651 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Feudtner, C., Feinstein, J. A., Zhong, W., Hall, M. & Dai, D. Pediatric complex chronic conditions classification system version 2: Updated for ICD-10 and complex medical technology dependence and transplantation. BMC Pediatr. 14, 199 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Climent Alcalá, F. J., García Fernández de Villalta, M., Escosa García, L., Rodríguez Alonso, A. & Albajara Velasco, L. A. Unidad de niños con patología crónica compleja. Un modelo necesario en nuestros hospitales. Anales de Pediatría 88, 12–18 (2018).PubMed 
Article 

Google Scholar 
Gesualdo, F. et al. Methicillin-resistant Staphylococcus aureus nasal colonization in a department of pediatrics: A cross-sectional study. Ital. J. Pediatr. 40, 3 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Yamamoto, M. et al. Effective surveillance to identify the surgical patients carrying methicillin-resistant Staphylococcus aureus on admission in a pediatric ward. Osaka City Med. J. 62, 1–9 (2016).PubMed 

Google Scholar 
Lukac, P. J., Bonomo, R. A. & Logan, L. K. Extended-spectrum-lactamase-producing Enterobacteriaceae in children: Old foe, emerging threat. Clin. Infect. Dis. https://doi.org/10.1093/cid/civ020 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Fedler, K. A., Biedenbach, D. J. & Jones, R. N. Assessment of pathogen frequency and resistance patterns among pediatric patient isolates: Report from the 2004 SENTRY Antimicrobial Surveillance Program on 3 continents. Diagn. Microbiol. Infect. Dis. 56, 427–436 (2006).CAS 
PubMed 
Article 

Google Scholar 
Caselli, D. et al. Incidence of colonization and bloodstream infection with carbapenem-resistant Enterobacteriaceae in children receiving antineoplastic chemotherapy in Italy. Infect. Dis. 48, 152–155 (2016).Article 

Google Scholar 
Logan, L. K. et al. Multidrug- and Carbapenem-Resistant Pseudomonas aeruginosa in Children, United States, 1999–2012. JPIDSJ piw064 (2016) https://doi.org/10.1093/jpids/piw064.Flokas, M. E., Alevizakos, M., Shehadeh, F., Andreatos, N. & Mylonakis, E. Extended-spectrum β-lactamase-producing Enterobacteriaceae colonisation in long-term care facilities: A systematic review and meta-analysis. Int. J. Antimicrob. Agents 50, 649–656 (2017).CAS 
PubMed 
Article 

Google Scholar 
Bharadwaj, R. et al. Drug-resistant Enterobacteriaceae colonization is associated with healthcare utilization and antimicrobial use among inpatients in Pune, India. BMC Infect. Dis. 18, 504 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Logan, L. K. Carbapenem-resistant Enterobacteriaceae: An emerging problem in children. Clin. Infect. Dis. 55, 852–859 (2012).CAS 
PubMed 
Article 

Google Scholar 
Badal, R. E. et al. Etiology, extended-spectrum β-lactamase rates and antimicrobial susceptibility of gram-negative bacilli causing intra-abdominal infections in patients in general pediatric and pediatric intensive care units—global data from the Study for Monitoring Antimicrobial Resistance Trends 2008 to 2010. Pediatr. Infect. Dis. J. 32, 636–640 (2013).PubMed 
Article 

Google Scholar 
Wang, Q. et al. Risk factors and clinical outcomes for carbapenem-resistant Enterobacteriaceae nosocomial infections. Eur. J. Clin. Microbiol. Infect. Dis. 35, 1679–1689 (2016).CAS 
PubMed 
Article 

Google Scholar 
Sahbudak Bal, Z. et al. The prospective evaluation of risk factors and clinical influence of carbapenem resistance in children with gram-negative bacteria infection. Am. J. Infect. Control 46, 147–153 (2018).CAS 
PubMed 
Article 

Google Scholar 
Simon, T. D. et al. Pediatric medical complexity algorithm: A new method to stratify children by medical complexity. Pediatrics 133, e1647–e1654 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Román, F. et al. Characterization of methicillin-resistant Staphylococcus aureus strains colonizing the nostrils of Spanish children. MicrobiologyOpen 10, e1235 (2021).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
EUCAST. European committee on antimicrobial susceptibility testing breakpoint tables for interpretation of MICs and zone diameters. The European Committee on Antimicrobial Susceptibility Testing. (2018).Oteo, J. et al. Prospective multicenter study of carbapenemase-producing Enterobacteriaceae from 83 hospitals in Spain reveals high in vitro susceptibility to colistin and meropenem. Antimicrob. Agents Chemother. 59, 3406–3412 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Maseda, E. et al. Risk factors for colonization by carbapenemase-producing enterobacteria at admission to a Surgical ICU: A retrospective study. Enferm. Infecc. Microbiol. Clin. 35, 333–337 (2017).PubMed 
Article 

Google Scholar 
Bassetti, M., Nicco, E. & Mikulska, M. Why is community-associated MRSA spreading across the world and how will it change clinical practice?. Int. J. Antimicrob. Agents 34, S15–S19 (2009).CAS 
PubMed 
Article 

Google Scholar 
El Cheikh, M. R., Barbosa, J. M., Caixêta, J. A. S. & Avelino, M. A. G. Microbiology of tracheal secretions: What to expect with children and adolescents with tracheostomies. Int. Arch. Otorhinolaryngol. 22, 50–54 (2018).PubMed 
Article 

Google Scholar 
González-Del Castillo, J. et al. BAHNG score: Predictive model for detection of subjects with the oropharynx colonized by uncommon microorganisms. Rev. Esp Quimioter. 30, 422–428 (2017).PubMed 

Google Scholar 
Hu, X. et al. Risk factors for methicillin-resistant Staphylococcus aureus colonization and infection in patients with human immunodeficiency virus infection: A systematic review and meta-analysis. J. Int. Med. Res. 50, 3000605211063019 (2022).CAS 
PubMed 

Google Scholar 
Gleeson, A., Larkin, P., Walsh, C. & O’Sullivan, N. Methicillin-resistant Staphylococcus aureus: Prevalence, incidence, risk factors, and effects on survival of patients in a specialist palliative care unit: A prospective observational study. Palliat. Med. 30, 374–381 (2016).PubMed 
Article 

Google Scholar 
Hogardt, M. et al. Current prevalence of multidrug-resistant organisms in long-term care facilities in the Rhine-Main district, Germany, 2013. Euro Surveill. 20, 21171 (2015).PubMed 
Article 

Google Scholar 
Warren, D. K. et al. Epidemiology of methicillin-resistant Staphylococcus aureus colonization in a surgical intensive care unit. Infect. Control Hosp. Epidemiol. 27, 1032–1040 (2006).PubMed 
Article 

Google Scholar 
Folgori, L. et al. Healthcare-associated infections in pediatric and neonatal intensive care units: Impact of underlying risk factors and antimicrobial resistance on 30-day case-fatality in Italy and Brazil. Infect. Control Hosp. Epidemiol. 37, 1302–1309 (2016).PubMed 
Article 

Google Scholar 
Béranger, A. et al. Early bacterial infections after pediatric liver transplantation in the era of multidrug-resistant bacteria: Nine-year single-center retrospective experience. Pediatr. Infect. Dis. J. 39, e169–e175 (2020).PubMed 
Article 

Google Scholar 
Bouras, D. et al. Staphylococcus aureus osteoarticular infections in children: An 8-year review of molecular microbiology, antibiotic resistance and clinical characteristics. J. Med. Microbiol. 67, 1753–1760 (2018).MathSciNet 
CAS 
PubMed 
Article 

Google Scholar 
Rodriguez, M., Hogan, P. G., Krauss, M., Warren, D. K. & Fritz, S. A. Measurement and impact of Staphylococcus aureus colonization pressure in households. J. Pediatr. Infect. Dis. Soc. 2, 147–154 (2013).Article 

Google Scholar 
Messina, N. L., Williamson, D. A., Robins-Browne, R., Bryant, P. A. & Curtis, N. Risk factors for carriage of antibiotic-resistant bacteria in healthy children in the community: A systematic review. Pediatr. Infect. Dis. J. 39, 397–405 (2020).PubMed 
Article 

Google Scholar 
Dualleh, N. et al. Colonization with multiresistant bacteria in acute hospital care: The association of prior antibiotic consumption as a risk factor. J. Antimicrob. Chemother. 75, 3675–3681 (2020).CAS 
PubMed 
Article 

Google Scholar 
Daskalaki, M. et al. Panton-Valentine leukocidin-positive Staphylococcus aureus skin and soft tissue infections among children in an emergency department in Madrid, Spain. Clin. Microbiol. Infect. 16, 74–77 (2010).CAS 
PubMed 
Article 

Google Scholar 
Aguilera-Alonso, D., Escosa-García, L., Saavedra-Lozano, J., Cercenado, E. & Baquero-Artigao, F. Carbapenem-resistant gram-negative bacterial infections in children. Antimicrob. Agents Chemother. 64, e02183-e2219 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Phichaphop, C. et al. High prevalence of multidrug-resistant gram-negative bacterial infection following pediatric liver transplantation. Medicine 99, e23169 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Tacconelli, E. et al. ESCMID guidelines for the management of the infection control measures to reduce transmission of multidrug-resistant Gram-negative bacteria in hospitalized patients. Clin. Microbiol. Infect. 20, 1–55 (2014).PubMed 
Article 

Google Scholar 
McConville, T. H., Sullivan, S. B., Gomez-Simmonds, A., Whittier, S. & Uhlemann, A.-C. Carbapenem-resistant Enterobacteriaceae colonization (CRE) and subsequent risk of infection and 90-day mortality in critically ill patients, an observational study. PLoS ONE 12, e0186195 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Tamma, P. D. et al. The likelihood of developing a carbapenem-resistant Enterobacteriaceae Infection during a hospital stay. Antimicrob. Agents Chemother. 63, e00757-e819 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Detsis, M., Karanika, S. & Mylonakis, E. ICU acquisition rate, risk factors, and clinical significance of digestive tract colonization with extended-spectrum beta-lactamase-producing Enterobacteriaceae: A systematic review and meta-analysis. Crit. Care Med. 45, 705–714 (2017).PubMed 
Article 

Google Scholar  More

Fossil dataWe downloaded all fossil occurrences recorded for the order Scleractinia at the species level from the Paleobiology Database (PBDB – paleobiodb.org; accessed on 3 August 2021). This is the most comprehensive repository for palaeontological data in reef corals to date. Due to the nature of the data, no ethics approval was required. To minimize identification issues, we excluded taxa with uncertain generic and species assignments (i.e., classified as aff. and cf.) and only selected species that had accepted names. We also selected the variables classification and palaeoenvironment from the output options to facilitate taxonomic and environmental filters applied in downstream analyses. The full dataset consisted of 24,011 occurrences across 4235 species, spanning over 250 Myr of coral evolution from the Triassic to the present. Although our focus here lies on the Cenozoic, we used the complete fossil dataset (i.e., including all of the occurrences) to have estimates of the diversification dynamics in scleractinian corals throughout the whole timespan of their evolution.Evolutionary ratesWith the full palaeontological dataset, we estimated evolutionary rates through time in scleractinian corals using the Bayesian framework of the program PyRate (v3.0)12,36,37. This program uses fossil occurrence data to calculate the temporal variation in rates of preservation, speciation and extinction, while incorporating multiple sources of uncertainty12. At its core implementation, PyRate jointly estimates the times of origination (Ts) and extinction (Te) for each fossil lineage; the fossilization and sampling parameters that determine preservation rates (q); and the overall rates of speciation (λ) and extinction (μ) through time36. Recently, the program has been upgraded to include a reversible jump Markov Chain Monte Carlo (rjMCMC) algorithm to estimate diversification rate heterogeneity, which provides more accurate and precise estimates than other commonly used methods12. Therefore, despite the inherent bias of the fossil record (i.e., estimates are conditioned on sampled lineages), PyRate is a robust method to quantify speciation and extinction rates, and their respective temporal shifts, from fossil occurrence data.Extant taxa can also be included in the PyRate framework as long as they are also represented in the fossil record. This is done to extend the fossil geologic ranges to the recent times. Hence, the first step in our analysis was to identify which species in our dataset is still alive at the present. To do this, we matched the accepted species names in the PBDB dataset with those from the extant species dataset of Huang et al.38. Subsequently, we split our dataset into eleven independent subsets, with the goal of keeping each subset with an equal number of species. Each data subset included a random selection of species with their respective occurrences, which was enough to calculate Ts and Te (see below). This was done to avoid convergence issues, given the large size of our dataset and the consequent complexity of the model37. For each of our subsets, we generated fifty replicates by resampling the fossil ages from their temporal ranges to account for the uncertainty associated with the age of occurrences. We then used the maximum-likelihood test in PyRate to compare between three models of fossil preservation12: the homogeneous Poisson process (HPP; q is constant through time); the nonhomogeneous Poisson process (NHPP; q varies throughout the lifespan of a species); and the time-variable Poisson process (TPP; q varies across geological epochs). The latter model (TPP) was selected across all of our data subsets (Supplementary Table 2).After selecting the preservation model, we first focused on assessing the estimates of times of origination and extinction in each data subset, rather than using the full dataset to jointly estimate all parameters at once as in the original implementation of PyRate. This further reduced the complexity of the model and allowed for more precise parameter estimates. For each replicate in all of our data subsets, we approximated the posterior distribution of Ts and Te through a 50 million generation run of the rjMCMC algorithm under the TPP, sampling parameters every 40 thousand iterations. At the end of each run, we discarded 20% of the samples as burn-in and assessed chain convergence through the effective sample sizes of posterior parameter estimates, using the software Tracer39 (v1.7.1).From the results of this first set of models, we extracted the median estimates of Ts and Te across replicates, and we merged the estimates from the eleven independent data subsets. This merged data frame contained estimated times of origination and extinction for all coral lineages within our fossil dataset. We then used this merged Ts and Te data frame as input for another rjMCMC chain to finally estimate overall λ and μ through time, by applying the option -d in PyRate. In this option, Ts and Te for all fossil lineages are given as fixed values and, therefore, are not estimated by the model. The chain for this model was run for 100 million generations, sampling parameters at every 40 thousand iterations. Once again, we excluded 20% of the initial samples as burn-in and checked model convergence using Tracer. Finally, we calculated net diversification rates through time by subtracting the post burn-in samples of μ from λ.To explore the taxonomic idiosyncrasies in the evolutionary rates of reef corals, we selected the most abundant families on present-day coral reefs in terms of the number of colonies per area18 (Acroporidae, Agariciidae, Merulinidae, Mussidae, Pocilloporidae, and Poritidae). Altogether, species within these families account for ~40% of the total extant diversity in Scleractinia. These families also account for most of the occurrences in the PBDB fossil dataset (excluding extinct families, which are generally older and had little temporal overlap with extant ones): Acroporidae (1457 occ. in 165 spp.); Agariciidae (722 occ. in 89 spp.); Merulinidae (2464 occ. in 229 spp.); Mussidae (1146 occ. in 100 spp.); Pocilloporidae (615 occ. in 64 spp.); and Poritidae (1149 occ. in 91 spp.). Therefore, from our full dataset, we selected six independent ones encompassing all species in each of the selected families. We also selected only species that are classified as reef-associated within these families, since we were specifically interested in these environments. This selection had a negligible effect on the size of the individual datasets, given that the vast majority of fossil species within these families are reef-associated. In each family, we followed the same modelling steps described above to estimate μ and λ, and diversity trajectories. However, this time it was not necessary to split the datasets into subsets, given that each family has far less occurrences than the full dataset. We started by comparing models of preservation, which showed the TPP as the best supported for all families (Supplementary Table 3). Then we created fifty replicates by resampling fossil ages to accommodate the uncertainty associated with the time of occurrences. For each replicate, we ran the rjMCMC algorithm for 50 million generations under the TPP model, with a sampling frequency of 40 thousand iterations. We discarded initial 20% of the samples as burn-in, and assessed convergence through Tracer. We then combined all replicates, resampling 100 random samples from each replicate to assess the estimates of μ and λ through time for each family. Finally, we extracted diversity trajectories in each family for all of the replicates by applying the -ltt option in PyRate, which generates a table with estimated range-through diversity at every 0.1 Myr. From these trajectories, we calculated the mean difference in diversity (slope in species per 0.1 Myr) between subsequent time samples backwards from the present (i.e., diversity in time t was subtracted from diversity in time t-1) using the diff function in R (v4.0.3).As an alternative to PyRate, we also calculated the diversity dynamics of reef coral fossils using the R package divDyn40, which combines a range of published methods for quantifying fossil diversification rates. Differently from PyRate, the metrics applied in divDyn require that the fossil occurrences are split into discrete time bins. Therefore, these metrics treat the origination and extinction rates as independent parameters in each bin, while PyRate is designed to detect rate heterogeneity through a continuous time setting12. Our goal here, however, was not to compare models but to assess the robustness of our rate patterns and diversity trajectories using alternative methods. We divided our dataset into one-million-year time bins to have enough temporal resolution for rate calculations. To account for the uncertainty in the assignment of fossil ages, we created 50 binned replicates by sampling the age of each occurrence from a random uniform distribution, with bounds defined by the age ranges provided in the PBDB dataset. We then used the divDyn function to calculate the per capita rates of origination and extinction through time (based on the rate equations by Foote41) for all scleractinians (Supplementary Fig. 5a) and for reef-associated acroporids alone (Supplementary Fig. 5b). We also used the same procedure to generate range-through diversity curves for each of the six families selected previously, to compare with the curves generated by PyRate (Fig. 2a). Although the rate results differed between the PyRate (Fig. 1) and the divDyn (Supplementary Fig. 5) approaches, the general patterns remained unchanged. Rates are more volatile through time in divDyn estimates, with larger confidence intervals, which is expected from the metrics applied in the package12,42. Yet, we found the same peaks in extinction for Scleractinia: at the Cretaceous-Paleogene and Eocene-Oligocene boundaries, and at the Pliocene-Pleistocene (Supplementary Fig. 5a). The recent peak in speciation in Acroporidae was also detected, although less strong (Supplementary Fig. 5b). Despite these slight differences in rate estimates, the diversity curves reconstructed through divDyn (Supplementary Fig. 6) mirrored almost exactly the ones found with PyRate (Fig. 2), demonstrating that the overall macroevolutionary trends described herein (Figs. 1 and 2) are robust to methodological choices.Diversity-dependent modelsTo assess the effects of diversity dependency on the evolution of reef coral lineages, we implemented the Multivariate Birth-Death model (MBD)11 within the PyRate framework. This method was first described as the Multiple Clade Diversity Dependence model (MCDD)19, in which rates of speciation and extinction are modelled as having linear correlations with the diversity trajectories of other clades. At its original implementation, the MCDD was developed to assess the effects of negative interactions, where increasing species diversity in one group can suppress speciation rates and/or promote extinction in itself or in other ecologically similar clades19. However, the model also incorporates the possibility of positive interactions, where increasing diversity in one clade can correlate with enhanced rates of speciation or buffered extinction. Through further model developments43, the MCDD was updated to also include a horseshoe prior44 on the diversity-dependence parameters, which helped controlling for overparameterization and enhanced the power of the model to recover true effects43. More recently, this model took its current form as the MBD11, with the additional possibilities of including environmental correlates and setting exponential, rather than just linear, correlations.We first applied the MBD to estimate the diversity-dependent effects of individual extant coral families (i.e., the ones selected in the previous analysis; see Evolutionary rates) in their combined diversity trajectories. From the rjMCMC model results for individual families, we extracted estimates of Ts and Te in each of the fifty replicates and merged them across families. This merged dataset with fifty replicates of Ts and Te was then used as input for the MBD model, where we set the relative diversity trajectories of each individual family as predictors. We also included three key environmental predictors—paleotemperature, sea level and rate of sea-level change—to assess their influence in overall evolutionary rates. The paleotemperature data was obtained from Westerhold et al.45, and consists of global mean temperature estimates for the last 66 million years, averaged across 0.1 Myr time bins. Eustatic sea-level data was downloaded from Miller et al.46, and contains estimates of sea level for the last 100 million years in comparison to present-day levels, also split in ~0.1 Myr time bins. With this dataset, we calculated the average rate of sea-level change per million years, as measured from the absolute difference between subsequent sea-level values backwards in time (i.e., sea level in time t was subtracted from sea level in time t-1). These environmental factors were rescaled between 0 and 1 to maintain all predictors on the same relative scale.Under our MBD model, the speciation and extinction rates of all families combined could change through time and through correlations with the relative diversity of individual families or environmental factors. The strength and directionality (positive or negative) of the correlations are also jointly estimated for each predictor within the model11. We ran both linear and exponential correlation models (see formulas in Lehtonen et al.11) in each of our fifty replicates for 25 million generations, sampling parameters at every 25 thousand iterations. We then compared the linear and exponential models through the posterior harmonic means of their log likelihoods, which supported the exponential one as having a better fit. From the posterior estimates, we summarized the speciation (Fig. 3c) and extinction (Fig. 3d) correlation parameters (i.e., the strength of the effect) by calculating their median and 95% Highest Posterior Density (HPD) interval across replicates. Finally, we also summarized the effect of families on lineage turnover (Fig. 3e), which we conceptualize as the sum of the effects on speciation and extinction.The MBD model also provides posterior samples of the weight of the correlation parameters, which is estimated through the horseshoe prior11. In essence, this prior is able to reliably distinguish correlation parameters that should be considered noise from those that represent a true signal in the data11. The parameterization of the horseshoe prior contains local and global Bayesian shrinkage parameters44 from which shrinkage weights (w) can be calculated (see formulas in Lehtonen et al.11). These shrinkage weights associated with each correlation parameter in the MBD model vary between 0 and 1, with values closer to 0 representing noise and values closer to 1 representing a true signal. Through simulations, it has been shown that values of w  > 0.5 indicate that the correlation parameter in question significantly differs from the background noise, being the correlation positive or negative43. However, as a conservative way to infer the weight of correlation parameters, here we use a value of w  > 0.7 to detect significance. This value was calculated for each diversity-dependence parameter (speciation, extinction and turnover) from the median values drawn from the model posteriors.The spatial distribution of reef-associated taxa varied considerably throughout the Cenozoic, with biodiversity hotspots moving halfway across the globe47. Therefore, the best way to capture this dynamic biogeographic history in reef corals is by analysing global diversity patterns like we did in our main MBD model. However, to assess the robustness of our diversity-dependent results against the influence of geographic scale and site co-occurrences, we repeated all the modelling steps described above with two data subsets. First, we selected only fossil species that have occurrences in the Indo-Pacific Ocean (i.e., 30°W–180°W) within the six families. Second, we excluded sites in which the Acroporidae did not co-occur with the other families. In each of these data subsets, we calculated diversity trajectories and used them as predictors in a separate MBD model. These models had a merged dataset of Ts and Te of all species included in each case (Indo-Pacific and co-occurrences) as a response variable.Finally, we followed the same modelling procedures described above to investigate the diversity-dependent effects in family pairwise analyses. We applied the MBD model to assess the effects of all other families in each individual family at a time, while also estimating correlations with the key environmental predictors. From the rjMCMC model results for individual families, we extracted the fifty replicates of estimated Ts and Te. Each replicate was then used as input for an MBD run using the relative diversity trajectories of each other individual family as predictors, along with the environmental variables. Once again, we ran 25 million generations of the MBD, with a sampling frequency of 25 thousand, using both linear and exponential correlation models in each age replicate. For all families, we found that the exponential model had a better fit. We then summarized the correlation parameters and the shrinkage weights (Supplementary Fig. 7) derived from the exponential models per family by calculating the median and 95% confidence intervals across replicates.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Ollerton, J. Pollinators & Pollination: Nature and Society (Pelagic Publishing, 2021).
Google Scholar 
Delaplane, K. S. Crop pollination By Bees: Evolution, Ecology, Conservation, and Management (CABI, 2021).
Google Scholar 
Aizen, M. A., Garibaldi, L. A., Cunningham, S. A. & Klein, A. M. Long-term global trends in crop yield and production reveal no current pollination shortage but increasing pollinator dependency. Curr. Biol. 18, 1572–1575 (2008).CAS 
PubMed 

Google Scholar 
Koh, I. et al. Modeling the status, trends, and impacts of wild bee abundance in the United States. Proc. Natl. Acad. Sci. 113, 140–145 (2016).ADS 
CAS 
PubMed 

Google Scholar 
Jordan, A., Patch, H. M., Grozinger, C. M. & Khanna, V. Economic dependence and vulnerability of United States agricultural sector on insect-mediated pollination service. Environ. Sci. Technol. 55, 2243–2253 (2021).ADS 
CAS 
PubMed 

Google Scholar 
Biddinger, D. J. & Rajotte, E. G. Integrated pest and pollinator management: Adding a new dimension to an accepted paradigm. Curr. Opin. Insect Sci. 10, 204–209 (2015).PubMed 

Google Scholar 
Egan, P. A., Dicks, L. V., Hokkanen, H. M. T. & Stenberg, J. A. Delivering integrated pest and pollinator management (IPPM). Trends Plant Sci. 25, 577–589 (2020).CAS 
PubMed 

Google Scholar 
Flöhr, A., Stenberg, J. A. & Egan, P. A. The Joint Economic Impact Level (jEIL): A Decision Metric for Integrated Pest and Pollinator Management. In Integrative Biological Control 17–38 (Springer, 2020).
Google Scholar 
Krupke, C. H., Hunt, G. J., Eitzer, B. D., Andino, G. & Given, K. Multiple routes of pesticide exposure for honey bees living near agricultural fields. PLoS ONE 7, e29268 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Traynor, K. S. et al. In-hive pesticide exposome: Assessing risks to migratory honey bees from in-hive pesticide contamination in the Eastern United States. Sci. Rep. 6, 1–16 (2016).
Google Scholar 
Mullin, C. A. et al. High levels of miticides and agrochemicals in North American apiaries: Implications for honey bee health. PLoS ONE 5, e9754 (2010).ADS 
PubMed 
PubMed Central 

Google Scholar 
Ravoet, J., Reybroeck, W. & de Graaf, D. C. Pesticides for apicultural and/or agricultural application found in belgian honey bee wax combs. Bull. Environ. Contam. Toxicol. 94, 543–548 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Liang, H. C., Bilon, N. & Hay, M. T. Analytical methods for pesticide residues. Water Environ. Res. 86, 2132–2155 (2014).CAS 

Google Scholar 
Calatayud-Vernich, P., Calatayud, F., Simó, E. & Picó, Y. Efficiency of QuEChERS approach for determining 52 pesticide residues in honey and honey bees. MethodsX 3, 452–458 (2016).PubMed 
PubMed Central 

Google Scholar 
Fernandez, M., Pico, Y. & Manes, J. Analytical methods for pesticide residue determination in bee products. J. Food Prot. 65, 1502–1511 (2002).CAS 
PubMed 

Google Scholar 
Vázquez, P. P., Lozano, A., Uclés, S., Ramos, M. M. G. & Fernández-Alba, A. R. A sensitive and efficient method for routine pesticide multiresidue analysis in bee pollen samples using gas and liquid chromatography coupled to tandem mass spectrometry. J. Chromatogr. A 1426, 161–173 (2015).PubMed 

Google Scholar 
Stoner, K. A., Cowles, R. S., Nurse, A. & Eitzer, B. D. Tracking pesticide residues to a plant genus using palynology in pollen trapped from honey bees (Hymenoptera: Apidae) at ornamental plant nurseries. Environ. Entomol. 48, 351–362 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Colwell, M. J., Williams, G. R., Evans, R. C. & Shutler, D. Honey bee-collected pollen in agro-ecosystems reveals diet diversity, diet quality, and pesticide exposure. Ecol. Evol. 7, 7243–7253 (2017).PubMed 
PubMed Central 

Google Scholar 
Favaro, R. et al. Botanical origin of pesticide residues in pollen loads collected by honeybees during and after apple bloom. Front. Physiol. 10, 1069 (2019).PubMed 
PubMed Central 

Google Scholar 
Tosi, S., Costa, C., Vesco, U., Quaglia, G. & Guido, G. A 3-year survey of Italian honey bee-collected pollen reveals widespread contamination by agricultural pesticides. Sci. Total Environ. 615, 208–218 (2018).ADS 
CAS 
PubMed 

Google Scholar 
Chaimanee, V., Chantawannakul, P., Khongphinitbunjong, K., Kamyo, T. & Pettis, J. S. Comparative pesticide exposure to Apis mellifera via honey bee-collected pollen in agricultural and non-agricultural areas of Northern Thailand. J. Apic. Res. 58, 720–729 (2019).
Google Scholar 
Friedle, C., Wallner, K., Rosenkranz, P., Martens, D. & Vetter, W. Pesticide residues in daily bee pollen samples (April–July) from an intensive agricultural region in Southern Germany. Environ. Sci. Pollut. Res. https://doi.org/10.1007/s11356-020-12318-2 (2021).Article 

Google Scholar 
Urbanowicz, C. et al. Low maize pollen collection and low pesticide risk to honey bees in heterogeneous agricultural landscapes. Apidologie 50, 379–390 (2019).
Google Scholar 
Stoner, K. A. & Eitzer, B. D. Using a hazard quotient to evaluate pesticide residues detected in pollen trapped from honey bees (Apis mellifera) in Connecticut. PLoS ONE 8, e77550 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
McArt, S. H., Fersch, A. A., Milano, N. J., Truitt, L. L. & Böröczky, K. High pesticide risk to honey bees despite low focal crop pollen collection during pollination of a mass blooming crop. Sci. Rep. 7, 1–10 (2017).
Google Scholar 
Calatayud-Vernich, P., Calatayud, F., Simó, E., Pascual Aguilar, J. A. & Picó, Y. A two-year monitoring of pesticide hazard in-hive: High honey bee mortality rates during insecticide poisoning episodes in apiaries located near agricultural settings. Chemosphere 232, 471–480 (2019).ADS 
CAS 
PubMed 

Google Scholar 
García-Valcárcel, A. I., Martínez-Ferrer, M. T., Campos-Rivela, J. M. & Hernando Guil, M. D. Analysis of pesticide residues in honeybee (Apis mellifera L.) and in corbicular pollen: Exposure in citrus orchard with an integrated pest management system. Talanta 204, 153–162 (2019).PubMed 

Google Scholar 
Fulton, C. A. et al. An assessment of pesticide exposures and land use of honey bees in Virginia. Chemosphere 222, 489–493 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Calatayud-Vernich, P., Calatayud, F., Simó, E. & Picó, Y. Pesticide residues in honey bees, pollen and beeswax: Assessing beehive exposure. Environ. Pollut. 241, 106–114 (2018).CAS 
PubMed 

Google Scholar 
Ostiguy, N. et al. Honey bee exposure to pesticides: A four-year nationwide study. Insects 10, 13 (2019).PubMed Central 

Google Scholar 
Seeley, T. D. The honey bee colony as a superorganism. Am. Sci. 77, 546–553 (1989).ADS 

Google Scholar 
Thompson, H. M. & Maus, C. The relevance of sublethal effects in honey bee testing for pesticide risk assessment. Pest Manag. Sci. 63, 1058–1061 (2007).CAS 
PubMed 

Google Scholar 
Sponsler, D. B. & Johnson, R. M. Mechanistic modeling of pesticide exposure: The missing keystone of honey bee toxicology. Environ. Toxicol. Chem. 36, 871–881 (2017).CAS 
PubMed 

Google Scholar 
Gradish, A. E. et al. Comparison of pesticide exposure in honey bees (Hymenoptera: Apidae) and Bumble Bees (Hymenoptera: Apidae): implications for risk assessments. Environ. Entomol. 48, 12–21 (2019).PubMed 

Google Scholar 
Tosi, S. & Nieh, J. C. Lethal and sublethal synergistic effects of a new systemic pesticide, flupyradifurone (Sivanto®), on honeybees. Proc. R. Soc. B 286, 20190433 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Iverson, A., Hale, C., Richardson, L., Miller, O. & McArt, S. Synergistic effects of three sterol biosynthesis inhibiting fungicides on the toxicity of a pyrethroid and neonicotinoid insecticide to bumble bees. Apidologie 50, 733 (2019).CAS 

Google Scholar 
Siviter, H. et al. Agrochemicals interact synergistically to increase bee mortality. Nature 596, 389–392 (2021).ADS 
CAS 
PubMed 

Google Scholar 
Graham, K. K. et al. Identities, concentrations, and sources of pesticide exposure in pollen collected by managed bees during blueberry pollination. Sci. Rep. 11, 1–13 (2021).
Google Scholar 
EFSA. Guidance on the risk assessment of plant protection products on bees (Apis mellifera, Bombus spp. and solitary bees). EFSA J. https://doi.org/10.2903/j.efsa.2013.3295 (2013).Article 

Google Scholar 
EPA. Guidance for assessing pesticide risks to bees. (2014).USDA APHIS. Wax Sampling Protocol for the National Honey Bee Disease Survey. (2018).European Committee for Standardization. Foods of plant origin – Multimethod for the determination of pesticide residues using GC- and LC-based analysis following acetonitrile extraction/partitioning and clean-up by dispersive SPE – Modular QuEChERS-method. (2018).Couvillon, M. J. et al. Honey bee foraging distance depends on month and forage type. Apidologie 46, 61–70 (2015).
Google Scholar 
Knight, M. E. et al. An interspecific comparison of foraging range and nest density of four bumblebee (Bombus) species. Mol. Ecol. 14, 1811–1820 (2005).CAS 
PubMed 

Google Scholar 
McArt, S. H., Urbanowicz, C., Mccoshum, S., Irwin, R. E. & Adler, L. S. Landscape predictors of pathogen prevalence and range contractions in US bumblebees. Proc. R. Soc. B 284, 20172181 (2017).PubMed 
PubMed Central 

Google Scholar 
USDA NASS. USDA National Agricultural Statistics Service Cropland Data Layer. (2018).R Core Team. R: A Language and Environment for Statistical Computing. (2019).GraphPad Software. GraphPad Prism. (2017).Hartig, F. DHARMa: Residual Diagnostics for Hierarchical (Multi-Level/Mixed) Regression Models. (2019).Barton, K. MuMIn: Multi-Model Inference. (2019).Fox, J. & Weisburg, S. An {R} Companion to Applied Regression. (2011).Lefcheck, J. S. piecewiseSEM: Piecewise structural equation modeling in R for ecology, evolution, and systematics. Methods Ecol. Evol. 7, 573–579 (2015).
Google Scholar 
Hothorn, T., Bretz, F. & Westfall, P. Simultaneous inference in general parametric models. Biom. J. 50, 346–363 (2008).MathSciNet 
PubMed 
MATH 

Google Scholar 
Fox, J. RcmdrMisc: R Commander Miscellaneous Functions. (2020).Bhattacharya, M., Primack, R. B. & Gerwein, J. Are roads and railroads barriers to bumblebee movement in a temperate suburban conservation area?. Biol. Conserv. 109, 37–45 (2003).
Google Scholar 
Fragoso, F. P. & Brunet, J. Patch fidelity of honey bees and bumble bees differs and is affected by spatial configuration. In Entomological Society of America Annual Meeting, Plant-Insect Ecosystems (2021).Javorek, S. K., Mackenzie, K. E. & Vander Kloet, S. P. Comparative (Hymenoptera: Apoidea) on Lowbush Blueberry (Ericaceae: Vaccinium angustifolium). Ann. Entomol. Soc. Am. 95, 345–351 (2002).
Google Scholar 
Sandrock, C. et al. Impact of chronic neonicotinoid exposure on honeybee colony performance and queen supersedure. PLoS ONE 9, e103592 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Wood, S. C., Kozii, I. V., Koziy, R. V., Epp, T. & Simko, E. Comparative chronic toxicity of three neonicotinoids on New Zealand packaged honey bees. PLoS ONE 13, e0190517 (2018).PubMed 
PubMed Central 

Google Scholar 
Degrandi-Hoffman, G., Graham, H., Ahumada, F., Smart, M. & Ziolkowski, N. The economics of honey bee (Hymenoptera: Apidae) management and overwintering strategies for colonies used to pollinate almonds. J. Econ. Entomol. https://doi.org/10.1093/jee/toz213 (2019).Article 
PubMed 

Google Scholar 
Biddinger, D. J. et al. Comparative toxicities and synergism of apple orchard pesticides to Apis mellifera (L.) and Osmia cornifrons (Radoszkowski). PLoS ONE 8, e72587 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Johnson, R. M., Dahlgren, L., Siegfried, B. D. & Ellis, M. D. Acaricide, fungicide and drug Interactions in Honey Bees (Apis mellifera). PLoS ONE 8, e54092 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Jaffe, B. D., Lois, A. N. & Guédot, C. Effect of fungicide on pollen foraging by honeybees (Hymenoptera: Apidae) in cranberry differs by fungicide type. J. Econ. Entomol. 112, 499–503 (2019).CAS 
PubMed 

Google Scholar 
Sanchez-Bayo, F. & Goka, K. Pesticide residues and bees: A risk assessment. PLoS ONE 9, e94482 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Wang, Y., Zhu, Y. C. & Li, W. Interaction patterns and combined toxic effects of acetamiprid in combination with seven pesticides on honey bee (Apis mellifera L.). Ecotoxicol. Environ. Saf. 190, 110100 (2020).CAS 
PubMed 

Google Scholar 
Wang, Y., Zhu, Y. C. & Li, W. Comparative examination on synergistic toxicities of chlorpyrifos, acephate, or tetraconazole mixed with pyrethroid insecticides to honey bees (Apis mellifera L.). Environ. Sci. Pollut. Res. 27, 6971–6980 (2019).
Google Scholar 
Sgolastra, F. et al. Synergistic mortality between a neonicotinoid insecticide and an ergosterol-biosynthesis-inhibiting fungicide in three bee species. Pest Manag. Sci. 73, 1236–1243 (2017).CAS 
PubMed 

Google Scholar 
Azpiazu, C. et al. Chronic oral exposure to field-realistic pesticide combinations via pollen and nectar: effects on feeding and thermal performance in a solitary bee. Sci. Reports 9, 1–11 (2019).CAS 

Google Scholar 
Becher, M. A., Hildenbrandt, H., Hemelrijk, C. K. & Moritz, R. F. A. Brood temperature, task division and colony survival in honeybees: A model. Ecol. Modell. 221, 769–776 (2010).
Google Scholar 
Zhu, W., Schmehl, D. R., Mullin, C. A. & Frazier, J. L. Four common pesticides, their mixtures and a formulation solvent in the hive environment have high oral toxicity to honey bee larvae. PLoS ONE 9, e77547 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Dai, P. et al. Chronic toxicity of clothianidin, imidacloprid, chlorpyrifos, and dimethoate to Apis mellifera L. larvae reared in vitro. Pest Manag. Sci. 75, 29–36 (2019).CAS 
PubMed 

Google Scholar 
Campbell, J. B. et al. The fungicide Pristine® inhibits mitochondrial function in vitro but not flight metabolic rates in honey bees. J. Insect Physiol. 86, 11–16 (2016).CAS 
PubMed 

Google Scholar 
DesJardins, N. S. et al. A common fungicide, Pristine®, impairs olfactory associative learning performance in honey bees (Apis mellifera). Environ. Pollut. 288, 117720 (2021).CAS 
PubMed 

Google Scholar 
Fisher, A. et al. Colony field test reveals dramatically higher toxicity of a widely-used mito-toxic fungicide on honey bees (Apis mellifera). Environ. Pollut. 269, 115964 (2021).CAS 
PubMed 

Google Scholar 
Mahefarisoa, K. L., Simon Delso, N., Zaninotto, V., Colin, M. E. & Bonmatin, J. M. The threat of veterinary medicinal products and biocides on pollinators: A one health perspective. One Heal. 12, 100237 (2021).CAS 

Google Scholar 
Christen, V., Schirrmann, M., Frey, J. E. & Fent, K. Global transcriptomic effects of environmentally relevant concentrations of the neonicotinoids clothianidin, imidacloprid, and thiamethoxam in the brain of honey bees (Apis mellifera). Environ. Sci. Technol. 52, 7534–7544 (2018).ADS 
CAS 
PubMed 

Google Scholar 
Tison, L., Duer, A., Púčiková, V., Greggers, U. & Menzel, R. Detrimental effects of clothianidin on foraging and dance communication in honey bees. PLoS ONE 15, e0241134 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Tison, L., Rößner, A., Gerschewski, S. & Menzel, R. The neonicotinoid clothianidin impairs memory processing in honey bees. Ecotoxicol. Environ. Saf. 180, 139–145 (2019).CAS 
PubMed 

Google Scholar 
Morfin, N., Goodwin, P. H., Correa-Benitez, A. & Guzman-Novoa, E. Sublethal exposure to clothianidin during the larval stage causes long-term impairment of hygienic and foraging behaviours of honey bees. Apidologie 50, 595–605 (2019).CAS 

Google Scholar 
Yao, J., Zhu, Y. C. & Adamczyk, J. Responses of honey bees to lethal and sublethal doses of formulated clothianidin alone and mixtures. J. Econ. Entomol. 111, 1517–1525 (2018).CAS 
PubMed 

Google Scholar 
Bortolotti, L. et al. Effects of sub-lethal imidacloprid doses on the homing rate and foraging activity of honey bees. Bull. Insectology 56, 63–67 (2003).
Google Scholar 
Yang, E. C., Chuang, Y. C., Chen, Y. L. & Chang, L. H. Abnormal foraging behavior induced by sublethal dosage of imidacloprid in the honey bee (Hymenoptera: Apidae). J. Econ. Entomol. 101, 1743–1748 (2008).CAS 
PubMed 

Google Scholar 
Karahan, A., Cakmak, I., Hranitz, J. M., Karaca, I. & Wells, H. Sublethal imidacloprid effects on honey bee flower choices when foraging. Ecotoxicology 24, 2017–2025 (2015).CAS 
PubMed 

Google Scholar 
Dively, G. P., Embrey, M. S., Kamel, A., Hawthorne, D. J. & Pettis, J. S. Assessment of chronic sublethal effects of imidacloprid on honey bee colony health. PLoS ONE 10, 1–25 (2015).
Google Scholar 
Meikle, W. G. et al. Sublethal effects of imidacloprid on honey bee colony growth and activity at three sites in the U.S.. PLoS ONE https://doi.org/10.1371/journal.pone.0168603 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Wu, Y.-Y. et al. Sublethal effects of imidacloprid on targeting muscle and ribosomal protein related genes in the honey bee Apis mellifera L.. Sci. Rep. https://doi.org/10.1038/s41598-017-16245-0 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Wu, M.-C., Chang, Y.-W., Lu, K.-H. & Yang, E.-C. Gene expression changes in honey bees induced by sublethal imidacloprid exposure during the larval stage. Insect Biochem. Mol. Biol. 88, 12–20 (2017).CAS 
PubMed 

Google Scholar 
Peng, Y.-C. & Yang, E.-C. Sublethal dosage of imidacloprid reduces the microglomerular density of honey bee mushroom bodies. Sci. Rep. https://doi.org/10.1038/srep19298 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
de Almeida-Rossi, C., Roat, T. C., Tavares, D. A., Cintra-Socolowski, P. & Malaspina, O. Brain morphophysiology of africanized bee Apis mellifera exposed to sublethal doses of imidacloprid. Arch. Environ. Contam. Toxicol. 65, 234–243 (2013).PubMed 

Google Scholar 
Tosi, S., Burgio, G. & Nieh, J. C. A common neonicotinoid pesticide, thiamethoxam, impairs honey bee flight ability. Sci. Reports 7, 1–8 (2017).
Google Scholar 
Coulon, M. et al. Interactions between thiamethoxam and deformed wing virus can drastically impair flight behavior of honey bees. Front. Microbiol. 0, 766 (2020).
Google Scholar 
Shi, T.-F., Wang, Y.-F., Liu, F., Qi, L. & Yu, L.-S. Sublethal effects of the neonicotinoid insecticide thiamethoxam on the transcriptome of the honey bees (Hymenoptera: Apidae). J. Econ. Entomol. 110, 2283–2289 (2017).CAS 
PubMed 

Google Scholar 
Tesovnik, T. et al. Exposure of honey bee larvae to thiamethoxam and its interaction with Nosema ceranae infection in adult honey bees. Environ. Pollut. 256, 113443 (2020).CAS 
PubMed 

Google Scholar 
Friol, P. S., Catae, A. F., Tavares, D. A., Malaspina, O. & Roat, T. C. Can the exposure of Apis mellifera (Hymenoptera, Apiadae) larvae to a field concentration of thiamethoxam affect newly emerged bees?. Chemosphere 185, 56–66 (2017).ADS 
CAS 
PubMed 

Google Scholar 
Coulon, M. et al. Influence of chronic exposure to thiamethoxam and chronic bee paralysis virus on winter honey bees. PLoS ONE 14, e0220703 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Stoner, K. A. & Eitzer, B. D. Movement of soil-applied imidacloprid and thiamethoxam into nectar and pollen of squash (Cucurbita pepo). PLoS ONE 7, e39114 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Rolke, D., Persigehl, M., Peters, B., Sterk, G. & Blenau, W. Large-scale monitoring of effects of clothianidin-dressed oilseed rape seeds on pollinating insects in northern Germany: residues of clothianidin in pollen, nectar and honey. Ecotoxicology 25, 1691–1701 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Wood, T. J., Kaplan, I., Zhang, Y. & Szendrei, Z. Honeybee dietary neonicotinoid exposure is associated with pollen collection from agricultural weeds. Proc. R. Soc. B 286, 1905 (2019).
Google Scholar 
Urlacher, E. et al. Measurements of chlorpyrifos levels in forager bees and comparison with levels that disrupt honey bee odor-mediated learning under laboratory conditions. J. Chem. Ecol. 42, 127–138 (2016).CAS 
PubMed 

Google Scholar 
Villalba, A., Maggi, M., Ondarza, P. M., Szawarski, N. & Miglioranza, K. S. B. Influence of land use on chlorpyrifos and persistent organic pollutant levels in honey bees, bee bread and honey: Beehive exposure assessment. Sci. Total Environ. 713, 136554 (2020).ADS 
CAS 
PubMed 

Google Scholar 
Traynor, K. S. et al. Pesticides in honey bee colonies: Establishing a baseline for real world exposure over seven years in the USA. Environ. Pollut. https://doi.org/10.1016/j.envpol.2021.116566 (2021).Article 
PubMed 

Google Scholar 
EPA Press Office. EPA Takes Action to Address Risk from Chlorpyrifos and Protect Children’s Health. (2021).Arena, M. & Sgolastra, F. A meta-analysis comparing the sensitivity of bees to pesticides. Ecotoxicology 233(23), 324–334 (2014).
Google Scholar 
Wright, G. A., Nicolson, S. W. & Shafir, S. Nutritional physiology and ecology of honey bees. Annu. Rev. Entomol 63, 327–344. https://doi.org/10.1146/annurev-ento-020117-043423 (2018).CAS 
Article 
PubMed 

Google Scholar 
Michener, C. D. The Social Behavior of the Bees: A Comparative Study (Belknap Press of Harvard University Press, 1974).
Google Scholar 
Stabler, D., Paoli, P. P., Nicolson, S. W. & Wright, G. A. Nutrient balancing of the adult worker bumblebee (Bombus terrestris) depends on the dietary source of essential amino acids. J. Exp. Biol. https://doi.org/10.1242/jeb.114249 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Böhme, F., Bischoff, G., Zebitz, C. P. W., Rosenkranz, P. & Wallner, K. Pesticide residue survey of pollen loads collected by honeybees (Apis mellifera) in daily intervals at three agricultural sites in South Germany. PLoS ONE 13, e0199995 (2018).PubMed 
PubMed Central 

Google Scholar 
Schilder, A. M. C., Hanson, E. J. & Hancock, J. F. An integrated approach to disease control in blueberries in Michigan. Acta Hortic. 715, 481–488 (2006).CAS 

Google Scholar 
Wise, J. C., Jenkins, P. E., Poppen, R. V. & Isaacs, R. Activity of broad-spectrum and reduced-risk insecticides on various life stages of cranberry fruitworm (Lepidoptera: Pyralidae) in Highbush Blueberry. J. Econ. Entomol. 103, 1720–1728 (2010).PubMed 

Google Scholar 
Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2009).MATH 

Google Scholar  More

In this picture, taken in February at Copenhagen Zoo, I’m holding a vacuum device equipped with a tiny fan and filter. The devices — we call them air samplers — are designed to collect DNA samples from the air. We deployed three samplers at the zoo: one in a stable with two okapi (Okapia johnstoni) and two duikers, one in a rainforest house and one outside, near an exhibit of animals that live in the African savannah.At best, we had hoped to detect nearby animals in small enclosures — an okapi in the stable, for instance. But as we reported in Current Biology, the devices outperformed our expectations (C. Lynggaard Curr. Biol. 32, 701–707; 2022). They picked up identifiable DNA from 49 vertebrates, including guppies in the rainforest pool, ostriches and giraffes in the savannah area, and even cats and dogs in the park next door. Interestingly, we didn’t get any signal from turtles in the rainforest house. Maybe turtles mostly keep their DNA to themselves.Our analysis ultimately found that the sampler could detect animals from nearly 200 metres away. The giraffe in the picture is standing much closer than is necessary for collection of a sample.Airborne DNA is all around us. Birds release skin cells when they flap their wings. Saliva from all sorts of animals can become airborne. Animals release DNA when they defecate. In November 2021, I received a grant to start a research group whose goal is to collect airborne DNA in nature. This approach could transform conservation biology and species monitoring. We could detect rare animals and get a better understanding of diversity without disturbing an environment.We have so many lines of inquiry in this work. The location of the samplers, the rates of air flow, the time, the best methods for sorting DNA from the sample — we’re still trying to work all of these out. We hadn’t expected that the zoo experiment would ever work, so we’re scrambling to plan the next steps. It’s an exciting time. More

Argonema gen. nov. Skoupý et Dvořák.Type species: Argonema galeatum.Morphology: Filamentous cyanobacterium, colonies macroscopic, growing in round bulbs and tufts. The filaments are dark green to blue-green, grey-green or brown-green in color. Cells are wider than they are long. Filaments sheathed, sheaths are colorless to light brown, distinct, and variable in length. The filament can protrude from the sheath or the sheath can exceed filament. Trichomes are cylindrical, not attenuated to slightly attenuated towards the end, slightly or not constricted at cell walls. The apical cell can be concave, dark brown, purple-brown to almost black. Cell content often granulated. Necridic cells present, reproduction by hormogonia. The morphological description was based on both culture and fresh material.Etymology: The genus epithet (Argonema) is derived from greek Argo – slow, latent (αργός) and nema – thread (νήμα).A. galeatum sp. nov. Skoupý et Dvořák.Morphology: The cells of A. galeatum are 6.5–9.1 µm (mean 7.81 µm) wide and 1.1–2.5 µm (mean 1.83 µm) long (Figs. 1–5). Filaments are straight, blue-green to gray-green in color. The sheaths are colorless to light brown, distinct, and variable in length. The filament can protrude from the sheath or the sheath can exceed filament. No true branching was observed. Trichomes are cylindrical, not attenuated or slightly attenuated towards the end, slightly or not constricted at cell walls. Some filaments have a concave apical cell that is dark brown, purple-brown to almost black (Fig. 11b). Cell content often granulated. Reproduction by necridic cells and subsequent breaking of the filaments into hormogonia (Fig. 11a,c). The morphological description was based on both culture and fresh material.Figures 1-8Microphotographs of Argonema galeatum (Figs 1–5) and Argonema antarcticum (Figs. 6–8) Trichomes of A. galeatum appear more straight (Fig 2), while trichomes of A. antarcticum form waves (Fig 6) and loops (Fig 7). Scale = 10 µm, wide arrow = necridic cells, arrowhead = granules, asterisk = colored apical cell, circle = empty sheath.Full size imageFigures 9 and 10Histograms of cell dimensions constructed using PAST software. Fig. 9 – Histogram of cell width frequencies in A. galeatum (blue) and A. antarcticum (red). Fig. 10 – Histogram of cell length frequencies in A. galeatum (blue) and A. antarcticum (red).Full size imageHolotype: 38,057, Herbarium of the Department of Botany (OL), Palacký University Olomouc, Czech Republic.Reference strain: Argonema galeatum A003/A1.Type locality: James Ross Island, Western Antarctica, 63.80589S, 57.92147 W.Habitat: Well-developed soil crust.Etymology: Species epithet A. galeatum was derived from latin galea – helmet.A. antarcticum sp. nov. Skoupý et Dvořák.Morphology: The cells are 7.6–9.2 µm (mean 8.52 µm) wide and 1.2–2.8 µm (mean 1.72 µm) long (Figs. 5–8). Filaments are wavy, gray-green to brown-green in color. The sheaths are colorless to light brown, distinct, and variable in length. The filament can protrude from the sheath or the sheath can exceed filament. No true branching was observed. Trichomes are cylindrical, not attenuated or slightly attenuated towards the end with a concave apical cell, slightly or not constricted at cell walls (Fig. 11d). Necridic cells present (Fig. 11e), reproduction by hormogonia. The morphological description was based on both culture and fresh material.Holotype: 38,058, Herbarium of the Department of Botany (OL), Palacký University, Olomouc, Czech Republic.Reference strain: Argonema antarcticum A004/B2.Type locality: James Ross Island, Western Antarctica, 63.89762S, 57.79743 W.Habitat: Well-developed soil crust.Etymology: Species epithet A. antarcticum was derived from the original sampling site.Morphological variabilityWe used light microscopy to assess the morphology of Argonema from soil crust samples and cultured strains. Argonema is morphologically similar to other Oscillatoriales, such as Lyngbya, Phormidium, and Oscillatoria. In culture, the morphology of A. galeatum and A. antarcticum differed slightly. Filaments of A. antarcticum are wider than cells of A. galeatum, averaging at 8.52 µm (A. galeatum – 7.81 µm). The average cell width/length ratio is 4.54 for A.galeatum and 4.89 for A. antarcticum. The cell width was significantly different between the two species (Nested ANOVA, p  More

Study siteThe study area covers the campus of Suranaree University of Technology (SUT) and its surrounding landscape in Muang, Nakhon Ratchasima, Thailand (14.879° N, 102.018° E; Fig. 1). The university campus covers about 11.2 km2, and comprises a matrix of human modified lands interspersed with mixed deciduous forest fragments (at the onset of this study we identified there were 37 mixed deciduous forest fragments on campus, mean = 7.36 ± 1.48 ha, range = 0.45–45.6 ha [note, “±” is used for standard error throughout the text]). More than 15,000 students are enrolled at SUT, and there are numerous multi-story classrooms, laboratory and workshop buildings, residential housing, parking areas, eating and sports facilities, an elementary school, and a large hospital on the university campus. During the first term of the 2019 school year, 7622 students, as well as numerous SUT staff, lived in on-campus residential areas. The landscape surrounding the university is primarily dominated by agriculture, though there are also patches of less-disturbed areas as well as several densely populated villages and suburban housing divisions among the monoculture plots of upland crops (e.g., cassava, maize, and eucalyptus).Figure 1Study site map illustrating the land-use types spanning the area where the Malayan kraits (Bungarus candidus) were tracked in Muang Nakhon Ratchasima, Nakhon Ratchasima province, Thailand. Map created using QGIS v.3.8.2 (https://qgis.org/) in combination with Inkscape v.1.1.0 (https://inkscape.org/).Full size imageThe study site is located within the Korat Plateau region with an altitude range of 205–285 m above sea level. Northeast Thailand has a tropical climate, and the average daily temperature from 1 January 2018 to 31 December 2020 in Muang Nakhon Ratchasima was 28.29 °C, with daily averages ranging from 19.3 to 34.1 °C38. The region receives an average annual rainfall ranging from 1270 to 2000 mm39. There are three distinct seasons in northeast Thailand: cold, wet, and hot, each are classified by annual changes in temperature and rainfall. Cold season is typically between mid-October and mid-February, hot season is generally from mid-February to May, while the highly unpredictable rainfall of the wet season is predominantly concentrated between the months May to October39,40.Due to the representation of agriculture, semi-urban, and suburban areas with patches of more natural areas all within a relatively small area, we determined the university campus provided an ideal setting to examine how land-use features and human activity influence the movements of B. candidus. Additionally, past studies have indicated northeast Thailand hosts the most bites by B. candidus in Thailand29,33, making sites like ours ideal.Study animalsWe opportunistically sampled Malayan kraits captured as a result of notifications from locals and ad-hoc encounters during transit due to low detectability in visual encounter surveys, in addition to those discovered through unstandardized visual encounter surveys. Upon capture, we collected morphometric data, including snout-vent length (SVL), tail length (TL), mass, and sex (Table 1, Supp. Table 1). We measured body lengths with a tape measure, measured body mass with a digital scale, and determined sex via cloacal probing, all while the snakes were anesthetized via inhaling vaporized isoflurane. We then housed individuals with an SVL > 645 mm and mass > 50 g in plastic boxes (with refugia and water) prior to surgical transmitter implantation by a veterinarian from the Nakhon Ratchasima Zoo. We attempted to minimize the time snakes were in captivity awaiting implantation; however, delays arose due to the veterinarian’s availability, the snake being mid-ecdysis, or the snake having a bolus that needed to pass through the digestive tract before implantation (n = 21 implantations, mean = 5.02 ± 0.61 days, range = 0.60–13.02 days). The Nakhon Ratchasima Zoo veterinarian implanted radio transmitters (1.8 g BD-2 or 3.6 g SB-2 Holohil Inc, Carp, Canada) into the coelomic cavity using procedures described by Reinert and Cundall41, while the snake was anesthetized. We assigned each individual an ID according to sex and individual detection number (e.g., M02 = a male was the second B. candidus individual documented during the study). We released the implanted individuals as close as possible to their capture locations (mean = 65.31 m ± 13.7 m, range = 0–226.42 m), though on six occasions we moved individuals ≥ 100 m because the individual came from either residential areas or a busy road (all but one were moved  800 mm; thus, nine of the males were adults and four were juveniles (though two of the males had an SVL > 720 mm, and therefore likely sub-adults). The single telemetered female was an adult.Individual tracking durations varied (mean = 106.46 ± 15.36 days, range = 28.5–222.77 days; Supp. Fig. 1), as many individuals were lost due to unexpected premature transmitter failures (n = 5) or unsuccessful recapture efforts due to individuals sheltering under large buildings as the transmitter reached the end of its battery life (n = 4). We only recorded one confirmed mortality in the study, M01, who was killed by a motorized vehicle when crossing a road (n = 1). Another three individuals were lost due to unknown reasons, which may have been due to premature transmitter failure, mortality, or the animal moving beyond radio signal despite extensive search efforts. Thus, we only successfully recaptured and re-implanted five individuals (M01 once, M02 twice, M07 once, M27 once, and M33 twice). Transmitter batteries generally lasted approximately 90–110 days, so we aimed to replace transmitters after ≥ 90 days of use. At the end of the study, only one individual was successfully recaptured to remove the transmitter.Data collectionWe used very high frequency radio-telemetry to locate each telemetered individual on average every 24.20 h (SE ± 0.41, 0.17–410.0 h; see Supp. Fig. 2 for distribution of tracking time lags). We aimed to locate each individual’s shelter locations once each day during the daylight (06:00–18:00 h); however, we were occasionally (n = 34 days) unable to locate a snake for several consecutive days when we were unable to obtain radio signal due to an individual having moved far away or deep underneath a large structure. There were also a few occasions where we were unable to track snakes due to prolonged and heavy rainfall (n = 4 days), as the moisture damages equipment, or other reasons (n = 4 days). We additionally located snakes nocturnally (18:00–06:00 h) ad hoc and in an attempt to observe nocturnal behaviors and movement pathways when animals were active. We defined fixes as any time a telemetered individual was located, and relocations (i.e., moves) as the occasions where we located an individual > 5 m from its previous known location.Each day we manually honed in on signal via a radio receiver to locate individuals (as described by Amelon et al.42, and recorded locations in Universal Transverse Mercator (UTM; 47 N World Geodetic System 84) coordinate reference system with a handheld global positioning system (GPS) unit (Garmin 64S GPS, Garmin International, Inc., Olathe, Kansas) directly above the sheltered snake. We generally approached within one meter of sheltering snakes during daylight to precisely record shelter locations and identify shelter type. Since we could not visually confirm snake locations, we methodically eliminated all possible locations where the snake could possibly be while at close range with the minimum possible gain on the radio receiver.Telemetered kraits tended to be inactive and sheltering underground during the daylight, thus we were confident that our diurnal location checks would not affect their movements. However, in some cases we resorted to determining an individual’s location via triangulation, where multiple lines cast from different vantage points towards the snake intersect on the snake’s location on the GPS, allowing us to determine the animal’s coordinate location from approximately 10–30 m away. This helped ensure that we recorded locations with greater accuracy when snakes sheltered underneath large buildings, as it allowed us to move away from large structures that hindered the GPS accuracy. This technique was also implemented during some nocturnal location checks when a snake was believed to be active among dense vegetation, in an attempt to prevent disturbance of the animals’ natural behavior. While we did hope to gain visual observations of active individuals during the night, we exercised more caution during nocturnal location checks, typically maintaining a minimum distance of approximately 5 m in attempt to lessen the chances of disturbing an active individual’s behavior. If the animal was active we recorded the animal’s observed behavioral state (i.e., moving, feeding, or foraging). When the radio signal was stable and the individual was not visible, we recorded the animal’s behavior as “sheltering”. We strived for an accuracy of  5 m difference), and land-use type (e.g., mixed deciduous forest, human-settlement, semi-natural area, agriculture, plantation; see Supp. Figs. 3 and 4 for photos of land-use types), behavior (e.g., sheltering, moving, foraging, or feeding), and shelter type (e.g., anthropogenic, burrow, or unknown, note we also recorded if we suspected the shelter to be part of a termite tunnel complex due to a close proximity to a visible termite mound; Supp. Fig. 5).During each location check we recorded the straight-line distance between the current and previous locations (distance moved/step length) with the GPS device. We then used step-lengths to summarize their movements by estimating the mean daily displacement (MDD; the total distance moved divided by the number of days the snake was located) and mean movement distance (MMD; the mean relocation distance, excludes distances ≤ 5). In order to limit biases due to some snakes being located multiple times within a given day/night, we limited our sample for estimating MMD and MDD to only include a single location per day. This was accomplished by manually removing “extra” nocturnal location checks that occurred within the same day, making sure to have all shelter relocations present within the dataset. When calculating MDD, we used the total number of daily location checks rather than the number of days between the individual’s tracking start and stop date since there were some days where individuals were not tracked. We also used the same one location check per day dataset to calculate movement/relocation probabilities and to examine each individual’s MMD, MDD, and relocation probability for the overall tracking duration as well as for each season.When feasible, we positioned a Bushnell (Bushnell Corporation, Overland Park, Kansas) time lapse field camera (Trophy Cam HD Essential E3, Model:119837) with infrared night capability on a tripod spaced 2–5 m from occupied shelter sites. We positioned the cameras so that we may gather photos of the focal snake as it exited the shelter site and/or behaviors exhibited near the shelter. We programmed the cameras using a combined setting, including field scan, which continuously captured one photo every minute, along with a motion sensor setting, which took photos upon movement trigger outside of the regular 1-min intervals.Space use and site fidelityAll analyses and most visualizations were done in R v.4.0.5 using RStudio v.1.4.1106 43,44. We attempted to estimate home ranges for the telemetered B. candidus individuals using autocorrelated kernel density estimates (AKDEs) using R package ctmm v.0.6.045,46 in order to better understand the spatial requirements of B. candidus. However, examination of the variograms revealed that the majority of the variograms had not fully stabilized (i.e., limited evidence of range stability in our sample), and many individuals had extremely low effective sample sizes (21.82 ± 9.75, range = 1.49–135.75; Supp. Table 4). Therefore, we do not report home ranges in this text, as the AKDE estimates would violate the assumption of range residency and either underestimate or misrepresent B. candidus spatial requirements. We also examined the speed estimates resulting from fitted movement models. Resulting variograms and tentative home range estimates are included in a supplementary file for viewing only (Supp. Fig. 6, Supp. Table 4). The original code is from Montaño et al.47.Since our data was not sufficient to estimate home range size for the telemetered B. candidus, we instead used Dynamic Brownian Bridge Movement Models (dBBMMs) with the R package move v.4.0.648 to estimate within study occurrence distributions. We caution readers that these are not home range estimates but instead modeling the potential movement pathways animals could have traversed49. Use of dBBMMs not only allows us to estimate occurrence distributions for each individual, thus helping us better understand the animal’s movement pathways and resource use, but it also allows us to examine movement patterns through dBBMM derived motion variance50,51. We selected a window size of 19 and margin size of 5, to catch short resting periods with the margin, while the window size of 19 is long enough to get a valid estimate of motion variance when the animals exhibit activity/movement. Contours however are somewhat arbitrary; therefore, we used three different contours levels (90%, 95%, 99%) to estimate dBBMM occurrence distributions (using R packages adehabitatHR v.0.4.19, and rgeos v.0.5.5), and show the sensitivity to contour choice52,53.All movement data, either including initial capture locations or beginning with the first location check ~ 24 h post release, was used for production of both the AKDEs and dBBMMs for each individual. We also estimated dBBMM occurrence distributions for each telemetered individual with the exception of M29, which only made three small moves within a burrow complex during the short time he was radio-tracked before transmitter failure.We compared space use estimates to two previously published B. candidus tracking datasets34,36, and one unpublished dataset shared on the Zenodo data repository54, all originating from the Sakaerat Biosphere Reserve (approximately 41 km to the south of our study site): two adult males from within the forested area of the reserve [one tracked every 27.8 ± 0.99 h over a period of 103 days, the other tracked every 38.63 ± 11.2 h over a period of 30.58 days]34,54, and a juvenile male from agriculture on a forest boundary [tracked every 50.19 ± h for 66.91 days]36. The previous studies on B. candidus only tracked the movements of a single individual each, had coarser tracking regimes, and used traditional—fundamentally flawed methods55,56—to estimate space use34,36. Therefore, we ran dBBMMs with these previous datasets using the same window (19) and margin size (5).To quantify site reuse and time spent at sites (residency time) we used recursive analysis with the R package recurse v.1.1.257. We defined each site as a circular area with a radius of 5 m around each unique location (matching the targeted GPS accuracy). Then we calculated each individual’s overall number of relocations, each individual’s total number of relocations to each site, and each individual’s site revisit frequency and residency time at each unique site. Then we plotted revisited locations on a land-use map with space use estimates (95% and 99% dBBMM) in an attempt to help identify and highlight activity centers for telemetered individuals (see Supp. Figs. 7–13). All maps were created using Quantum Geographic Information System (QGIS v.3.8.2).Habitat selectionWe used Integrated Step Selection Function models (ISSF) to examine the influence of land-use features on the movements of B. candidus at both the individual and population levels. We included movement data from all male individuals that used more than one habitat feature in our ISSF analysis. Therefore, we excluded F16 and M29 who both only used settlement habitat. Excluding M29 was justified by the individual having been tracked for the shortest duration (19 days) and had the fewest number of moves (n = 3), thus there were not enough relocations for ISSF models to work effectively. Using modified code from Smith et al.51 that used ISSF with Burmese python radio-telemetry data, we used the package amt v.0.1.458 to run ISSF for each individual, with Euclidean distance to particular land-use features (natural areas, agriculture, settlement, buildings, and roads) to determine association or avoidance of features. Cameron Hodges created all land-use shape files in QGIS by digitizing features from satellite imagery and verified all questionable satellite land-use types via on-ground investigation.The semi-natural areas, plantations, mixed deciduous forest and water bodies (such as irrigation canals and ponds which have densely vegetated edges) were all combined into a single layer of less-disturbed habitats which we refer to as “natural areas”. All feature raster layers were then converted into layers with a gradient of continuous values of Euclidean distances to the land-use features, and were inverted in order to avoid zero-inflation of distance to feature values and to make the resulting model directional effects easier to more intuitive. We were able to generate 200 random steps per each observed step (following Smith et al.51), due to the coarse temporal resolution of manually collected radio-telemetry data (i.e., we were not computational limited when deciding the number of random locations). Higher numbers of random steps are preferable as they can aid in detecting smaller effects and rarer landscape features59.To investigate individual selection, we created nine different models testing for association to habitat features, with one being a null model which solely incorporated step-length and turning angle to predict movement60, five examining land-use features individually (agriculture, buildings, settlement, natural areas, roads), and the other three being multi-factor models. Each model considers distance to a land-use variable, step-length, and turn-angle as an aspect of the model. After running each of the nine models for each individual, we then examined the AIC for each model, point estimates (with lower and upper confidence intervals), and p-values in order to identify the best models for each individual and determine the strongest relationships and trends among the samples. We considered models with ∆ AIC  More