Ecology

1.
Hewitt, G. M. Some genetic consequences of ice ages. Biol. J. Lin. Soc. 58, 247–276 (1996).
Article  Google Scholar 
2.
Husemann, M., Schmitt, T., Zachos, F. E., Ulrich, W. & Habel, J. C. Palearctic biogeography revisited: evidence for the existence of a North African refugium for Western Palaearctic biota. J. Biogeogr. 41, 81–94. https://doi.org/10.1111/jbi.12180 (2014).
Article  Google Scholar 

3.
Hewitt, G. M. Postglacial recolonization of European biota. Biol. J. Lin. Soc. 68, 87–112 (1999).
Article  Google Scholar 

4.
Hewitt, G. M. Genetic consequences of climatic oscillations in the Quaternary. Philos. Trans. R. Soc. B Biol. Sci. 359, 183–195. https://doi.org/10.1098/rstb.2003.1388 (2004).
CAS  Article  Google Scholar 

5.
de Lattin, G. Grundriß der Zoogeographie (Gustav Fischer, Jena, 1967).
Google Scholar 

6.
Schmitt, T. & Varga, Z. Extra-Mediterranean refugia: The rule and not the exception?. Front. Zool. 9, 1–12. https://doi.org/10.1186/1742-9994-9-22 (2012).
Article  Google Scholar 

7.
Drees, C. et al. Molecular analyses and species distribution models indicate cryptic northern mountain refugia for a forest-dwelling ground beetle. J. Biogeogr. 43, 2223–2236. https://doi.org/10.1111/jbi.12828 (2016).
Article  Google Scholar 

8.
Juřičková, L., Horáčková, J. & Ložek, V. Direct evidence of central European forest refugia during the last glacial period based on mollusc fossils. Quat. Res. 82, 222–228. https://doi.org/10.1016/j.yqres.2014.01.015 (2014).
Article  Google Scholar 

9.
Pinceel, J., Jordaens, K., Pfenninger, M. & Backeljau, T. Rangewide phylogeography of a terrestrial slug in Europe: Evidence for Alpine refugia and rapid colonization after the Pleistocene glaciations. Mol. Ecol. 14, 1133–1150. https://doi.org/10.1111/j.1365-294X.2005.02479.x (2005).
CAS  Article  PubMed  Google Scholar 

10.
Gratton, P., Konopiński, M. K. & Sbordoni, V. Pleistocene evolutionary history of the Clouded Apollo (Parnassius mnemosyne): Genetic signatures of climate cycles and a “time-dependent” mitochondrial substitution rate. Mol. Ecol. 17, 4248–4262. https://doi.org/10.1111/j.1365-294X.2008.03901.x (2008).
CAS  Article  PubMed  Google Scholar 

11.
Hofman, S. et al. Phylogeography of the fire-bellied toads Bombina: Independent Pleistocene histories inferred from mitochondrial genomes. Mol. Ecol. 16, 2301–2316. https://doi.org/10.1111/j.1365-294X.2007.03309.x (2007).
CAS  Article  PubMed  Google Scholar 

12.
Junker, M. et al. Three in one-multiple faunal elements within an endangered european butterfly species. PLoS ONE 10, 1–24. https://doi.org/10.1371/journal.pone.0142282 (2015).
CAS  Article  Google Scholar 

13.
Ursenbacher, S., Carlsson, M., Helfer, V., Tegelström, H. & Fumagalli, L. Phylogeography and Pleistocene refugia of the adder (Vipera berus) as inferred from mitochondrial DNA sequence data. Mol. Ecol. 15, 3425–3437. https://doi.org/10.1111/j.1365-294X.2006.03031.x (2006).
CAS  Article  PubMed  Google Scholar 

14.
Magri, D. Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica). J. Biogeogr. 35, 450–463. https://doi.org/10.1111/j.1365-2699.2007.01803.x (2008).
Article  Google Scholar 

15.
Svenning, J.-C., Normand, S. & Kageyama, M. Glacial refugia of temperate trees in Europe: Insights from species distribution modelling. J. Ecol. 96, 1117–1127. https://doi.org/10.1111/j.1752-4598.2012.00212.x (2008).
Article  Google Scholar 

16.
Cheddadi, R. et al. Imprints of glacial refugia in the modern genetic diversity of Pinus sylvestris. Global Ecol. Biogeogr. 15, 271–282. https://doi.org/10.1111/j.1466-822x.2006.00226.x (2006).
Article  Google Scholar 

17.
Willis, K. J. & Van Andel, T. H. Trees or no trees? The environments of central and eastern Europe during the Last Glaciation. Quatern. Sci. Rev. 23, 2369–2387. https://doi.org/10.1016/j.quascirev.2004.06.002 (2004).
ADS  Article  Google Scholar 

18.
Rudner, Z. E. & Sümegi, P. Recurring Taiga forest-steppe habitats in the Carpathian Basin in the Upper Weichselian. Quatern. Int. 76, 177–189. https://doi.org/10.1016/S1040-6182(00)00101-4 (2001).
Article  Google Scholar 

19.
van Swaay, C., Warren, M. & Grégoire, L. Biotope use and trends of European butterflies. J. Insect Conserv. 10, 189–209. https://doi.org/10.1007/s10841-006-6293-4 (2006).
Article  Google Scholar 

20.
Slamova, I., Klecka, J. & Konvicka, M. Diurnal behavior and habitat preferences of Erebia aethiops, an aberrant lowland species of a mountain butterfly clade. J. Insect Behav. 24, 230–246. https://doi.org/10.1007/s10905-010-9250-8 (2011).
Article  Google Scholar 

21.
Burnaz, S. & Balazs, S. Contributions to the knowledge of diurnal Lepidoptera fauna of the North-Eastern part of Ţarcu Mountains (Southern Carpathians, Romania). Buletin de Informare Entomologica 22, 41–52 (2011).
Google Scholar 

22.
Slamova, I., Klecka, J. & Konvicka, M. Woodland and grassland mosaic from a butterfly perspective: Habitat use by Erebia aethiops (Lepidoptera: Satyridae). Insect Conserv. Divers. 6, 243–254. https://doi.org/10.1111/j.1752-4598.2012.00212.x (2013).
Article  Google Scholar 

23.
Tshikolovets, V. V. Butterflies of Europe & the Mediterranean area (Tshikolovets Publications, Pardubice, 2011).
Google Scholar 

24.
Varga, Z. Die Erebien der Balkanhalbinsel und Karpaten V. Übersicht der subspezifischen Gliederung und der Verbreitung der Erebia Dalman, 1816 -Arten (Lepidoptera, Nymphalidae, Saryrinae) in der Balkanhalbinsel und in den Karpaten (II. Teil). Entomol Roman 6, 5–39 (2001).
Google Scholar 

25.
Imbrie, J., et al. Milankovitch theory, the two shorter cycles can be explained radiation cycle (arising from of Geological Sciences, Brown 2 Institut d’ Astronomie et de Geophysique of Earth, Atmospheric, and Planetary Earth Observatory, Columbia for Climatic and Sp. Cycle 8, 699–735 (1993).

26.
Nei, M. Genetic Distance between Populations 106, 283–29 (1972).

27.
Hausdorf, B. & Hennig, C. Species delimitation using dominant and codominant multilocus. Markers 59, 491–503. https://doi.org/10.1093/sysbio/syq039 (2010).
CAS  Article  Google Scholar 

28.
Nakatani, T., Usami, S. & Itoh, T. Molecular phylogenetic analysis of the Erebia aethiops group (Lepidoptera, Nymphalidae). Lepidopterol. Soc. Jpn. 58, 387–404 (2007).
Google Scholar 

29.
Lukhtanov, V. & Lukhtanov, A. Die Tagfalter Nordwestasiens (Dr. Ulf Eitschberger, Marktleuthen, 1994).
Google Scholar 

30.
Peña, C., Witthauer, H., Klečková, I., Fric, Z. & Wahlberg, N. Adaptive radiations in butterflies: Evolutionary history of the genus Erebia (Nymphalidae: Satyrinae). Biol. J. Lin. Soc. 116(2), 449–467. https://doi.org/10.1111/bij.12597 (2015).
Article  Google Scholar 

31.
Hampe, A. & Petit, R. J. Conserving biodiversity under climate change: The rear edge matters. Ecol. Lett. 8, 461–467. https://doi.org/10.1111/j.1461-0248.2005.00739.x (2005).
Article  PubMed  Google Scholar 

32.
Schmitt, T., Hewitt, G. M. & Müller, P. Disjunct distributions during glacial and interglacial periods in mountain butterflies: Erebia epiphron as an example. J. Evol. Biol. 19, 108–113. https://doi.org/10.1111/j.1420-9101.2005.00980.x (2006).
CAS  Article  PubMed  Google Scholar 

33.
Schmitt, T. & Seitz, A. Intraspecific allozymatic differentiation reveals the glacial refugia and the postglacial expansions of European Erebia medusa (Lepidoptera: Nymphalidae). Biol. J. Linn. Soc. 74, 429–458. https://doi.org/10.1006/bijl.2001.0584 (2001).
Article  Google Scholar 

34.
van Husen, D. Geological processes during the Quaternary. Mitteilungen Der Österreichischen Geologischen Gesellschaft 92, 135–156 (2000).
Google Scholar 

35.
Schmitt, T. Biogeographical and evolutionary importance of the European high mountain systems. Front. Zool. 6, 1–10. https://doi.org/10.1186/1742-9994-6-9 (2009).
Article  Google Scholar 

36.
Albre, J., Gers, C. & Legal, L. Molecular phylogeny of the Erebia tyndarus (Lepidoptera, Rhopalocera, Nymphalidae, Satyrinae) species group combining CoxII and ND5 mitochondrial genes: A case study of a recent radiation. Mol. Phylogenet. Evol. 47, 196–210. https://doi.org/10.1016/j.ympev.2008.01.009 (2008).
CAS  Article  PubMed  Google Scholar 

37.
Fišer Pečnikar, Ž, Balant, M., Glasnović, P. & Surina, B. Seed dormancy and germination of the rare, high elevation Balkan endemic Cerastium dinaricum (Caryophyllaceae). Biologia 73, 937–943. https://doi.org/10.2478/s11756-018-0115-5 (2018).
CAS  Article  Google Scholar 

38.
Giachino, P. M. A new species of Aphaobiella Pretner, 1949 from Grintavec Mt., Slovenia (Coleoptera: Cholevidae, Leptodirinae). Fragm. Entomol. 48, 19–23. https://doi.org/10.4081/fe.2016.165 (2016).
Article  Google Scholar 

39.
Canestrelli, D., Salvi, D., Maura, M., Bologna, M. A. & Nascetti, G. One species, three Pleistocene evolutionary histories: Phylogeography of the Italian crested newt, Triturus carnifex. PLoS ONE 7, e41754. https://doi.org/10.1371/journal.pone.0041754 (2012).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

40.
Previšić, A., Walton, C., KuČiniĆ, M., Mitrikeski, P. T. & Kerovec, M. Pleistocene divergence of Dinaric Drusus endemics (Trichoptera, Limnephilidae) in multiple microrefugia within the Balkan Peninsula. Mol. Ecol. 18, 634–647. https://doi.org/10.1111/j.1365-294X.2008.04046.x (2009).
CAS  Article  PubMed  Google Scholar 

41.
Deffontaine, V. et al. Beyond the Mediterranean peninsulas: Evidence of central European glacial refugia for a temperate forest mammal species, the bank vole (Clethrionomys glareolus). Mol. Ecol. 14, 1727–1739. https://doi.org/10.1111/j.1365-294X.2005.02506.x (2005).
CAS  Article  PubMed  Google Scholar 

42.
Krebs, P., Pezzatti, G. B., Beffa, G., Tinner, W. & Conedera, M. Revising the sweet chestnut (Castanea sativa Mill.) refugia history of the last glacial period with extended pollen and macrofossil evidence. Quat. Sci. Rev. 206, 111–128. https://doi.org/10.1016/j.quascirev.2019.01.002 (2019).
ADS  Article  Google Scholar 

43.
Varga, Z. Geographische Isolation und Subspeziation bei den Hochgebirgslepidopteren der Balkanhalbinsel. Acta Entomol Jugoslaviae. (1975).

44.
Bhagwat, S. A. & Willis, K. J. Species persistence in northerly glacial refugia of Europe: A matter of chance or biogeographical traits?. J. Biogeogr. 35, 464–482. https://doi.org/10.1111/j.1365-2699.2007.01861.x (2008).
Article  Google Scholar 

45.
Filipi, K., Marková, S., Searle, J. B. & Kotlík, P. Mitogenomic phylogenetics of the bank vole Clethrionomys glareolus, a model system for studying end-glacial colonization of Europe. Mol. Phylogenet. Evol. 82, 245–257. https://doi.org/10.1016/j.ympev.2014.10.016 (2015).
Article  PubMed  Google Scholar 

46.
Hammouti, N., Schmitt, T., Seitz, A., Kosuch, J. & Veith, M. Combining mitochondrial and nuclear evidences: A refined evolutionary history of Erebia medusa (Lepidoptera: Nymphalidae: Satyrinae) in Central Europe based on the COI gene. J. Zool. Syst. Evolut. Res. 48, 115–125. https://doi.org/10.1111/j.1439-0469.2009.00544.x (2010).
Article  Google Scholar 

47.
Vila, M., Marí-Mena, N., Guerrero, A. & Schmitt, T. Some butterflies do not care much about topography: a single genetic lineage of Erebia euryale (Nymphalidae) along the northern Iberian mountains. J. Zool. Syst. Evolut. Res. 49, 119–132. https://doi.org/10.1111/j.1439-0469.2010.00587.x (2011).
Article  Google Scholar 

48.
Vodă, R., Dapporto, L., Dincă, V. & Vila, R. Cryptic matters: Overlooked species generate most butterfly beta-diversity. Ecography 38, 405–409. https://doi.org/10.1111/ecog.00762 (2015).
Article  Google Scholar 

49.
Hebert, P. D. N., Ratnasingham, S. & DeWaard, J. R. Barcoding animal life: Cytochrome c oxidase subunit 1 divergences among closely related species. Proc. R. Soc. B Biol. Sci. 270, 1–4. https://doi.org/10.1098/rsbl.2003.0025 (2003).
CAS  Article  Google Scholar 

50.
Weller, S. J., Pashley, D. P., Martin, J. A. & Constable, J. L. Phylogeny of noctuoid moths and the utility of combining independent nuclear and mitochondrial genes. Syst. Biol. 43, 194–211. https://doi.org/10.1093/sysbio/43.2.194 (1994).
Article  Google Scholar 

51.
Hebert, P. D. N., & Beaton, M. J. Methodologies for allozyme analysis using cellulose acetate electrophoresis. Zoology 32 (1993).

52.
Kearse, M. et al. Geneious basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649. https://doi.org/10.1093/bioinformatics/bts199 (2012).
Article  PubMed  PubMed Central  Google Scholar 

53.
Thompson, J. D., Higgins, D. G. & Gibson, T. J. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22, 4673–4680. https://doi.org/10.1093/nar/22.22.4673 (1994).
CAS  Article  PubMed  PubMed Central  Google Scholar 

54.
Hall, T. A. BioEdit. Nucleic Acids Symp. Ser. 41, 95–98 (1999).
CAS  Google Scholar 

55.
Rozas, J. et al. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol. Biol. Evol. 34, 3299–3302. https://doi.org/10.1093/molbev/msx248 (2017).
CAS  Article  Google Scholar 

56.
Bandelt, H. J., Forster, P. & Röhl, A. Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. https://doi.org/10.1093/oxfordjournals.molbev.a026036 (1999).
CAS  Article  Google Scholar 

57.
Leigh, J. W. & Bryant, D. POPART: Full-feature software for haplotype network construction. Methods Ecol. Evol. 6, 1110–1116. https://doi.org/10.1111/2041-210X.12410 (2015).
Article  Google Scholar 

58.
Bouckaert, R. et al. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 15, e1006650. https://doi.org/10.1371/journal.pcbi.1006650 (2019).
CAS  Article  PubMed  PubMed Central  Google Scholar 

59.
Lanfear, R., Frandsen, P. B., Wright, A. M., Senfeld, T. & Calcott, B. Partitionfinder 2: New methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol. Biol. Evol. 34, 772–773. https://doi.org/10.1093/molbev/msw260 (2017).
CAS  Article  PubMed  Google Scholar 

60.
Papadopoulou, A., Anastasiou, I. & Vogler, A. P. Revisiting the insect mitochondrial molecular clock: The mid-aegean trench calibration. Mol. Biol. Evol. 27, 1659–1672. https://doi.org/10.1093/molbev/msq051 (2010).
CAS  Article  PubMed  Google Scholar 

61.
Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901–904. https://doi.org/10.1093/sysbio/syy032 (2018).
CAS  Article  PubMed  PubMed Central  Google Scholar 

62.
Drummond, A. J., & Ho, S. Y. W. A rough guide to BEAST 1.4. Edinburgh, 1–41, http://beast-mcmc.googlecode.com/files/BEAST14_Manual_17May2007.pdf (2007).

63.
Yu, Y., Harris, A. J., Blair, C. & He, X. RASP (Reconstruct Ancestral State in Phylogenies): A tool for historical biogeography. Mol. Phylogenet. Evol. 87, 46–49. https://doi.org/10.1016/j.ympev.2015.03.008 (2015).
Article  PubMed  Google Scholar 

64.
Siegismund, H. R. G-STAT, Version 3, Genetical Statistical Programs for the Analysis of Population Data. (The Arboretum, Royal Veterinary and Agricultural University, Horsholm, Denmark, 1993).

65.
Excoffier, L., Laval, G. & Schneider, S. Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evol. Bioinform. 1, 117693430500100. https://doi.org/10.1177/117693430500100003 (2005).
Article  Google Scholar 

66.
Rice, W. R. The sequential Bonferroni test. Evolution 43, 235 (1989).
Article  Google Scholar 

67.
Saitou, N. & Nei, M. The neighbor-joining method: A new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 97, 407–425 (1987).
Google Scholar 

68.
Felsenstein, J. PHYLIP (Phylogeny Inference Package), Version 3.5.c. (Department of Genetics, University of Washington, Seattle, Washington, 2000).

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

70.
Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 14, 2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x (2005).
CAS  Article  PubMed  Google Scholar 

71.
Guillot, G., Estoup, A., Mortier, F., & Cosson, J. F. A spatial statistical model for landscape genetics. Genetics 170, 1261–1280. https://doi.org/10.1534/genetics.104.033803 (2005).
CAS  Article  PubMed  PubMed Central  Google Scholar 

72.
Team, Q. D. QGIS Geographic Information System. Open Source Geospatial Foundation Project (2020). http://qgis.osgeo.org. Accessed 05 Nov 2020. More

1.
Ohmann, B. & Babiuk, L. A. Viral infections in domestic animals as models for studies of viral immunology and pathogenesis. J. Gen. Virol. 66, 1–25 (1986).
Article  Google Scholar 
2.
Woolhouse, M., Scott, F., Hudson, Z., Howey, R. & Chase-Topping, M. Human viruses: discovery and emergence. Phil. Trans. R. Soc. B. 367, 2864–2871 (2012).
PubMed  Article  PubMed Central  Google Scholar 

3.
Johnson, C. K. et al. Spillover and pandemic properties of zoonotic viruses with high host plasticity. Sci. Rep. 5, 14830 (2015).
ADS  Article  CAS  Google Scholar 

4.
Glennon, E. E. et al. Domesticated animals as hosts of henipaviruses and filoviruses: a systematic review. Veterinary J. 233, 25–34 (2018).
Article  Google Scholar 

5.
Letarov, A. & Kulikov, E. The bacteriophages in human- and animal body-associated microbial communities. J. Appl. Microbiol. 107, 1–13 (2009).
CAS  PubMed  Article  PubMed Central  Google Scholar 

6.
Shkoporov, A. N. & Hill, C. Bacteriophages of the human gut: the “known unknown” of the microbiome. Cell Host Microbe 25, 95–209 (2019).
Article  CAS  Google Scholar 

7.
Kwok, K. T. T., Nieuwenhuijse, D. F., Phan, M. V. T. & Koopmans, M. P. G. Virus metagenomics in farm animals: a systematic review. Viruses 12, E107 (2020).
PubMed  Article  CAS  PubMed Central  Google Scholar 

8.
Grasis, J. A. et al. Species-specific viromes in the ancestral holobiont Hydra. PLoS ONE 9, e109952 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

9.
Bosch, T. C., Grasis, J. A. & Lachnit, T. Microbial ecology in Hydra: Why viruses matter. J. Microbiol. 53, 193–200 (2015).
PubMed  Article  PubMed Central  Google Scholar 

10.
Weynberg, K. D. et al. Coral-associated viral communities show high levels of diversity and host auxiliary functions. PeerJ 5, e4054 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

11.
Mahmoud, H. & Jose, L. Phage and nucleocytoplasmic large viral sequences dominate coral viromes from the Arabian Gulf. Front. Microbiol. 8, 2063 (2017).
PubMed  PubMed Central  Article  Google Scholar 

12.
López-Madrigal, S. & Duarte, E. H. Titer regulation in arthropod-Wolbachia symbioses. FEMS Microbiol. Lett. 366, fnz232 (2019).
PubMed  Article  CAS  Google Scholar 

13.
Bordenstein, S. R., Marshall, M. L., Fry, A. J., Kim, U. & Wernegreen, J. J. The tripartite associations between bacteriophage, Wolbachia, and arthropods. PLoS Pathog. 2, e43 (2006).
PubMed  PubMed Central  Article  Google Scholar 

14.
Perlmutter, J. I. et al. The phage gene wmk is a candidate for male killing by a bacterial endosymbiont. PLoS Pathog. 15, e1007936 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

15.
Johnson, P. T. Viral diseases of marine invertebrates. Helgolander Meeresunters 37, 65–98 (1984).
Article  Google Scholar 

16.
Weinbauer, M. G. Ecology of prokaryotic viruses. FEMS Microbiol. Rev. 28, 127–181 (2004).
CAS  PubMed  Article  Google Scholar 

17.
Munn, C. B. Viruses as pathogens of marine organisms—from bacteria to whales. J. Mar. Biol. Ass. 86, 453–467 (2006).
Article  Google Scholar 

18.
Lang, A. S., Rise, M. L., Culley, A. I. & Steward, G. F. RNA viruses in the sea. FEMS Microbiol. Rev. 33, 295–323 (2009).
CAS  PubMed  Article  Google Scholar 

19.
Rosario, K., Schenck, R. O., Harbeitner, R. C., Lawler, S. N. & Breitbart, M. Novel circular single-stranded DNA viruses identified in marine invertebrates reveal high sequence diversity and consistent predicted intrinsic disorder patterns within putative structural proteins. Front. Microbiol. 6, 696 (2015).
PubMed  PubMed Central  Article  Google Scholar 

20.
Reuter, M. Viruslike particles in Gyratrix hermaphroditus (Turbellaria: Rhabdocoela). J. Invertebr. Pathol. 25, 79–95 (1975).
CAS  PubMed  Article  PubMed Central  Google Scholar 

21.
Vijayan, K. K. et al. Polychaete worms—a vector for white spot syndrome virus (WSSV). Dis. Aquat. Org. 63, 107–111 (2005).
CAS  Article  Google Scholar 

22.
Nobiron, I. et al. Genome and polypeptides characterization of Tellina virus 1 reveals a fifth genetic cluster in the Birnaviridae family. Virology 371, 350–361 (2008).
CAS  PubMed  Article  PubMed Central  Google Scholar 

23.
Crespo-González, C. et al. Virus-like particles in Urastoma cyprinae, a turbellarian parasite of Mytilus galloprovincialis. Dis. Aquat. Organ. 79, 83–86 (2008).
PubMed  Article  PubMed Central  Google Scholar 

24.
Marhaver, K. L., Edwards, R. A. & Rohwer, F. Viral communities associated with healthy and bleaching corals. Environ. Microbiol. 10, 2277–2286 (2008).
CAS  PubMed  PubMed Central  Article  Google Scholar 

25.
Claverie, J.-M. et al. Mimivirus and Mimiviridae: giant viruses with an increasing number of potential hosts, including corals and sponges. J. Invertebr. Pathol. 101, 172–180 (2009).
CAS  PubMed  Article  PubMed Central  Google Scholar 

26.
Jackson, E. W., Bistolas, K. S. I., Button, J. B. & Hewson, I. Novel circular single-stranded DNA viruses among an asteroid, echinoid and holothurian (Phylum: Echinodermata). PLoS ONE 11, e0166093 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

27.
Suttle, C. A. Viruses in the sea. Nature 437, 356–361 (2005).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

28.
Brum, J. R., Schenck, R. O. & Sullivan, M. B. Global morphological analysis of marine viruses shows minimal regional variation and dominance of non-tailed viruses. The ISME J 7, 1738–1751 (2013).
CAS  PubMed  Article  PubMed Central  Google Scholar 

29.
Middelboe, M. & Brussaard, C. Marine viruses: key players in marine ecosystems. Viruses 9, E302 (2017).
PubMed  Article  PubMed Central  Google Scholar 

30.
Vacelet, J. & Gallissian, M.-F. Virus-like particles in the cells of the sponge Verongia cavernicola (Demospongiae, Dictyoceratida) and accompanying tissues changes. J. Invertebr. Pathol. 31, 246–254 (1978).
Article  Google Scholar 

31.
Luter, H. M., Whalan, S. & Webster, N. S. Exploring the role of microorganisms in the disease-like syndrome affecting the sponge Lanthella basta. Appl. Environ. Microbiol. 76, 5736–5744 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

32.
Pascelli, C., Laffy, P. W., Kupresanin, M., Ravasi, T. & Webster, N. S. Morphological characterization of virus-like particles in coral reef sponges. PeerJ 6, e5625 (2018).
PubMed  PubMed Central  Article  CAS  Google Scholar 

33.
Lohr, J. E., Chen, F. & Hill, R. T. Genomic analysis of bacteriophage JL001: insights into its interaction with a sponge-associated alpha-Proteobacterium. Appl. Environ. Microbiol. 71, 1598–1609 (2005).
CAS  PubMed  PubMed Central  Article  Google Scholar 

34.
Lohr, J., Munn, C. B. & Wilson, W. H. Characterization of a latent virus-like infection of symbiotic zooxanthellae. Appl. Environ. Microbiol. 73, 2976–2981 (2007).
CAS  PubMed  PubMed Central  Article  Google Scholar 

35.
Patten, N. L., Harrison, P. L. & Mitchell, J. G. Prevalence of virus-like particles within a staghorn scleractinian coral (Acropora muricata) from the Great Barrier Reef. Coral Reefs 27, 569–580 (2008).
ADS  Article  Google Scholar 

36.
van Oppen, M. H., Leong, J.-A. & Gates, R. D. Coral-virus interactions: A double-edged sword?. Symbiosis 47, 1–8 (2009).
Article  Google Scholar 

37.
Vega Thurber, R. L. & Correa, A. M. S. Viruses of reef-building scleractinian corals. J. Exp. Mar. Biol. Ecol. 408, 102–113 (2011).
Article  Google Scholar 

38.
Leruste, A., Bouvier, T. & Bettarel, Y. Enumerating viruses in coral mucus. Appl. Environ. Microbiol. 78, 6377–6379 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

39.
Pollock, F. J. et al. Abundance and morphology of virus-like particles associated with the coral Acropora hyacinthus differ between healthy and white syndrome-infected states. Mar. Ecol. Prog. Ser. 510, 39–43 (2014).
ADS  Article  Google Scholar 

40.
Correa, A. M. S. et al. Viral outbreak in corals associated with an in situ bleaching event: atypical herpes-like viruses and a new megavirus infecting Symbiodinium. Front. Microbiol. 7, 127 (2016).
PubMed  PubMed Central  Article  Google Scholar 

41.
Farley, C. A. Viruses and viruslike lesions in marine mollusks. Mar. Fish. Rev. 40, 18–20 (1978).
Google Scholar 

42.
Elston, R. A. Bivalve mollusc viruses. World J. Microbiol. Biotechnol. 13, 393–403 (1997).
Article  Google Scholar 

43.
Renault, T. & Novoa, B. Viruses infecting bivalve molluscs. Aquat. Living Resour. 17, 397–409 (2004).
Article  Google Scholar 

44.
Richards, G. P. et al. Bacteriophages against pathogenic vibrios in Delaware Bay oysters (Crassostrea virginica) during a period of high levels of pathogenic Vibrio parahaemolyticus. Food Environ. Virol. 11, 101–112 (2019).
PubMed  Article  PubMed Central  Google Scholar 

45.
Leigh, B., Karrer, C., Cannon, J. P., Breitbart, M. & Dishaw, L. J. Isolation and characterization of a Shewanella phage–host system from the gut of the tunicate Ciona intestinalis. Viruses 9, E60 (2017).
PubMed  Article  CAS  PubMed Central  Google Scholar 

46.
Winston, J. E. Feeding in marine bryozoans. In Biology of Bryozoans (eds Woollacott, R. M. & Zimmer, R. L.) 233–271 (Academic Press, London, 1977).
Google Scholar 

47.
Winston, J. E. Polypide morphology and feeding behaviour in marine ectoprocts. Bull. Mar. Sci. 28, 1–31 (1978).
Google Scholar 

48.
Shunatova, N. N. & Ostrovsky, A. N. Individual autozooidal behaviour and feeding in marine bryozoans. Sarsia 86, 113–142 (2001).
Article  Google Scholar 

49.
Shunatova, N. & Ostrovsky, A. Group autozooidal behaviour and chimneys in marine bryozoans. Mar. Biol. 140, 503–518 (2002).
Article  Google Scholar 

50.
Schwaha, T. F., Ostrovsky, A. N. & Wanninger, A. Key novelties in the evolution of Bryozoa: evidence from the soft-body morphology. Biol. Rev. 95, 696–729 (2020).
PubMed  Article  PubMed Central  Google Scholar 

51.
Ryland, J. S. Bryozoans (Hutchinson University Library, London, 1970).
Google Scholar 

52.
Ryland, J. S. Bryozoa: an introductory overview. In Moostiere (Bryozoa). Denisia Vol. 16 (ed. Woess, E.) 9–22 (Springer, Linz, 2005).
Google Scholar 

53.
McKinney, F. K. & Jackson, J. D. C. Bryozoan Evolution (Unwin Hyman, Boston, MA, 1989).
Google Scholar 

54.
Nielsen, C. Bryozoa (Ectoprocta: ‘Moss’ animals). eLS 1, 1–6. https://doi.org/10.1002/9780470015902.a0001613.pub2 (2013).
Article  Google Scholar 

55.
Lutaud, G. Sur la présence de microorganismes spécifiques dans les glandes vestibulaires et dans l’aviculaire de Palmicellaria skenei (Ellis et Solander) Bryozoaire Chilostome. Cah. Biol. Mar. 6, 181–190 (1965).
Google Scholar 

56.
Lutaud, G. La nature des corps funiculaires des cellularines Bryozoaires Chilostomes. Arch. Zool. Exp. Gen. 110, 2–30 (1969).
Google Scholar 

57.
Lutaud, G. L’infestation du myoépithélium de l’oesophage par des microorganismes pigmentés et la structure des organes à bactéries du vestibule chez le Bryozoaire Chilostome Palmicellaria skenei (E. et S.) Can. J. Zool. 64, 1842–1851 (1986).
Google Scholar 

58.
Woollacott, R. M. & Zimmer, R. L. A simplified placenta-like system for the transport of extraembryonic nutrients during embryogenesis of Bugula neritina (Bryozoa). J. Morphol. 147, 355–377 (1975).
PubMed  Article  Google Scholar 

59.
Dyrynda, P. E. J. & King, P. E. Sexual reproduction in Epistomia bursaria (Bryozoa: Cheilostomata), an endozooidal brooder without polypide recycling. J. Zool. 198, 337–352 (1982).
Article  Google Scholar 

60.
Moosbrugger, M., Schwaha, T., Walzl, M. G., Obst, M. & Ostrovsky, A. N. The placental analogue and the pattern of sexual reproduction in the cheilostome bryozoan Bicellariella ciliata (Gymnolaemata). Front. Zool. 9, 29 (2012).
PubMed  PubMed Central  Article  Google Scholar 

61.
Mathew, M., Schwaha, T., Ostrovsky, A. N. & Lopanik, N. B. Symbiont-dependent sexual reproduction in marine colonial invertebrate: Morphological and molecular evidence. Mar. Biol. 165, 14 (2018).
Article  CAS  Google Scholar 

62.
Karagodina, N. P., Vishnyakov, A. E., Kotenko, O. N., Maltseva, A. L. & Ostrovsky, A. N. Ultrastructural evidence for nutritional relationships between a marine colonial invertebrate (Bryozoa) and its bacterial symbionts. Symbiosis 75, 155–164 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

63.
Woollacott, R. M. Association of bacteria with bryozoan larvae. Mar. Biol. 65, 155–158 (1981).
Article  Google Scholar 

64.
Zimmer, R. L. & Woollacott, R. M. Mycoplasma-like organisms: occurrence with the larvae and adults of a marine bryozoan. Science 220, 208–210 (1983).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

65.
Zimmer, R. L. & Woollacott, R. M. Larval morphology of the bryozoan Watersipora arcuata (Cheilostomata: Ascophora). J. Morphol. 199, 125–150 (1989).
PubMed  Article  PubMed Central  Google Scholar 

66.
Boyle, P. J., Maki, J. S. & Mitchell, R. Mollicute identified in novel association with aquatic invertebrate. Curr. Microbiol. 15, 85–89 (1987).
CAS  Article  Google Scholar 

67.
Lim, G. E. & Haygood, M. G. “Candidatus Endobugula glebosa”, a specific bacterial symbiont of the marine bryozoan Bugula simplex. Appl. Environ. Microbiol. 70, 4921–4929 (2004).
CAS  PubMed  PubMed Central  Article  Google Scholar 

68.
Sharp, K. H., Davidson, S. K. & Haygood, M. G. Localization of ‘Candidatus Endobugula sertula’ and the bryostatins throughout the life cycle of the bryozoan Bugula neritina. ISME J. 1, 693–702 (2007).
PubMed  Article  PubMed Central  Google Scholar 

69.
Lim-Fong, G. E., Regali, L. A. & Haygood, M. G. Evolutionary relationships of “Candidatus Endobugula” bacterial symbionts and their Bugula bryozoan hosts. Appl. Environ. Microbiol. 74, 3605–3609 (2008).
CAS  PubMed  PubMed Central  Article  Google Scholar 

70.
Richardson, K. C., Jarrett, L. & Finke, E. H. Embedding in epoxy resins for ultrathin sectioning in electron microscopy. Stain. Technol. 35, 313–323 (1960).
CAS  PubMed  Article  Google Scholar 

71.
Mäntynen, S., Sundberg, L.-R., Oksanen, H. M. & Poranen, M. M. Half a century of research on membrane-containing bacteriophages: bringing new concepts to modern virology. Viruses 11, E76 (2019).
PubMed  Article  CAS  Google Scholar 

72.
Cuozzo, S. A., Castellano, P., Sesma, F. J., Vignolo, G. M. & Raya, R. R. Differential roles of the two-component peptides of lactocin 705 in antimicrobial activity. Curr. Microbiol. 46, 180–183 (2003).
CAS  PubMed  Article  Google Scholar 

73.
Chibani-Chennoufi, S., Bruttin, A., Dillmann, M. L. & Brüssow, H. Phage-host interaction: An ecological perspective. J. Bacteriol. 186, 3677–3686 (2004).
CAS  PubMed  PubMed Central  Article  Google Scholar 

74.
Haygood, M. G. & Davidson, S. K. Small-subunit rRNA genes and in situ hybridization with oligonucleotides specific for the bacterial symbionts in the larvae of the bryozoan Bugula neritina and proposal of “Candidatus Endobugula sertula”. Appl. Environ. Microbiol. 63, 4612–4616 (1997).
CAS  PubMed  PubMed Central  Article  Google Scholar 

75.
Ostrovsky, A. N., Gordon, D. P. & Lidgard, S. Independent evolution of matrotrophy in the major classes of Bryozoa: Transitions among reproductive patterns and their ecological background. Mar. Ecol. Prog. Ser. 378, 113–124 (2009).
ADS  Article  Google Scholar 

76.
Ostrovsky, A. N. Evolution of Sexual Reproduction in Marine Invertebrates: Example of Gymnolaemate Bryozoans (Springer, Dordrecht, Heidelberg, NewYork, London, 2013).
Google Scholar 

77.
Ostrovsky, A. N. From incipient to substantial: evolution of placentotrophy in a phylum of aquatic colonial invertebrates. Evolution 67, 1368–1382 (2013).
PubMed  PubMed Central  Article  Google Scholar 

78.
Miller, I. J., Vanee, N., Fong, S. S., Lim-Fong, G. E. & Kwan, J. C. Lack of overt genome reduction in the bryostatin-producing bryozoan symbiont “Candidatus Endobugula sertula”. Appl. Environ. Microbiol. 82, 6573–6583 (2016).
CAS  PubMed  PubMed Central  Article  Google Scholar 

79.
Miller, I. J., Weyna, T. R., Fong, S. S., Lim-Fong, G. E. & Kwan, J. C. Single sample resolution of rare microbial dark matter in a marine invertebrate metagenome. Sci. Rep. 6, 34362 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

80.
Ostrovsky, A. N. & Porter, J. S. Pattern of occurrence of supraneural coelomopores and intertentacular organs in Gymnolaemata (Bryozoa) and its evolutionary implications. Zoomorphology 130, 1–15 (2011).
Article  Google Scholar 

81.
Shikuma, N. J., Pilhofer, M., Weiss, G. L., Hadfield, M. G. & Jensen, G. J. Marine tubeworm metamorphosis induced by arrays of bacterial phage tail–like structures. Science 343, 529–533 (2014).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

82.
Leiman, P. G. et al. Type VI secretion apparatus and phage tail-associated protein complexes share a common evolutionary origin. Proc. Nat. Acad. Sci. 106, 4154–4159 (2009).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

83.
Taylor, N. M. I., van Raaij, M. J. & Leiman, P. G. Contractile injection systems of bacteriophages and related systems. Mol. Microbiol. 108, 6–15 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

84.
Lawrence, S. A., Wilson, W. H., Davy, J. E. & Davy, S. K. Latent virus-like infections are present in a diverse range of Symbiodinium spp. (Dinophyta). J. Phycol. 50, 984–997 (2014).
PubMed  Article  PubMed Central  Google Scholar 

85.
Ericson, C. F. et al. A contractile injection system stimulates tubeworm metamorphosis by translocating a proteinaceous effector. eLife 8, e46845 (2019).
PubMed  PubMed Central  Article  Google Scholar 

86.
Ackermann, H. W. Basic phage electron microscopy. Methods Mol. Biol. 501, 113–126 (2009).
CAS  PubMed  Article  PubMed Central  Google Scholar 

87.
Moran, N. A., Degnan, P. H., Santos, S. R., Dunbar, H. E. & Ochman, H. The players in a mutualistic symbiosis: insects, bacteria, viruses, and virulence genes. Proc. Natl. Acad. Sci. USA 102, 16919–16926 (2005).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

88.
Weldon, S. R. & Oliver, K. M. Diverse bacteriophage roles in an aphid-bacterial defensive mutualism. In Advances in Environmental Microbiology Vol. 2 (ed. Hurst, C. J.) 173–206 (Springer, Berlin, 2016).
Google Scholar 

89.
Rohwer, F. & Thurber, R. V. Viruses manipulate the marine environment. Nature 459, 207–212 (2009).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

90.
Roossinck, M. J. The good viruses: viral mutualistic symbioses. Nat. Rev. Microbiol. 9, 99–108 (2011).
CAS  PubMed  Article  PubMed Central  Google Scholar 

91.
Sharp, J. H., Winson, M. K. & Porter, J. S. Bryozoan metabolites: an ecological perspective. Nat. Prod. Rep. 24, 659–673 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

92.
Maltseva, A. L. et al. Novel brominated metabolites from Bryozoa: a functional analysis. Nat. Prod. Res. 31, 1840–1848 (2016).
PubMed  Article  CAS  PubMed Central  Google Scholar  More

Clothianidin exposure in field experiments
The experimental bee yard was situated in a suburban area in Kirchhain (Germany) with sufficient natural pollen sources in late summer, where the honeybees were allowed to forage freely. The field experiment took place with 20 Apis mellifera carnica colonies in Kirchhainer mating nuclei (25.5 × 19.8 × 17 cm), established with sister queens (A. m. carnica breeding line, mother: 17:27:20:11) mated on an island mating station (Norderney) and 180 g of worker bees each. We tried to reduce the variance between colonies by using young sister queens from an island mating station. Furthermore, we used small hives to minimize food competition so that we were able to set up the colonies within a rather small area of 20 by 20 m. In detail, the worker bees originated from four healthy colonies of the institute, located at an apiary of the institute. In order to get a homogenous composition, worker bees of brood combs and honeycombs of two colonies were shaken into a box and mixed thoroughly, before they were uniformly distributed into 10 mating nuclei. This procedure was repeated to fill 20 mating nuclei. Subsequently, the mating boxes were placed in a cool room (12 °C) for three day, to acclimatize, before they were transported to the island mating station Norderney. The queens started to lay eggs at June 25, 2014. The clothianidin exposure started at July 28. Thus, the colonies containing freshly mated, egg-laying queens were allowed to establish and build up an intact brood nest for four weeks. Two days before clothianidin exposure (sampling day 0, S0; SFig. 1), the colony strength was assessed and the treatment groups were randomly assigned to the hives such that differences between treatment groups were minimized with respect to the strength of their colonies. Every week, from July 30 to September 10, each colony received 400 mL of Apiinvert (Südzucker AG, Mannheim, Germany) sugar syrup (39% w/v fructose, 31% w/v saccharose and 30% w/v glucose) spiked with clothianidin (1, 10 or 100 µg/L). Control colonies received sugar syrup containing the same concentration of the solvent (water) as the clothianidin-treated groups. The syrup was fed in zip-look bags placed inside the food chamber containing a climbing aid. After 1 week, the leftovers were removed and weighed to record food consumption. For all four experimental groups, the analyzed clothinidin levels were close to the target concentrations (Suppl. Table 1). Environmental data were recorded during the study period using a USB data logger (EL-USB-2, Lascar Electronics Ltd., temperature accuracy ± 0.5 °C, relative humidity accuracy ± 3%) located under the colonies.
Sampling
During weeks 3 and 7, random samples of worker bees, larvae and worker jelly were taken from each colony. In detail, 20 randomly chosen worker bees located on a brood comb and five larvae (larval stage: day 7 or 8 after egg laying) were immediately frozen for chemical analysis. To collect worker jelly, extra thick blotting paper (Protean, Xi size; Bio-Rad, Hercules, CA, USA) was cut into strips, cleaned in pure ethanol and acetone (Carl Roth, Karlsruhe, Germany) and dried in a heating cabinet at 80 °C. Each strip was inserted into five brood cells containing a small larva (developmental day 4–5) to suck up the worker jelly and the strips were immediately frozen.
To document brood development, each side of each comb was photographed within an empty hive box transformed into a photo box containing a digital camera (Canon PowerShot A1000 IS, Tokyo, Japan) and a ring-flash (Aputure Amaran AHL-C60 LED, Shenzhen, China). Brood documentation took place every week. The colonies were sampled starting with the control and then from the lowest to highest concentration of clothianidin to minimize the risk of carryover.
HPG size measurements
To obtain worker bees of a defined age, single frames of late-stage capped brood (Binder, Tuttlingen, Germany) were brought to the laboratory and incubated in the dark at 32 °C, with humidity provided by open water jars. The frames with worker brood were collected from two full-sized colonies, which were regularly inspected for symptoms of disease and tested for Chronic bee paralysis virus, Deformed wing virus, Acute bee paralysis virus, and Sac brood virus65. Newly-emerged bees (≤ 24 h) were collected, color marked, and transferred to the experimental colonies on July 28 (15 marked bees per colony). During the second week of exposure, the marked bees (Suppl. Table 2) were removed from the colonies after 12 days in the hive and immobilized on ice. The HPGs were dissected in ice-cold phosphate-buffered saline (PBS, pH 7.4). The specimens were fixed in formaldehyde (4% in PBS, Carl Roth), rinsed three times in PBS, and mounted in Aquapolymount (Polysciences, Eppelheim, Germany). Three pictures of each gland (only one gland per bee) were photographed at 400x magnification using a phase contrast/fluorescence microscope (Leica DMIL, Leica camera DFC 420C) and LAS v4.4 image-capturing software (Leica Microsystems, Wetzlar, Germany). To measure the size of each gland, the diameters of 30 acini per bee were measured using ImageJ v1.49o (http://rsb.info.nih.gov/ij/index.html).
High-performance thin-layer chromatography
HPTLC was chosen as sensitive analytical method in order to detect qualitative and semi-quantitative differences in the composition of brood food and larvae of dosed hives. Individual larvae (one per colony) differing in their weights (200–500 mg/larva) were solely macerated and extracted in 1 mL n-hexane ( > 99% pure, Rotisolv, Carl Roth) in an ultrasonic bath for 1 min and then vortexed for 1 min. For the analysis of worker jelly (from three cells per colony), adsorptive filter strips (Sugi strips, Kettenbach, Eschenburg, Germany) were cut in half and one part was dipped in brood combs until maximal absorption of the material, whereas the other strip was used as background control strip, and both strips were extracted as above. The supernatants of larvae and worker jelly were transferred to a fresh vial 1000 μL isopropylacetate/methanol (3/2, v/v). After maceration for 1 min, the mix was vortexed for 1 min and the supernatant was transferred to another vial. Between extractions, the samples were cooled on ice and stored at –20 °C.
The isopropylacetate/methanol extracts (20 µL/band for worker jelly and 7 µL/band for larvae) were sprayed onto the silica gel 60 F254 HPTLC plate (Merck, Darmstadt, Germany) using an Automatic TLC Sampler 4. The plate of worker jelly was developed with a mobile phase consisting of chloroform/methanol/water/ammonia (30/17/2/1, v/v/v/v, all Carl Roth66 and the plate of larvae with an 8-step gradient development based on methanol, chloroform, toluene and n-hexane66,67. After drying in a stream of cold air for 2 min, the plate images were documented at UV 366 nm using the TLC Visualizer. For derivatisation, the chromatogram was dipped into the primuline solution (100 mg primuline in 200 mL acetone/water, 4/1 v/v, Sigma-Aldrich, Steinheim, Germany) at an immersion speed of 2.5 cm/s and an immersion time of 1 s, dried and derivatised as before. The data were processed using winCATS version 1.4.2.8121 (all instrumentation from CAMAG).
The bacterium Aliivibrio fischeri (NRRL-B11177, strain 7151), obtained from the German Collection of Microorganisms and Cell Cultures (DSMZ, Leibniz Institute, Berlin, Germany), was used to assess a non-targeted, broad range of effective substances within the worker jelly. HPTLC plates were developed as described above, neutralized and immersed in a bacterial suspension, prepared according to DIN EN ISO 11,348–1 845 at an immersion speed of 3 cm/s and an immersion time of 2 s. The bioluminescence of the wet bioautogram was recorded in an interval of 3 min over 30 min using the BioLuminizer (CAMAG).
Brood assessment
The free extension of the open source program ImageJ, Fiji (http://fiji.sc/Fiji) was used to count all cells with eggs, larvae, or sealed brood, for every colony on every sampling day. The photos of a single comb from different sampling days were aligned in the program and all brood cells were counted. Because the colonies had no wax foundations, some cells on the edges of these naturally-built combs were at an unfavorable angle. Therefore, most but presumably not all cells with brood were visible in the photos. This uncertainty was similar in all hives. To estimate the brood survival, we first tracked the development of individual eggs, but found a high mortality rate even in control colonies. Therefore, we tracked the development of individual larvae (day 4–5 after egg laying) over 4 weeks (= sampling weeks). Brood survival was estimated for two periods during the experiment (weeks 1–4 and 3–7).
Modeling of demographic compensation
The aim of the model was to estimate the number of new larvae that must hatch each week in order to maintain a stable number of larvae up to 7 weeks of age given the survivorship schedule of larvae week by week. Let h denote the number of eggs that hatch into larvae each week, and let the probability that any individual larva dies in each successive week be mi, where i takes values in the set {1, 2, …, 7} to indicate each of seven successive weeks, after which we assume that surviving larvae pupate. We used a demographic matrix model68 to describe the state of the population of larvae each week as follows. Let lx denote the number of larvae in the colony that are aged x weeks post-hatching and let mx denote their per capita weekly mortality rate. The population of larvae is distributed into seven age classes and we also assign a class to queens, which give rise to larvae by producing eggs. Using the Lefcovitch matrix approach, the larval population of a colony can thus be viewed as a state vector nt whose week-by-week change is the product of a matrix A and the population state vector nt, as shown in Eq. (1):

$$An_{t} = left[ {begin{array}{*{20}c} 0 & 0 & 0 & ldots & 0 & {varvec{h}} \ {(1 – {varvec{m}}_{1} )} & 0 & 0 & ldots & 0 & 0 \ 0 & {left( {1 – {varvec{m}}_{2} } right)} & 0 & ldots & 0 & 0 \ vdots & vdots & vdots & ddots & vdots & vdots \ 0 & 0 & 0 & ldots & {left( {1 – {varvec{m}}_{7} } right)} & 0 \ 0 & 0 & 0 & ldots & 0 & 1 \ end{array} } right]left[ {begin{array}{*{20}c} {{varvec{l}}_{1} } \ {{varvec{l}}_{2} } \ {{varvec{l}}_{3} } \ vdots \ {{varvec{l}}_{7} } \ {varvec{Q}} \ end{array} } right]$$
(1)

The lowermost element of nt is the number of queens (Q) which is here set to Q = 1 for all models, and the total number of larvae in the colony at any time is ({varvec{L}}_{{varvec{t}}} = mathop sum nolimits_{{varvec{x}}} {varvec{l}}_{{varvec{x}}}). Given the observed mortality schedule (mx), we used the model to solve for the value of h that produces a stable value for Lt, as shown in Eq. (2):

$$An_{t} = n_{t}$$
(2)

We determined the mortality schedule by observing the survivorship of a larval cohort. For example, if a cohort of St larvae was originally marked at time t and, of these, St+1 survived until the following week, the per capita weekly mortality rate would be estimated as shown in Eq. (3):

$${varvec{m}}_{{varvec{t}}} = 1 – frac{{{varvec{S}}_{{{varvec{t}} + 1}} }}{{{varvec{S}}_{{varvec{t}}} }}$$
(3)

BEEHAVE simulations
To predict the impact of clothianidin on colony development in standard hives we used BEEHAVE, a honeybee model that simulates colony dynamics and agent-based foraging in realistic landscapes44(http://beehave-model.net/). Although BEEHAVE does not explicitly allow the incorporation of pesticides, the effect of pesticides on behavior and mortality can nevertheless be addressed59,60. BEEHAVE simulates the development of a single honeybee colony, starting with 10,000 foragers on January 1. Colony dynamics are based on a daily egg laying rate, with the developmental stages eggs, larvae, pupae, and adults (in drones) or in-hive workers and foragers (in workers). The brood needs to be tended by in-hive bees, and the larvae additionally need to be fed with nectar and pollen. Foragers can scout for new food sources or collect nectar and pollen from sources already known. Successful foragers can recruit nestmates to the food source. Mortality rates depend on the developmental stage and the time spent on foraging. The colony dies if it either runs out of honey or if the colony size falls below 4000 bees at the end of the year. Swarming may take place when the brood nest grows to more than 17,000 bees before July 18. Under default conditions, two food sources are present at distances of 500 or 1,500 m. Daily foraging conditions are based on weather data from Rothamsted, UK.
We ran two sets of simulations: (A) We first set up BEEHAVE to mimic our experimental nucleus colonies. We then determined how the protein content of the jelly produced by nurses (ProteinFactorNurses) had to be modified to replicate the larval mortality we observed in our empirical data. (B) We then set up BEEHAVE under default conditions but modified ProteinFactorNurses according to the results from the previous simulations to assess the impact of clothianidin exposure under more realistic conditions.
To mimic the experimental nucleus colonies (simulation A), the maximum honey store (MAX_HONEY_STORE_kg) was reduced to 0.77 kg and the maximum size of the brood nest (MAX_BROODCELLS) was reduced to 250. Furthermore, HoneyIdeal was set to ‘true’, so that even though the honey store was small it was filled every day, reflecting the feeding of the experimental bees. In contrast, PollenIdeal was set to ‘false’, because the experimental colonies still had to forage for pollen. On day 209 (July 28), we set the number of pupae to 100 and the number of workers to 650, similar to the experimental colony sizes. During the exposure to clothianidin between days 211 (July 30) and 253 (September 10), the protein content of the jelly fed to the larvae (ProteinFactorNurses) was modified by the new variable ProteinNursesModifier_Exposed. We tested for ProteinNursesModifier_Exposed values from 0.6 to 1 in steps of 0.01. The main output of the simulation was the survival of the brood, calculated from the brood cohort sizes aged 19 days divided by the sizes of these cohorts when they were 3 days old. We calculated the mean brood survival over 30 replicates, using the last 10 cohorts only (i.e. those reaching the age of 19 days between September 1 and 10). Those parameter values for ProteinNursesModifier_Exposed resulting in brood survival most similar to the experimental brood survival were then chosen to represent the clothianidin concentrations of 100, 10 and 1 µg/L.
To assess the impact of clothianidin on standard colonies under more realistic conditions (simulation B), we ran BEEHAVE under default settings but reduced the protein content of the jelly (ProteinFactorNurses) during times of exposure. We assumed that colonies would be exposed when rapeseed plants are flowering, defined in the model as the period between days 95 (April 5) and 130 (May 10). ProteinNursesModifier_Exposed values representing the tested concentrations of clothianidin were derived from the previous set of simulations, and for the control we set ProteinNursesModifier_Exposed to 1 (i.e. no effect of clothianidin). Swarming was either prevented or allowed, in which case the simulation followed the colony remaining in the hive.
Statistical methods
Statistical analysis was carried out using R v3.4.269, including the add-on packages lme470 for linear mixed-effects models, pbkrtest71 for testing fixed effects in mixed-effects models, parallel69 to increase computational power, RLRsim72 for testing random effects in mixed-effects model, multcomp73 for multiple comparisons, and lattice74 for various graphical displays. We used linear mixed-effects models for one- and two-factorial analysis of variance (ANOVA) or regressions as indicated and where necessary. In those models, the colonies were modeled as random effects to reflect the (longitudinal) grouping structure in the data. For the analysis of larval survival, we used a logarithmic transformation of the proportion of surviving larvae. When testing fixed effects in mixed-effects models, we used the Kenward-Roger method (and double-checked the results by comparing them with parametric bootstrap values). Where sufficient, we simplified the analysis using linear (fixed-effects) models. Model diagnostics were performed for all fitted models using qualitative tools such as normal q-q-plots for residuals and plots of residuals versus fitted values to assess the validity of model assumptions like homoscedastic normality. Dunnett’s test (or customized contrasts as appropriate) was used for multiple comparisons in post hoc analysis, and P values were appropriately adjusted for multiple testing within well-defined test families using the single-step or Westfall’s method. Model and analysis details, model diagnostic graphs, and further information are available in the supplemental statistical report. More

A total of 779 persons opened the survey, 486 persons started responding to it, and 308 persons from 40 countries submitted their responses. Seven of these 308 respondents had never heard about biases, and the remaining 301 persons (i.e. 38.6% of those who opened the survey) answered our questions about their attitude to biases.
Nearly all (98%) scientists who responded to our survey were aware of the importance of biases in science. Among these, 33% reported ‘very well’ for their awareness on this topic, 52% classified their awareness as ‘well’ and 13% as ‘poor’. Most of respondents learned about biases from their university courses (36%), from contacts with colleagues (22%) or from scientific literature (20%). Among the different kinds of biases, the best known was observer/observation bias (82%), followed by publication bias (71%) and selection bias (70%); confirmation, reporting/presentation, researcher, measurement, geographic and funding biases were known to 50–60% of respondents (see Appendix S1). Among the seven suggested definitions of biases (see Appendix S1), ‘the tendency to search for, interpret, and publish information in a way that confirms one’s pre-existing beliefs or hypotheses’ corresponded to the understanding of biases by 80% of the respondents, and ‘preferential publication of statistically significant results’ was ranked the second, being selected by 61% of respondents.
In the opinion of our respondents, the stages of scientific research differed in their susceptibility to biases (χ220 = 90.0, P  planning/designing the study  > publishing the outcomes  > reporting the outcomes  > analysing the results  > implementing the study. Similarly, the types of publications were considered prone to biases to different extents (χ216 = 86.8, P  studies based on observational data  > studies based on modelling  > studies based on experiments  > meta-analyses.
Most scientists thought that the severity of the impact of biases on science in general and on their particular research field was medium or high, and a few considered that it was negligible. At the same time, our respondents estimated the impact of biases on their own studies as high almost three times less frequently and as negligible seven times more frequently when compared with their estimates of the impact of biases on other studies within their own research field (Fig. 1).
Figure 1

Impact of biases on science in general, on one’s own research field and on one’s own studies, as estimated by 301 respondents. Bars marked with different letters differ from each other at P  More

A unified framework model of the H/P ratio
Based on the Lotka–Volterra equations24,25, the biomass dynamics of primary producers (P) and herbivores (H) are described as follows:

$$begin{array}{l}dP/dt = gleft( P right)P – xP – fleft( P right)PH,\ dH/dt = kfleft( P right)PH – mHend{array}$$
(1)

where g(P) is biomass-specific primary production (gC gC−1 d−1) and may be a function of P (gC m−2) owing to density-dependent growth, f(P) is per capita grazing rate of herbivores (m2 gC−1 d−1) and also may be function of P depending on the functional response, x is biomass-specific loss rate of primary producers other than due to grazing loss (gC gC−1 d−1), k is the conversion efficiency of herbivores as a fraction of ingested food converted into herbivore biomass (dimensionless: 0~1), and m is per capita mortality rate of herbivores owing to predation and other factors (gC gC−1 d−1). A list of model variables is listed in Supplementary Table 1. If we assume both g(P) and f(P) are constants, Eq. (1) is basically the Lotka–Volterra model, whereas it is an expansion of the Rosenzweig–MacArthur model if we assume logistic growth for g(P) and Michaelis–Menten (Holling type II) functional responses for f(P)24. At the equilibrium state, i.e., dP/dt = 0 and dH/dt = 0, the abundance of producers (P*) and consumers (H*) can be represented as:

$$H ast /P ast = left{ {left[ {gleft( {P ast } right) – x} right]} right./left. {left. {fleft( {P ast } right)} right]} right}/left{ {{mathrm{m/}}left[ {kfleft( {P ast } right)} right]} right}$$
(2)

Thus, the relationship between H and P is not affected by the types of the functional response in herbivores (f(P)). If we set g(P) as the biomass-specific primary production rate at equilibrium, as in simple Lotka–Volterra equations (i.e., g(P*) = g), then the H/P ratio can be expressed with log transformation as:

$${mathrm{log}}left( {H ast /P ast } right) = {mathrm{log}}left( k right) + {mathrm{log}}(g – x) – {mathrm{log}}left( m right)$$
(3)

At equilibrium, (g − x)P* is the amount of primary production that herbivores consume per unit of time (f(P*)P*H*). Thus, if this amount is divided by primary production per unit of time (P*g), it corresponds to the fraction of primary production that herbivores consume (0~1). We define it as β (= 1 − x/g). Large β values imply that producers are efficiently grazed at the equilibrium state. Thus, β is a gauge of inefficiency in the producers’ defensive traits. Using these parameters, the H*/P* ratio can be expressed as:

$${mathrm{log}}left( {H ast /P ast } right) = {mathrm{log}}left( k right) + {mathrm{log}}left( beta right) + {mathrm{log}}left( g right) – {mathrm{log}}left( m right)$$
(4)

This equation implies that the H*/P* biomass ratio on a log scale is affected additively by the specific primary production rate (log(g)), the grazeable fraction of primary production (log(β)), the conversion efficiency (log(k)), and the mortality rate of herbivores (log(m)). According to this equation, communities with relatively low carnivore abundance would have a correspondingly low value of m and will exhibit high herbivore biomass relative to producer biomass (H*/P*), whereas those with low primary production (with low value of g) owing to, for example, low light supply will have a low H*/P* ratio. An increase in defended producers such as armored plants or a decrease in edible producers will decrease β by increasing the loss rate x owing to the cost of defense, and will result in a decreased H*/P* ratio. Finally, when the nutritional value of producers decreases, the conversion efficiency of herbivores (k) should be low, which in turn decreases the H*/P*biomass ratio.
The model for a test with plankton communities
To apply Eq. (4) to a natural community, some modifications are necessary. Here, we consider a plankton community composed of algae and zooplankton. A theory of ecological stoichiometry suggests that the carbon content of primary producers relative to their nutrient content such as nitrogen or phosphorus is an important property affecting growth efficiency in herbivores4. Supporting the theory, a number of studies have shown that growth rate in terms of carbon accumulation relative to ingestion rate strongly depends on the carbon contents of the food relative to nutrients26,27,28,29. Thus, k can be expressed as:

$$k = q_1 times a_{{mathrm{nut}}}^{varepsilon 1}$$
(5)

where αnut is carbon content relative to nutrient content of primary producers and q1 is the conversion factor adjusting to biomass units. In this study, we applied a power function with coefficient of ε1 as a first order approximation because effects of this factor on the H*/P* biomass ratio may not be proportionally related to plant nutrient content. For example, if ε1 is much smaller than zero, it means that negative effects of the carbon to phosphorus ratio of algal food on an herbivore’s k are more substantial when the carbon to phosphorus ratio is high compared with the case when the carbon to phosphorus ratio is low. However, if this factor does not affect the H*/P* biomass ratio, ε1 = 0 and k is constant.
As herbivorous zooplankton cannot efficiently graze on larger phytoplankton due to gape limitation30, the feeding efficiency of herbivores or the defense efficiency of the producers’ resistance traits, β, would be related to the fraction of edible algae in terms of size as follows:

$$beta = q_2 times a_{{mathrm{edi}}}^{varepsilon 2}$$
(6)

where αedi is a trait determining producer edibility, q2 is a factor for converting the traits to edible efficiency, and ε2 is how effective the trait is in defending against grazing. We expect ε2 = 0 if this factor does not matter in regulating the H*/P* biomass ratio but ε2  > 0 if it has a role. Similarly, g can be described as

$$g = q_3 times mu ^{varepsilon 3}$$
(7)

where μ is the specific growth rate of producers, q3 is a conversion factor, and ε3 is the effect of µ on growth rate. Again, we expect that ε3 ≠ 0 if g has a role in determining the H*/P* ratio. Finally, assuming a Holling type I functional response of carnivores, the mortality rate of herbivores, m, is expressed as:

$$m = q_4 times theta ^{varepsilon 4}$$
(8)

where θ is abundance of carnivores, q4 is specific predation rate, and ε4 is the effect size of carnivore abundance on m.
By inserting Eqs. (5–8) to Eq. (4), effects of factors on the H*/P* biomass ratio is formulated as:

$${!}{mathrm{log}}left( {H ast /P ast } right) {!}= varepsilon _1{mathrm{log}}(a_{{mathrm{nut}}}) + varepsilon _2{mathrm{log}}(a_{{mathrm{edi}}}) + varepsilon _3{mathrm{log}}(mu ) – varepsilon _4{mathrm{log}}(theta ) + gamma$$
(9)

where γ is log(q1) + log(q2) + log(q3) − log(q4). If differences in the H*/P* ratio among communities are regulated by growth rate (μ), edibility (αedi), and nutrient contents (αnut) of producers as well as by predation by carnivores (θ), we expect non-zero values for ε1–ε4. Thus, Eq. (9) can be used to evaluate the relative importance of the four hypothesized agents if all of them simultaneously affect the H*/P* ratio. Here we undertake this analysis using data for natural plankton communities in experimental ponds where primary production rate was manipulated with different abundance of carnivore fish.
Experimental test by plankton communities
The experiment was carried out at two ponds (pond ID 217 and 218) located at the Cornell University Experimental Ponds Facility in Ithaca, NY, USA during 4 June to 28 August 2016 (Fig. 1). Each pond has a 0.09 ha surface area (30 × 30 m) and is 1.5 m deep. To initiate the experiment, we equally divided each of the two ponds into four sections using vinyl-coated canvas curtains, and randomly assigned the four sections to either high-shade (64% shading), mid-shade (47% shading), low-shade (33% shading), or no-shade treatments (no shading). Shading in each treatment was made using opaque floating mats (6 m diameter; Solar-cell SunBlanket, Century Products, Inc., Georgia, USA)31. The floating mats were deployed silvered side up to reflect sunlight and blue side down to avoid pond heating. Sampling was performed biweekly for water chemistry and abundance of phytoplankton and zooplankton with measurements of vertical profiles of water temperature, dissolved oxygen (DO) concentration and photosynthetic active radiation (PAR).
Fig. 1: Ponds used in the experimental test.

Pond 217 (a) and Pond 218 (b) in the Cornell University Experimental Ponds Facility divided into four sections by vinyl-canvas curtains and partially shaded by floating mats to regulate primary production rate. Floating docks were placed at the center of the ponds for sampling.

Full size image

PAR in the water column was lower in the sections with larger shaded areas throughout the experiment (Fig. 2b). Water temperature varied from 18 to 25 °C during the experiment but showed no notable differences in mean values of the water columns or the vertical profiles among the four treatments of the two ponds regardless of the shading treatment (Supplementary Fig. 1a, Supplementary Fig. 2). In all treatments, pH values gradually decreased towards the end of experiment and were higher in pond 218 (Supplementary Fig. 1b). DO concentration varied among the treatments and between the ponds, but were within the range of 5–12 mg L−1 (Supplementary Fig. 1c).
Fig. 2: Relationships between zooplankton biomass and phytoplankton biomass, and between H/P mass ratio and photosynthetic active radiation.

Biplots of zooplankton biomass (μg C L−1) and phytoplankton biomass (mg C L−1) (a) and H/P mass ratio and photosynthetic active radiation (PAR, mol photon m−2 d−1) (b) in the water column during the experiment in no-shade (blue), low-shade (orange), mid-shade (red), and high-shade (gray) treatments in pond 217 (circles) and 218 (squares). In each panel, small symbols denote values at each sampling date, and large symbols denote the mean values among the sampling dates. Bars denote standard errors on the means (n = 7 sampling date in each section). Correlation coefficients (r) with p values between the mean values are inserted in each panel.

Full size image

Phytoplankton biomass (mg C L−1) correlated significantly with chlorophyll a (µg L−1) (r = 0.702, p  More

The concentration of signal molecule ComX increases by the square of bacterial density
In order to determine the dynamics of signal molecule (ComX) production, we have used experimental and mathematical modeling approaches. We quantified the ComX concentration over time in the spent medium of PS-216 (ΔcomP) with the biosensor strain BD2876 (for strain-description see Supplementary Table 1), which produces β-galactosidase in response to the exogenous addition of ComX34. The assay included proper controls and calibrations to assure the biosensor-derived ComX concentrations are accurate (for details see Materials and methods). We found that the ComX concentration correlated positively with population density of PS-216 (ΔcomP) and remained constant at 10 nM after entering the stationary phase (Fig. 2a). Importantly, the representation of ComX concentration versus cell density (OD650) (Fig. 2b) showed a non-linear trend between the two parameters. The experimental data were fitted by an allometric function:

$$SMleft( t right) = aNleft( t right)^b$$
(1)

Where SM(t) is a signal molecule (ComX) concentration in time, N(t) is bacterial cell density in time, expressed as optical density of the bacterial suspension (OD650). The fitting results for parameters a and b were (9.6 ± 0.6) nM a.u.−2.09 and 2.09 ± 0.10, respectively. The value of parameter a means that at OD650 = 1.0 a.u., which corresponds to the stationary growth phase in our experimental conditions and the bacterial density of 4 × 108 cells mL‒1, the ComX concentration is about 10 nM. In the early exponential growth phase the concentration was about 0.1 nM. Considering parameter b, the value obtained (2.09 ± 0.10) indicates that the ComX concentration increases by the square of bacterial density. This means that with increasing population density, the SM concentration (ComX) increases by the second power, while the amount of ComX per cell increases linearly. This relationship suggests that ComQXPA has an ultra-sensitive encoder module9, where signal molecule production is very sensitive to cell density. The same mathematical relationship can be obtained by assuming that SM production rate per cell corresponds to the product of a specific cell growth rate and a cell density (i.e., population growth rate, for details, see Supplementary Methods, Derivation of ComQXPA communication system model). The dependence of the SM production rate per cell on (a) cell density and (b) the specific cell growth rate can be seen as an alternative way to obtain the ultra-sensitivity of encoders, which is usually achieved by SM dependent positive feedback in many QS systems9. This makes ComX a true indicator of population density, which also encodes information about the cell growth rate.
Fig. 2: The accumulation of signal molecule (SM) during the growth of B. subtilis and fitness cost of SM production.

a The growth curve (OD650) of B. subtilis PS-216 ΔcomP (no signal receptor) producing SM (ComX) that is accumulating in the growth medium of fermenter working in the batch mode of n = 3 biologically independent replicates is presented; b The experimental data n = 3 biologically independent replicates where the data ≥ limit of detection of SM was fitted by Eq. (1); the error bars for SM concertation are standard errors calculated from 8 technical replicates for each biological replicate; c The comparison of growth curves of B. subtilis PS-216 with no signal molecule receptor (ΔcomP) and no signal molecule production and receptor (ΔcomQXP) of n = 3 biologically independent replicates; the OD650 at t = 0 h was corrected with respect to the measured CFU of the inoculum. The slopes of the fitted lines in c correspond to the growth rate divided by log 2; the exponential phase points in the most reliable OD650 region ( >0.1 a.u. and  1.75 h, P  1 values obtained for all fits indicate positive cooperativity (i.e., ultrasensitivity14) in the binding of the transcriptional factor ComA to the srfA promoter43,44. This agrees with the research showing that two molecules of the ComA homodimer cooperatively bind to the two promoter regions located upstream of the RNAP binding sites of srfA13,20,45. The inactivation of the second promoter region decreases the promoter activity of srfA by 100-fold (ref. 13), which underscores the importance of the second binding region, explains n ≥ 2 and the sharp increase in srfA promoter activity with ComX concentration. In addition, we show here that the critical concentration of ComX required to induce quantifiable response (designated here as lower limit response (LLR)) is 0.2–0.5 nM. These results, therefore, suggest that the response per cell depends cooperatively on the ComX concentration and that the cells respond to very low concentrations of ComX.
Fig. 4: Influence of oxygen on the presence of SM (ComX).

a Signal molecule deficient B. subtilis PS-216 (∆comQ, PsrfAA-yfp) was incubated in the presence of SM for 4 h and the maximum normalized response was determined from the activity of the srfA promoter. n = 4 biologically independent replicates were performed. Best, concatenated fit to the model in Eq. (2) is presented together with 95% confidence level. b the logistic fit to one of the three growth curves (n = 3, biologically independent replicates) measured as OD650 of the culture B. subtilis PS 216 (srfA-lacZ) producing signal molecule, SM that accumulated in the growth medium of batch system is shown. The response per cell data, obtained as the β-galactosidase activity of srfA promoter of B. subtilis PS-216 (srfA-lacZ) was fitted by Eq. (3), (R2  > 0.99). The time interval at which SM concentration is high enough to cause the measurable response, i.e., lower limit of response (LLR) as predicted from data in experiments in (a) is given as dashed window in (b). One of the five qualitatively and quantitatively similar experimental results is presented. Error bars represent SD of 8 technical replicates. For fits of additional replicates and data variability refer to Supplementary Tables 3 and 4.

Full size image

Fully functional ComQXPA communication system does not require a positive feedback loop–the validation of the ComQXPA communication system
The response curve in Fig. 4a is a function of the SM concentration only. In the more natural setting (i.e., during growth) the cells encounter growth-dependent changes in SM concentrations as well as changes in bacterial density and growth rate over time. We have therefore asked whether the response curve based on the modeling and results presented in Figs. 2b and  4a could fit the response data in the SM producing and responding strain exposed to changes in these three parameters.
We cultivated the SM producing and responding PS-216 strain carrying the response reporter (PsrfA-lacZ) in a large volume bioreactor system (Fig. 4b). This allowed sterile sampling of spent medium and cells (for response quantification) at several time points, without affecting growth conditions. Immediately after the inoculation of the fresh medium by overnight culture the β-galactosidase activity of PS-216 (PsrfA-lacZ) was high. We assumed that this was a consequence of the accumulation of the expressed PsrfA reporter (β-galactosidase, RM) during the overnight growth. As a consequence of the dilution of the intracellular β-galactosidase (RM) due to cell division, the activity of the β-galactosidase decreased sharply after 2 h incubation (Fig. 4b). Simultaneously, as predicted by (Eq. 1), the concentration of SM (ComX) in the medium was increased exponentially during growth, and soon reached a critical concentration to activate the srfA promoter. In particular, as elucidated by the fits of (Eq. 2) to the data in Fig. 4a, the lower limit of the response (LLR) is reached shortly before upturn of the cell response curve in Fig. 4b. At this point the culture is in exponential growth phase at the cell density of 3 to 8 × 107 cells mL−1. The steep slope of the response curve indicates that the rate by which the response molecule (RM) is synthesized now exceeds the dilution due to the cell division rate. From now on, the response per cell correlates approximately linearly with OD650, suggesting a strong coupling to cell growth. Taking these facts into account and considering that the response molecule (RM) concentration is sensitive to the concentration of the signal, SM, (Eq. 2) and that SM can be expressed in terms of cell density (Eq. 1), the concentration of a response molecule per cell, RM(t)/N(t), can be analytically described (see also Supplementary Methods, Derivation of ComQXPA communication system model) as:

$$frac{{RM(t)}}{{N(t)}} = frac{{RM0}}{{N(t)}} + frac{{RM1(t)}}{{N(t)}}$$
(3)

where RM0/N(t) is the response per cell of overnight culture, i.e., the overnight accumulated β-galactosidase. The second term, RM1(t)/N(t) accounts for the synthesis of the β-galactosidase after inoculation of a fresh medium and comprises the parameters describing the sensitivity of the response to a signal molecule, Wmax, Km, n (Eq. 2), the signal production, a (Eq. 1), cell density, N(t) and proportionality constant, k that gives the magnitude of the response per cell when the potential to respond to the signal is maximally fulfilled (i.e., at Wmax) and the specific growth rate is 1 h−1. The definition of RM1(t) is given in Supplementary Methods (eq S11). Note that for the derivation of Eq. (3) we assumed no degradation of SM occurs, as our experiments suggest SM was stable under the conditions studied (Supplementary Figs. 6a and 7, see also Supplementary Methods, Derivation of ComQXPA communication). All the parameters in (Eq. 3), except k in RM1(t) and RM0, were taken as constants obtained in the independent experiments by fits of (Eq. 1) and (Eq. 2). With k and RM0, as the only fitting parameters, we applied the mathematical model in (Eq. 3) (for details of the model equation see Supplementary Methods, Derivation of ComQXPA communication system model) to fit the experimental cell response data (Fig. 4b). The successful fit (see Supplementary Table 4 for details) indicates that the relationship assumed among cell density, cell growth, signal concentration and response in (Eq. 3) is valid and yields (760 ± 120) M.U. for k and (5.5 ± 1.5) M.U. a.u. for RM0. Again, we did not need to incorporate the SM feedback loop into our model, which is consistent with published results suggesting that this communication system lacks a feedback loop10,11,12,13.
The ComQXPA dependent signaling and response at the cellular level
So far, we have focused on the population averages, which is a traditional approach in studies on microbial communication systems17,46. We here report results on communication dynamics of B. subtilis at the single cell level using fluorescence-based molecular tools. This approach provides the means to track temporal changes in expression of genes involved in signal synthesis (signaling) and in response and thus provides the insight into a phenotypic heterogeneity within the population.
We used the double-labeled fluorescencent strain B. subtilis PS-216 (comQ-yfp, srfA-cfp), in which fluorescent reporters were fused to the comQ and srfA promoters. The two genes code for the ComX signal-processing protein and the communication-activated operon, respectively. Since comQ and comX share the same promoter and their genetic sequences often overlap15 expression level of comQ corresponds to the expression level of comX. The fluorescence of individual cells was observed under the microscope in different growth phases and quantitatively analyzed (Fig. 5).
Fig. 5: Single cell quantitative fluorescence microscopy of B. subtilis PS-216 (comQ-yfp, srfA-cfp).

The fluorescence microscopy images were taken periodically during incubation of B. subtilis PS-216 (comQ-yfp, srfA-cfp) in a batch fermenter by the YFP filter (a), CFP filter (b) or DIC (c). The example YFP and CFP images, taken after 3 h of incubation were pseudo-colored. The scale bar represents 10 µm. n = 3 biologically independent experiments were performed. d % of population of cells that are hyper-expressing comQ-yfp or srfA-cfp is depicted. Gene expression level determined by single cell fluorescence microscopy was measured as Na-fluoresceinate standard normalized mean fluorescence intensity per cells expressing comQ-yfp (e) and srfA-cfp (f); ON is the overnight culture. One of the three qualitatively similar cell distributions is shown in (g) and (h) for comQ-yfp and srfA-cfp, respectively; areas under the curves are the same in all time points. For additional replicate see Supplementary Fig. 8.

Full size image

The observed expression pattern for the signaling gene (comQ-yfp) follows lognormal distribution (Fig. 5g). A small number of outliers in the comQ expression (on average, 10x brighter than the majority) were easily detected in the qualitative image analysis (Fig. 5a). These hyperproducers were not present in the overnight culture and began to occur during exponential growth, after 1-hour incubation in fresh medium, (Fig. 5d). In general, hyperproducers accounted for about 0.1–1% of the population, and their frequency increased during the first 6 hours. These data suggest that bulk of the ComX is nevertheless produced by the majority of the population as expected for the true QS communication systems. The contribution of the srfA hyperproducers to the total surfactin production is even less pronounced since their occurrence did not exceed 0.1% of the population (Fig. 5b, d).
The most heterogeneous expression of the communication signaling gene (comQ-yfp) was observed in overnight culture, immediately after its transfer to the fresh medium (Fig. 5g, Supplementary Fig. 8a), but hyperproducers where not detectable at this time (Fig. 5d). Once the cells begun to divide, the distribution shifted to lower fluorescence intensities with a simultaneous decrease in heterogeneity, but from 3 to 4 h onwards single cell fluorescence gradually increased, along with an increase in population heterogeneity (Fig. 5e, g, Supplementary Fig. 8a). This suggests that the expression rate is now higher than the division rate (i.e., the production overpowers the dilution due to cell division). A similar pattern was observed in the communication response (srfA-cfp) (Fig. 5f, h, Supplementary Fig. 8b) with two major differences. The minimum level in comQ-Yfp fluorescence was reached 1 h later than srfA-Cfp fluorescence, which suggests the expression of ComX is in first hours low compared to the srfA expression. Nevertheless, the entire cell population, with the exception of hyperproducers, which represented only a fraction of the comQ/srfA expressing cells, followed unimodal lognormal distribution expression pattern. This suggests that the ComQXPA communication phenomenon in B. subtilis, at least under the conditions in our experiment, is not restricted to individuals and can be studied at the population level, i.e., the averages represent well the population.
The comQ-yfp and srfA-cfp expression co-localization analysis (Supplementary Table 5) reveals the correlation coefficient of about 0.5, which is significantly (P = 0.01) higher than in the overnight culture. The presence of the correlation suggests that on average, cells that produce the signal more intensively, also respond to signal more intensively, supporting the idea of self-sensing47. However, the correlation coefficient strength was only moderate, suggesting that sensing of the external signal (sensing-of others) still works as expected for a typical QS system.
The induced response in ComQXPA communication system is graded and almost switch-like
The perfect QS system does not produce a response until the threshold bacterial density is reached and then immediately switches to a full response. This minimum to maximum transition may be either a perfect switch or a graded induction. By combining the information from Figs. 2b and 4a in the form of eq S9 results in the normalized ComQXPA response curve as a function of bacterial density (Fig. 6a) that resembles a graded switch like induction (compare to Fig. 1a). The perfect switch like communication system is unrealistic, because it requires that all the cells are perfectly synchronized and immediately switch to maximal response, leaving no time for the adaptation to the signal stimuli. It is reasonable to expect that for the true quorum sensing system (QS) most of the response has to occur within the same generation of dividing cells (ngen  More

1.
USDAFS. Press Release: Survey Finds 18 Million Trees Died in California in 2018. https://www.fs.usda.gov/Internet/FSE_DOCUMENTS/FSEPRD609321.pdf (USDAFS, 2019).
2.
Griffin, D. & Anchukaitis, K. J. How unusual is the 2012-2014 California drought? Geophys. Res. Lett. 41, 9017–9023 (2014).
ADS  Article  Google Scholar 

3.
Robeson, S. M. Revisiting the recent California drought as an extreme value. Geophys. Res. Lett. 42, 6771–6779 (2015).
ADS  Article  Google Scholar 

4.
Asner, G. P. et al. Progressive forest canopy water loss during the 2012-2015 California drought. Proc. Natl Acad. Sci. USA 113, E249–E255 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
Brodrick, P. G. & Asner, G. P. Remotely sensed predictors of conifer tree mortality during severe drought. Environ. Res. Lett. 12, 115013 (2017).
ADS  Article  Google Scholar 

6.
Fettig, C. J. in Managing Sierra Nevada Forests. PSW-GTR-237 Ch. 2 (USDA Forest Service, 2012).

7.
Kolb, T. E. et al. Observed and anticipated impacts of drought on forest insects and diseases in the United States. For. Ecol. Manag. 380, 321–334 (2016).
Article  Google Scholar 

8.
Waring, R. H. & Pitman, G. B. Modifying lodgepole pine stands to change susceptibility to mountain pine beetle attack. Ecology 66, 889–897 (1985).
Article  Google Scholar 

9.
Restaino, C. et al. Forest structure and climate mediate drought-induced tree mortality in forests of the Sierra Nevada, USA. Ecol. Appl. 0, e01902 (2019).
Article  Google Scholar 

10.
USDAFS. Press Release: Record 129 million dead trees in California. https://www.fs.usda.gov/Internet/FSE_DOCUMENTS/fseprd566303.pdf (USDAFS, 2017).

11.
Young, D. J. N. et al. Long-term climate and competition explain forest mortality patterns under extreme drought. Ecol. Lett. 20, 78–86 (2017).
PubMed  Article  PubMed Central  Google Scholar 

12.
Raffa, K. F. et al. Cross-scale drivers of natural disturbances prone to anthropogenic amplification: The dynamics of bark beetle eruptions. BioScience 58, 501–517 (2008).
Article  Google Scholar 

13.
Boone, C. K., Aukema, B. H., Bohlmann, J., Carroll, A. L. & Raffa, K. F. Efficacy of tree defense physiology varies with bark beetle population density: a basis for positive feedback in eruptive species. Can. J. Res. 41, 1174–1188 (2011).
Article  Google Scholar 

14.
Fettig, C. J., Mortenson, L. A., Bulaon, B. M. & Foulk, P. B. Tree mortality following drought in the central and southern Sierra Nevada, California, U.S. For. Ecol. Manag. 432, 164–178 (2019).
Article  Google Scholar 

15.
Stephenson, N. L., Das, A. J., Ampersee, N. J. & Bulaon, B. M. Which trees die during drought? The key role of insect host-tree selection. J. Ecol. 75, 2383–2401 (2019).
Article  Google Scholar 

16.
Senf, C., Campbell, E. M., Pflugmacher, D., Wulder, M. A. & Hostert, P. A multi-scale analysis of western spruce budworm outbreak dynamics. Landsc. Ecol. 32, 501–514 (2017).
Article  Google Scholar 

17.
Seidl, R. et al. Small beetle, large-scale drivers: how regional and landscape factors affect outbreaks of the European spruce bark beetle. J. Appl Ecol. 53, 530–540 (2016).
Article  Google Scholar 

18.
Fettig, C. J. in Insects and Diseases of Mediterranean Forest Systems (eds Lieutier, F. & Paine, T. D.) 499–528 (Springer International Publishing, 2016).

19.
Raffa, K. F. & Berryman, A. A. The role of host plant resistance in the colonization behavior and ecology of bark beetles (Coleoptera: Scolytidae). Ecol. Monogr. 53, 27–49 (1983).
Article  Google Scholar 

20.
Logan, J. A., White, P., Bentz, B. J. & Powell, J. A. Model analysis of spatial patterns in mountain pine beetle outbreaks. Theor. Popul. Biol. 53, 236–255 (1998).
CAS  PubMed  MATH  Article  PubMed Central  Google Scholar 

21.
Wallin, K. F. & Raffa, K. F. Feedback between individual host selection behavior and population dynamics in an eruptive herbivore. Ecol. Monogr. 74, 101–116 (2004).
Article  Google Scholar 

22.
Franceschi, V. R., Krokene, P., Christiansen, E. & Krekling, T. Anatomical and chemical defenses of conifer bark against bark beetles and other pests. N. Phytol. 167, 353–376 (2005).
CAS  Article  Google Scholar 

23.
Raffa, K. F., Grégoire, J.-C. & Staffan Lindgren, B. Natural History and Ecology of Bark Beetles 1–40 (Elsevier, 2015).

24.
Bentz, B. J. et al. Climate change and bark beetles of the western United States and Canada: direct and indirect effects. BioScience 60, 602–613 (2010).
Article  Google Scholar 

25.
DeRose, R. J. & Long, J. N. Drought-driven disturbance history characterizes a southern Rocky Mountain subalpine forest. Can. J. Res. 42, 1649–1660 (2012).
Article  Google Scholar 

26.
Hart, S. J., Veblen, T. T., Schneider, D. & Molotch, N. P. Summer and winter drought drive the initiation and spread of spruce beetle outbreak. Ecology 98, 2698–2707 (2017).
PubMed  Article  PubMed Central  Google Scholar 

27.
Netherer, S., Panassiti, B., Pennerstorfer, J. & Matthews, B. Acute drought Is an important driver of bark beetle infestation in Austrian Norway spruce stands. Front. For. Glob. Change 2, 39 (2019).

28.
Kaiser, K. E., McGlynn, B. L. & Emanuel, R. E. Ecohydrology of an outbreak: mountain pine beetle impacts trees in drier landscape positions first. Ecohydrology 6, 444–454 (2013).
Article  Google Scholar 

29.
Marini, L. et al. Climate drivers of bark beetle outbreak dynamics in Norway spruce forests. Ecography 40, 1426–1435 (2017).
Article  Google Scholar 

30.
Sambaraju, K. R., Carroll, A. L. & Aukema, B. H. Multiyear weather anomalies associated with range shifts by the mountain pine beetle preceding large epidemics. For. Ecol. Manag. 438, 86–95 (2019).
Article  Google Scholar 

31.
Hayes, C. J., Fettig, C. J. & Merrill, L. D. Evaluation of multiple funnel traps and stand characteristics for estimating western pine beetle-caused tree mortality. J. Econ. Entomol. 102, 2170–2182 (2009).
PubMed  Article  PubMed Central  Google Scholar 

32.
Thistle, H. W. et al. Surrogate pheromone plumes in three forest trunk spaces: composite statistics and case studies. For. Sci. 50, 610–625 (2004).

33.
Miller, J. M. & Keen, F. P. Biology and Control of the Western Pine Beetle: A Summary of The First Fifty Years of Research (US Department of Agriculture, 1960).

34.
Chubaty, A. M., Roitberg, B. D. & Li, C. A dynamic host selection model for mountain pine beetle, Dendroctonus ponderosae Hopkins. Ecol. Model. 220, 1241–1250 (2009).
Article  Google Scholar 

35.
Graf, M., Reid, M. L., Aukema, B. H. & Lindgren, B. S. Association of tree diameter with body size and lipid content of mountain pine beetles. Can. Entomol. 144, 467–477 (2012).
Article  Google Scholar 

36.
Geiszler, D. R. & Gara, R. I. in Theory and Practice of Mountain Pine Beetle Management in Lodgepole Pine Forests: Symposium Proceedings (eds Berryman, A. A., Amman, G. D. & Stark, R. W.) (1978).

37.
Klein, W. H., Parker, D. L. & Jensen, C. E. Attack, emergence, and stand depletion trends of the mountain pine beetle in a lodgepole pine stand during an outbreak. Environ. Entomol. 7, 732–737 (1978).
Article  Google Scholar 

38.
Mitchell, R. G. & Preisler, H. K. Analysis of spatial patterns of lodgepole pine attacked by outbreak populations of the mountain pine beetle. For. Sci. 37, 1390–1408 (1991).
Google Scholar 

39.
Preisler, H. K. Modelling spatial patterns of trees attacked by bark-beetles. Appl. Stat. 42, 501 (1993).
MATH  Article  Google Scholar 

40.
Jactel, H. & Brockerhoff, E. G. Tree diversity reduces herbivory by forest insects. Ecol. Lett. 10, 835–848 (2007).
PubMed  Article  Google Scholar 

41.
Faccoli, M. & Bernardinelli, I. Composition and elevation of spruce forests affect susceptibility to bark beetle attacks: implications for forest management. Forests 5, 88–102 (2014).
Article  Google Scholar 

42.
Berryman, A. A. in Bark Beetles in North American Conifers: A System for the Study of Evolutionary Biology 264–314 (University of Texas Press, 1982).

43.
Fettig, C. J. et al. The effectiveness of vegetation management practices for prevention and control of bark beetle infestations in coniferous forests of the western and southern United States. For. Ecol. Manag. 238, 24–53 (2007).
Article  Google Scholar 

44.
Moeck, H. A., Wood, D. L. & Lindahl, K. Q. Host selection behavior of bark beetles (Coleoptera: Scolytidae) attacking Pinus ponderosa, with special emphasis on the western pine beetle, Dendroctonus brevicomis. J. Chem. Ecol. 7, 49–83 (1981).
CAS  PubMed  Article  Google Scholar 

45.
Evenden, M. L., Whitehouse, C. M. & Sykes, J. Factors influencing flight capacity of the mountain pine beetle (Coleoptera: Curculionidae: Scolytinae). Environ. Entomol. 43, 187–196 (2014).
CAS  PubMed  Article  Google Scholar 

46.
Raffa, K. F. & Berryman, A. A. Accumulation of monoterpenes and associated volatiles following inoculation of grand fir with a fungus transmitted by the fir engraver, Scolytus ventralis (Coleoptera: Scolytidae). Can. Entomol. 114, 797–810 (1982).
CAS  Article  Google Scholar 

47.
Anderegg, W. R. L. et al. Tree mortality from drought, insects, and their interactions in a changing climate. N. Phytol. 208, 674–683 (2015).
Article  Google Scholar 

48.
Kane, V. R. et al. Assessing fire effects on forest spatial structure using a fusion of Landsat and airborne LiDAR data in Yosemite National Park. Remote Sens. Environ. 151, 89–101 (2014).
ADS  Article  Google Scholar 

49.
Larson, A. J. & Churchill, D. Tree spatial patterns in fire-frequent forests of western North America, including mechanisms of pattern formation and implications for designing fuel reduction and restoration treatments. For. Ecol. Manag. 267, 74–92 (2012).
Article  Google Scholar 

50.
Morris, J. L. et al. Managing bark beetle impacts on ecosystems and society: Priority questions to motivate future research. J. Appl. Ecol. 54, 750–760 (2017).
Article  Google Scholar 

51.
Shiklomanov, A. N. et al. Enhancing global change experiments through integration of remote-sensing techniques. Front. Ecol. Environ. 17, 215–224 (2019).

52.
Jeronimo, S. M. A. et al. Forest structure and pattern vary by climate and landform across active-fire landscapes in the montane Sierra Nevada. For. Ecol. Manag. 437, 70–86 (2019).
Article  Google Scholar 

53.
Roussel, J.-R., Auty, D., De Boissieu, F. & Meador, A. S. lidR: Airborne LiDAR Data Manipulation and Visualization for Forestry Applications (2019).

54.
McDowell, N. et al. Mechanisms of plant survival and mortality during drought: why do some plants survive while others succumb to drought? N. Phytol. 178, 719–739 (2008).
Article  Google Scholar 

55.
Seybold, S. J. et al. Management of western North American bark beetles with semiochemicals. Annu. Rev. Entomol. 63, 407–432 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

56.
Fettig, C. J., McKelvey, S. R. & Huber, D. P. W. Nonhost angiosperm volatiles and verbenone disrupt response of western pine beetle, Dendroctonus brevicomis (Coleoptera: Scolytidae), to attractant-baited traps. J. Econ. Entomol. 98, 2041–2048 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

57.
Fettig, C. J., Dabney, C. P., McKelvey, S. R. & Huber, D. P. W. Nonhost angiosperm volatiles and verbenone protect individual ponderosa pines from attack by western pine beetle and red turpentine beetle (Coleoptera: Curculionidae, Scolytinae). West J. Appl. 23, 40–45 (2008).
Google Scholar 

58.
Fettig, C. J. et al. Efficacy of ‘Verbenone Plus’ for protecting ponderosa pine trees and stands from Dendroctonus brevicomis (Coleoptera: Curculionidae) attack in British Columbia and California. J. Econ. Entomol. 105, 1668–1680 (2012).
PubMed  Article  PubMed Central  Google Scholar 

59.
Oliver, W. W. Is self-thinning in ponderosa pine ruled by Dendroctonus bark beetles? In Forest Health Through Silviculture: Proceedings of the 1995 National Silviculture Workshop 6 (1995).

60.
Fettig, C. & McKelvey, S. Resiliency of an interior ponderosa pine forest to bark beetle infestations following fuel-reduction and forest-restoration treatments. Forests 5, 153–176 (2014).
Article  Google Scholar 

61.
Fettig, C. J. & Hilszczański, J. Bark Beetles 555–584. https://doi.org/10.1016/B978-0-12-417156-5.00014-9 (Elsevier, 2015).

62.
Chesson, P. Mechanisms of maintenance of species diversity. Annu. Rev. Ecol. Syst. 31, 343–366 (2000).
Article  Google Scholar 

63.
Fricker, G. A. et al. More than climate? Predictors of tree canopy height vary with scale in complex terrain, Sierra Nevada, CA (USA). For. Ecol. Manag. 434, 142–153 (2019).
Article  Google Scholar 

64.
Ma, S., Concilio, A., Oakley, B., North, M. & Chen, J. Spatial variability in microclimate in a mixed-conifer forest before and after thinning and burning treatments. For. Ecol. Manag. 259, 904–915 (2010).
Article  Google Scholar 

65.
Stovall, A. E. L., Shugart, H. & Yang, X. Tree height explains mortality risk during an intense drought. Nat. Commun. 10, 1–6 (2019).
Article  CAS  Google Scholar 

66.
Stephenson, N. L. & Das, A. J. Height-related changes in forest composition explain increasing tree mortality with height during an extreme drought. Nat. Commun. 11, 3402 (2020).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

67.
Stovall, A. E. L., Shugart, H. H. & Yang, X. Reply to ‘Height-related changes in forest composition explain increasing tree mortality with height during an extreme drought’. Nat. Commun. 11, 3401 (2020).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

68.
Person, H. L. Tree selection by the western pine beetle. J. For. 26, 564–578 (1928).
Google Scholar 

69.
Person, H. L. Theory in explanation of the selection of certain trees by the western pine beetle. J. For. 29, 696–699 (1931).
CAS  Google Scholar 

70.
Pile, L. S., Meyer, M. D., Rojas, R., Roe, O. & Smith, M. T. Drought impacts and compounding mortality on forest trees in the southern Sierra Nevada. Forests 10, 237 (2019).
Article  Google Scholar 

71.
Frey, J., Kovach, K., Stemmler, S. & Koch, B. UAV photogrammetry of forests as a vulnerable process. A sensitivity analysis for a structure from motion RGB-image pipeline. Remote Sens. 10, 912 (2018).
ADS  Article  Google Scholar 

72.
James, M. R. & Robson, S. Mitigating systematic error in topographic models derived from UAV and ground-based image networks. Earth Surf. Process. Landf. 39, 1413–1420 (2014).
ADS  Article  Google Scholar 

73.
Gray, P. C. et al. A convolutional neural network for detecting sea turtles in drone imagery. Methods Ecol. Evol. 10, 345–355 (2019).
Article  Google Scholar 

74.
Millar, C. I., Stephenson, N. L. & Stephens, S. L. Climate change and forests of the future: Managing in the face of uncertainty. Ecol. Appl. 17, 2145–2151 (2007).
PubMed  Article  PubMed Central  Google Scholar 

75.
Vose, J. M. et al. in Impacts, Risks, and Adaptation in the United States: The Fourth National Climate Assessment Vol II (eds Reidmiller, D. R., et al.) 232–267. https://nca2018.globalchange.gov/chapter/6/https://doi.org/10.7930/NCA4.2018.CH6 (2018).

76.
Bedard, W. D. et al. Western pine beetle: field response to its sex pheromone and a synergistic host terpene, myrcene. Science 164, 1284–1285 (1969).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

77.
Byers, J. A. & Wood, D. L. Interspecific inhibition of the response of the bark beetles, Dendroctonus brevicomis and Ips paraconfusus, to their pheromones in the field. J. Chem. Ecol. 6, 149–164 (1980).
CAS  Article  Google Scholar 

78.
Shepherd, W. P., Huber, D. P. W., Seybold, S. J. & Fettig, C. J. Antennal responses of the western pine beetle, Dendroctonus brevicomis (Coleoptera: Curculionidae), to stem volatiles of its primary host, Pinus ponderosa, and nine sympatric nonhost angiosperms and conifers. Chemoecology 17, 209–221 (2007).
CAS  Article  Google Scholar 

79.
DJI. Zenmuse X3 – Creativity Unleashed. DJI Official https://www.dji.com/zenmuse-x3/info (2015).

80.
Micasense. MicaSense. https://support.micasense.com/hc/en-us/articles/215261448-RedEdge-User-Manual-PDF-Download- (2015).

81.
DJI. DJI – The World Leader in Camera Drones/Quadcopters for Aerial Photography. DJI Official https://www.dji.com/matrice100/info (2015).

82.
Wyngaard, J. et al. Emergent challenges for science sUAS data management: fairness through community engagement and best practices development. Remote Sens. 11, 1797 (2019).
ADS  Article  Google Scholar 

83.
Rouse, W., Haas, R. H., Deering, W. & Schell, J. A. Monitoring the Vernal Advancement and Retrogradation (Green Wave Effect) of Natural Vegetation (Remote Sensing Center, Texas A&M Univ., 1973).

84.
DronesMadeEasy. Map Pilot for DJI on iOS. App Store https://itunes.apple.com/us/app/map-pilot-for-dji/id1014765000?mt=8 (2018).

85.
Farr, T. G. et al. The shuttle radar topography mission. Rev. Geophys. 45, RG2004 (2007).

86.
Zhang, W. et al. An easy-to-use airborne LiDAR data filtering method based on cloth simulation. Remote Sens. 8, 501 (2016).
ADS  Article  Google Scholar 

87.
Hijmans, R. J. et al. Raster: Geographic Data Analysis and Modeling (2019).

88.
Gitelson, A. & Merzlyak, M. N. Spectral reflectance changes associated with autumn senescence of Aesculus hippocastanum L. and Acer platanoides L. leaves. Spectral features and relation to chlorophyll estimation. J. Plant Physiol. 143, 286–292 (1994).
CAS  Article  Google Scholar 

89.
Coops, N. C., Johnson, M., Wulder, M. A. & White, J. C. Assessment of QuickBird high spatial resolution imagery to detect red attack damage due to mountain pine beetle infestation. Remote Sens. Environ. 103, 67–80 (2006).
ADS  Article  Google Scholar 

90.
Clevers, J. G. P. W. & Gitelson, A. A. Remote estimation of crop and grass chlorophyll and nitrogen content using red-edge bands on Sentinel-2 and -3. Int. J. Appl. Earth Observ. Geoinf. 23, 344–351 (2013).
ADS  Article  Google Scholar 

91.
Li, W., Guo, Q., Jakubowski, M. K. & Kelly, M. A new method for segmenting individual trees from the LiDAR point cloud. Photogramm. Eng. Remote Sens. 78, 75–84 (2012).
Article  Google Scholar 

92.
Jakubowski, M. K., Li, W., Guo, Q. & Kelly, M. Delineating individual trees from LiDAR data: a comparison of vector- and raster-based segmentation approaches. Remote Sens. 5, 4163–4186 (2013).
ADS  Article  Google Scholar 

93.
Shin, P., Sankey, T., Moore, M. & Thode, A. Evaluating unmanned aerial vehicle images for estimating forest canopy fuels in a ponderosa pine stand. Remote Sens. 10, 1266 (2018).
ADS  Article  Google Scholar 

94.
Roussel, J.-R. lidRplugins: Extra Functions and Algorithms for lidR Package (2019).

95.
Eysn, L. et al. A benchmark of LiDAR-based single tree detection methods using heterogeneous forest data from the alpine space. Forests 6, 1721–1747 (2015).
Article  Google Scholar 

96.
Vega, C. et al. PTrees: a point-based approach to forest tree extraction from LiDAR data. Int. J. Appl. Earth Observ. Geoinf. 33, 98–108 (2014).
ADS  Article  Google Scholar 

97.
Plowright, A. ForestTools: Analyzing Remotely Sensed Forest Data (2018).

98.
Pau, G., Fuchs, F., Sklyar, O., Boutros, M. & Huber, W. EBImage: An R package for image processing with applications to cellular phenotypes. Bioinformatics 26, 979–981 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

99.
Meyer, F. & Beucher, S. Morphological segmentation. J. Vis. Commun. Image Represent. 1, 21–46 (1990).
Article  Google Scholar 

100.
Hunziker, P. Velox: Fast Raster Manipulation and Extraction (2017).

101.
Kuhn, M. Building predictive models in R using the caret package. J. Stat. Softw. 28, 1–26 (2008).
Article  Google Scholar 

102.
Wang, Y. et al. Is field-measured tree height as reliable as believed – a comparison study of tree height estimates from field measurement, airborne laser scanning and terrestrial laser scanning in a boreal forest. ISPRS J. Photogramm. Remote Sens. 147, 132–145 (2019).
ADS  Article  Google Scholar 

103.
Weinstein, B. G., Marconi, S., Bohlman, S., Zare, A. & White, E. Individual tree-crown detection in RGB imagery using semi-supervised deep learning neural networks. Remote Sens. 11, 1309 (2019).
ADS  Article  Google Scholar 

104.
dos Santos, A. A. et al. Assessment of CNN-based methods for individual tree detection on images captured by RGB cameras attached to UAVs. Sensors (Basel) 19, 3595 (2019).

105.
Stephenson, N. Actual evapotranspiration and deficit: biologically meaningful correlates of vegetation distribution across spatial scales. J. Biogeogr. 25, 855–870 (1998).
Article  Google Scholar 

106.
Flint, L. E., Flint, A. L., Thorne, J. H. & Boynton, R. Fine-scale hydrologic modeling for regional landscape applications: The California Basin Characterization Model development and performance. Ecol. Process. 2, 25 (2013).
Article  Google Scholar 

107.
Millar, C. I. et al. Forest mortality in high-elevation whitebark pine (Pinus albicaulis) forests of eastern California, USA: influence of environmental context, bark beetles, climatic water deficit, and warming. Can. J. For. Res. 42, 749–765 (2012).
Article  Google Scholar 

108.
Baldwin, B. G. et al. Species richness and endemism in the native flora of California. Am. J. Bot. 104, 487–501 (2017).
PubMed  Article  PubMed Central  Google Scholar 

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

110.
Hoffman, M. D. & Gelman, A. The no-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J. Mach. Learn. Res. 15, 31 (2014).
MathSciNet  MATH  Google Scholar 

111.
Carpenter, B. et al. Stan: a probabilistic programming language. J. Stat. Softw. 76, 1–32 (2017).
Article  Google Scholar 

112.
Brooks, S. P. & Gelman, A. General methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 7, 434 (1998).
MathSciNet  Google Scholar 

113.
Gabry, J., Simpson, D., Vehtari, A., Betancourt, M. & Gelman, A. Visualization in Bayesian workflow. J. R. Stat. Soc. Ser. A 182, 389–402 (2019).
MathSciNet  Article  Google Scholar 

114.
Koontz, M. J., Latimer, A. M., Mortenson, L. A., Fettig, C. J. & North, M. P. Drone-derived data supporting “Cross-scale interaction of host tree size and climatic water deficit governs bark beetle-induced tree mortality”. https://doi.org/10.17605/OSF.IO/3CWF9 (2020).

115.
Baldwin, B. G. et al. Master spatial file for native California vascular plants used by Baldwin et al. (2017 Amer. J. Bot.), Dryad, Dataset, 2017. https://doi.org/10.6078/D16K5W.

116.
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2018).

117.
Koontz, M. J., Latimer, A. M., Mortenson, L. A., Fettig, C. J. & North, M. P. Local-structure-wpb-severity. https://doi.org/10.17605/OSF.IO/WPK5Z (2019). More

Study design
This study began in 2013, before P. destructans invaded Michigan and Wisconsin, and continued through 2020, after all Michigan and Wisconsin hibernacula were invaded by the pathogen. Throughout this period, we visited 22 hibernacula twice per winter during bat hibernation to quantify bat colony sizes, individual bat roosting temperatures, fungal loads (to which we added a constant 0.0001 before transforming to the log10 scale), and recapture probabilities (Fig. 2). We sampled bats during early hibernation in November, when more than 95% of swarm activity was expected to be finished51, and again during late hibernation in March, when less than 1% of spring emergence activity was expected to have begun51,52. For each sampling event, we counted all bats of all species within the site. We focused our analyses on the little brown bat (Myotis lucifugus). This species suffered severe declines due to WNS24,53, and was the only species abundant enough to provide sufficient sample sizes for our analyses. We divided bat counts by sections within the hibernaculum (i.e., “rooms”) that potentially varied in microclimate, and we used HOBO Pro v2 data loggers to continuously record temperatures every 3 h in these sections (Supplementary Figs. 2 and 7).
After counting all bats, we haphazardly sampled 20–25 individual little brown bats stratified across sections and roughly in proportion to the number of bats in each section (Fig. 2). For each sampled bat, we used a Fluke 62 MAX IR laser thermometer to quantify the temperature of the substrate directly adjacent to the bat (2 years is defined here as the post-invasion period) and that we had sampled during both the pre-invasion period (up to and including year 0, when WNS was first detected) and invasion period (years 1–2 of invasion). From pre-invasion to post-invasion, bat temperature distributions might have changed across these 12 hibernacula via two mechanisms: (1) bat preferences could have remained the same while populations declined in hibernacula and/or microsites with unfavorable temperatures, causing distribution shifts via selection; and/or (2) bat preferences could have changed towards avoidance of hibernacula and/or microsites with high temperature-mediated disease mortality, causing distribution shifts via selection or learning. In both cases, we would expect the average regional bat roosting temperatures to change over time. In contrast, we would only expect the average bat roosting temperature in any given site (or section) to change over time if the site provided a wide range of temperatures to select from, which was not always the case in the hibernacula we studied. Therefore, we quantified how the regional mean early hibernation roosting temperatures for all counted bats changed across all 12 of these hibernacula from pre-invasion to post-invasion (2013–2020, up to four years after the fungus was first detected) using a Gaussian regression, without including a random site effect. We visualized this change using density plots with 1 °C smoothing bin widths.
Though we stratified bat roosting temperature sampling by section, roughly in proportion to the total number of bats present in a section, we only sampled up to 25 little brown bats per site, and we could not sample exactly in proportion to the total bats present. Therefore, calculating the mean temperatures based on only sampled bats might be misleading. Instead, we used the observed temperatures from our sampled bats to estimate the roosting temperatures for bats that were counted but not sampled. In particular, we randomly generated an additional set of temperatures for all unmeasured, counted bats based on the mean and standard deviation of the roosting temperatures of sampled bats in the same section during the same survey. We then used both known (sampled bats) and estimated (counted bats) temperatures in our Gaussian regression to determine how average temperatures changed from pre-invasion to post-invasion. We also confirmed that qualitative results were the same if we only used the temperatures from sampled bats, without estimation. By considering how bats’ early hibernation temperature distributions changed throughout invasion, regardless of infection status, we were quantifying bat roosting preferences within and between hibernacula, rather than behavioral responses to infection within a season (i.e., “behavioral fever”).
Finally, we used the number of bats counted per hibernaculum during early hibernation (November) and late hibernation (March) in all hibernacula that still had bats in post-invasion years to calculate two demographic parameters: over summer population growth rates (reflecting immigration and recruitment to sites between March and November) and over winter mortality rates (population growth rates from November to March). We used a Gaussian linear model to determine whether the log10 population growth rates were correlated with average early hibernation roosting temperatures for sampled bats in each hibernaculum, the demographic season (over summer versus over winter), and the interaction between temperature and demographic season. A significant interaction term would indicate that temperature-correlated bat immigration and recruitment rates were decoupled from temperature-mediated population declines.
All analyses described above were performed in R version 3.5.161 with packages ‘ggplot2’62, ‘lme4’63, ‘loo’64, ‘pROC’65, and ‘R2jags56. Visual inspection of residuals and predictions plots confirmed acceptable model fits for all regression models.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More