Philippa Kaur

1.
Amato P, Brisebois E, Draghi M, Duchane C, Fröhloch-Nowoisky J, Huffman JA, et al. Main biological aerosols, specificities, abundance, and diversity. In: Delort A-M, Amato P, editors. Microbiology of aerosols. 1st ed. John Wiley & Sons, Inc.; Hoboken, New Jersey, USA. 2018. p. 3–11.
2.
Burrows SM, Butler T, Jöckel P, Tost H, Kerkweg A, Pöschl U, et al. Bacteria in the global atmosphere—part 2: modeling of emissions and transport between different ecosystems. Atmos Chem Phys. 2009;9:9281–97.
CAS  Article  Google Scholar 

3.
Wilkinson DM, Koumoutsaris S, Mitchell EAD, Bey I. Modelling the effect of size on the aerial dispersal of microorganisms. J Biogeogr. 2012;39:89–97.
Article  Google Scholar 

4.
Bryan NC, Christner BC, Guzik TG, Granger DJ, Stewart MF. Abundance and survival of microbial aerosols in the troposphere and stratosphere. ISME J. 2019;13:2789–99.
CAS  PubMed  Article  Google Scholar 

5.
Amato P, Joly M, Schaupp C, Attard E, Möhler O, Morris CE, et al. Survival and ice nucleation activity of bacteria as aerosols in a cloud simulation chamber. Atmos Chem Phys. 2015;15:6455–65.
CAS  Article  Google Scholar 

6.
Brown JKM, Hovmøller MS. Aerial dispersal of pathogens on the global and continental scales and its impact on plant disease. Science. 2002;297:537–41.
CAS  PubMed  Article  Google Scholar 

7.
Douwes J, Thorne P, Pearce N, Heederik D. Bioaerosol health effects and exposure assessment: progress and prospects. Ann Occup Hyg. 2003;47:187–200.
CAS  PubMed  Google Scholar 

8.
Fisher MC, Henk DA, Briggs CJ, Brownstein JS, Madoff LC, McCraw SL, et al. Emerging fungal threats to animal, plant and ecosystem health. Nature. 2012;484:186–94.
CAS  PubMed  Article  Google Scholar 

9.
Huffman JA, Prenni AJ, DeMott PJ, Pöhlker C, Mason RH, Robinson NH, et al. High concentrations of biological aerosol particles and ice nuclei during and after rain. Atmos Chem Phys. 2013;13:6151–64.
Article  CAS  Google Scholar 

10.
Möhler O, Demott PJ, Vali G, Levin Z. Microbiology and atmospheric processes: the role of biological particles in cloud physics. Biogeosciences. 2007;1059–71.

11.
Morris CE, Conen F, Alex Huffman J, Phillips V, Pöschl U, Sands DC. Bioprecipitation: a feedback cycle linking Earth history, ecosystem dynamics and land use through biological ice nucleators in the atmosphere. Glob Change Biol. 2014;20:341–51.
Article  Google Scholar 

12.
Sands DC, Langhans VE, Scharen AL, de Smet G. The association between bacteria and rain and possible resultant meteorological implications. J Hung Meteorol Ser. 1982;86:148–52.
Google Scholar 

13.
Burrows SM, Elbert W, Lawrence MG, Pöschl U. Bacteria in the global atmosphere—part 1: review and synthesis of literature data for different ecosystems. Atmos Chem Phys. 2009;9:9263–80.
CAS  Article  Google Scholar 

14.
Lindow SE, Brandl MT. Microbiology of the phyllosphere. Appl Environ Microbiol. 2003;69:1875–83.
CAS  PubMed  PubMed Central  Article  Google Scholar 

15.
Lindemann J, Constantinidou HA, Barchet WR, Upper CD. Plants as sources of airborne bacteria, including ice nucleation-active bacteria. Appl Environ Microbiol. 1982;44:1059–63.

16.
Butterworth J, McCartney HA. The dispersal of bacteria from leaf surfaces by water splash. J Appl Bacteriol. 1991;71:484–96.

17.
Wickman HH. Deposition, adhesion, and release of bioaerosols. Atmos Microb Aerosols. 1994;5:99–165.

18.
Aylor DE, Parlange JY. Ventilation required to entrain small particles from leaves. Plant Physiol. 1975;56:97–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

19.
Lighthart B, Shaffer BT, Marthi B, Ganio LM. Artificial wind-gust liberation of microbial bioaerosols previously deposited on plants. Aerobiologia. 1993;9:189–96.
Article  Google Scholar 

20.
Jones AM, Harrison RM. The effects of meteorological factors on atmospheric bioaerosol concentrations—a review. Sci Total Environ. 2004;326:151–80.
CAS  PubMed  Article  PubMed Central  Google Scholar 

21.
Andreae MO, Merlet P. Emission of trace gases and aerosols from biomass burning. Glob Biogeochem Cycles. 2001;15:955–66.
CAS  Article  Google Scholar 

22.
Yang Y, Chan C, Tao J, Lin M, Engling G, Zhang Z, et al. Observation of elevated fungal tracers due to biomass burning in the Sichuan Basin at Chengdu City, China. Sci Total Environ. 2012;431:68–77.
CAS  PubMed  Article  PubMed Central  Google Scholar 

23.
Rajput P, Anjum MH, Gupta T. One year record of bioaerosols and particles concentration in Indo-Gangetic Plain: implications of biomass burning emissions to high-level of endotoxin exposure. Environ Pollut. 2017;224:98–106.
CAS  PubMed  Article  PubMed Central  Google Scholar 

24.
Kobziar LN, Pingree MRA, Larson H, Dreaden TJ, Green S, Smith JA. Pyroaerobiology: the aerosolization and transport of viable microbial life by wildland fire. Ecosphere. 2018;9:e02507.
Article  Google Scholar 

25.
Mims SA, Mims IIIFM. Fungal spores are transported long distances in smoke from biomass fires. Atmos Environ. 2004;38:651–5.
CAS  Article  Google Scholar 

26.
Bonde MR, Prescott JM, Matsumoto TT, Peterson GL. Possible dissemination of teliospores of Tilletia-indica by the practice of burning wheat stubble. Am Phytopathol Soc. 1987;77:639.

27.
Hu W, Wang Z, Huang S, Ren L, Yue S, Li P, et al. Biological aerosol particles in polluted regions. Curr Pollut Rep. 2020;6:65–89.

28.
Urbanski SP, Reeves MC, Corley RE, Silverstein RP, Hao WMM. Contiguous United States wildland fire emission estimates during 2003-2015 | Rocky Mountain Research Station. Earth Syst Sci Data. 2018;10:2241–74.
Article  Google Scholar 

29.
Kobziar LN, Pingree MRA, Watts AC, Nelson KN, Dreaden TJ, Ridout M. Accessing the Life in Smoke: A New Application of Unmanned Aircraft Systems (UAS) to Sample Wildland Fire Bioaerosol Emissions and Their Environment. 2019;15.

30.
Petters MD, Parsons MT, Prenni AJ, DeMott PJ, Kreidenweis SM, Carrico CM, et al. Ice nuclei emissions from biomass burning. J Geophys Res. 2009;114:D07209.
Google Scholar 

31.
McCluskey CS, DeMott PJ, Prenni AJ, Levin EJT, McMeeking GR, Sullivan AP, et al. Characteristics of atmospheric ice nucleating particles associated with biomass burning in the US: prescribed burns and wildfires. J Geophys Res Atmos. 2014;119:10458–70.
Article  Google Scholar 

32.
Maki LR, Galyan EL, Chang-Chien M-M, Caldwell DR. Ice nucleation induced by Pseudomonas syringae. Appl Environ Microbiol. 1974;28:456–59.

33.
DeMott PJ, Prenni AJ, Liu X, Kreidenweis SM, Petters MD, Twohy CH, et al. Predicting global atmospheric ice nuclei distributions and their impacts on climate. Proc Natl Acad Sci USA. 2010;107:11217–22.
CAS  PubMed  Article  Google Scholar 

34.
Urbanski S. Wildland fire emissions, carbon, and climate: emission factors. For Ecol Manag. 2013;317:51–60.

35.
Schnell RC, Vali G. Atmospheric ice nuclei from decomposing vegetation. Nature. 1972;236:163–5.
Article  Google Scholar 

36.
Beattie GA, Lindow SE. The secret life of foliar bacterial pathogens on leaves. Ann Rev Phytopathol. 1995;33:145–72.

37.
Baldrian P. Forest microbiome: diversity, complexity and dynamics. FEMS Microbiol Rev. 2016;41:fuw040.
Article  CAS  Google Scholar 

38.
Vali G, Christensen M, Fresh RW, Galyan EL, Maki LR, Schnell RC. Biogenic ice nuclei—2. Bacterial Sources. J Atmos Sci. 1976;33:1565–70.
Article  Google Scholar 

39.
FNAI. Guide to the natural communities of Florida: 2010 edition. Tallahassee; Florida Natural Areas Inventory. 2010.

40.
Moore RA, Hanlon R, Powers C, Schmale DG, Christner BC. Scavenging of sub-micron to micron-sized microbial aerosols during simulated rainfall. Atmosphere. 2020;11:80.
Article  Google Scholar 

41.
Vali G. Quantitative evaluation of experimental results an the heterogeneous freezing nucleation of supercooled liquids. J Atmos Sci. 1971;28:402–9.
Article  Google Scholar 

42.
Murray BJ, O’Sullivan D, Atkinson JD, Webb ME. Ice nucleation by particles immersed in supercooled cloud droplets. Chem Soc Rev. 2012;41:6519.
CAS  PubMed  Article  Google Scholar 

43.
DeMott PJ, Prenni AJ. New directions: need for defining the numbers and sources of biological aerosols acting as ice nuclei. Atmos Environ. 2010;44:1944–5.
CAS  Article  Google Scholar 

44.
Turner S, Pryer KM, Miao VPW, Palmer JD. Investigating deep phylogenetic relationships among cyanobacteria and plastids by small subunit rRNA sequence analysis. J. Eukaryot. Microbiol. 1999;46:327–338.

45.
Lane DJ. 16S/23S rRNA sequencing. In: Nucleic acid techniques in bacterial systematics. (Eds Stackebrandt E and Goodfellow M) John Wiley and Sons: New York, NY; 1991.

46.
Reinhardt ED, Keane RE, Brown JK. First order fire effects model: FOFEM 4.0, user’s guide. Missoula; U.S. Department of Agriculture, Forest Service, Intermountain Research Station. 1997.

47.
Kreye JK, Varner JM, Kobziar LN. Long-duration soil heating resulting from forest floor duff smoldering in longleaf pine ecosystems. For Sci. 2020;66:1–13.

48.
Kreye JK, Brewer NW, Morgan P, Varner JM, Smith AMS, Hoffman CM, et al. Fire behavior in masticated fuels: a review. Ecol Manag. 2014;314:193–207.
Article  Google Scholar 

49.
R Core Team. A language and environment for statistical computing. R Core Team; Vienna, Austria. 2018.

50.
Wickham H. ggplot2: elegant graphics for data analysis. Springer-Verlag: New York; 2016.

51.
Yu G, Smith DK, Zhu H, Guan Y, Lam TT-Y. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol. 2017;8:28–36.
Article  Google Scholar 

52.
Radzi bin Abas M, Oros DR, Simoneit BRT. Biomass burning as the main source of organic aerosol particulate matter in Malaysia during haze episodes. Chemosphere. 2004;55:1089–95.
CAS  PubMed  Article  PubMed Central  Google Scholar 

53.
Rajput P, Sarin M, Sharma D, Singh D. Characteristics and emission budget of carbonaceous species from post-harvest agricultural-waste burning in source region of the Indo-Gangetic Plain. Tellus B Chem Phys Meteorol. 2014;66:21026.
Article  Google Scholar 

54.
Reid JS, Koppmann R, Eck TF, Eleuterio DP. A review of biomass burning emissions part II: intensive physical properties of biomass burning particles. Atmos Chem Phys. 2005;5:799–825.
CAS  Article  Google Scholar 

55.
Hungershoefer K, Zeromskiene K, Iinuma Y, Helas G, Trentmann J, Trautmann T, et al. Modelling the optical properties of fresh biomass burning aerosol produced in a smoke chamber: results from the EFEU campaign. Atmos Chem Phys. 2008;8:3427–39.
CAS  Article  Google Scholar 

56.
Robertson KM, Hsieh YP, Bugna GC. Fire environment effects on particulate matter emission factors in Southeastern US pine-grasslands. Atmos Environ. 2014;99:104–11.
CAS  Article  Google Scholar 

57.
Santos-Burgoa C, Rosas I, Yela A. Occurrence of airborne enteric bacteria in Mexico city. Aerobiologia. 1994;10:39–45.
Article  Google Scholar 

58.
Ziemba LD, Beyersdorf AJ, Chen G, Corr CA, Crumeyrolle SN, Diskin G, et al. Airborne observations of bioaerosol over the Southeast United States using a Wideband Integrated Bioaerosol Sensor. J Geophys Res Atmos. 2016;121:8506–24.
CAS  Article  Google Scholar 

59.
Yousefi V, Rama DBK. Monitoring of air for microbial and metal contamination at selected sites in the vicinity of Johannesburg, South Africa. Sci Total Environ. 1992;116:159–67.
CAS  PubMed  Article  Google Scholar 

60.
Bauer H, Kasper-Giebl A, Löflund M, Giebl H, Hitzenberger R, Zibuschka F, et al. The contribution of bacteria and fungal spores to the organic carbon content of cloud water, precipitation and aerosols. Atmos Res. 2002;64:109–19.
CAS  Article  Google Scholar 

61.
Hara K, Zhang D. Bacterial abundance and viability in long-range transported dust. Atmos Environ. 2012;47:20–5.
CAS  Article  Google Scholar 

62.
Vela GR. Survival of Azotobacter in dry soil. Appl Microbiol. 1974;28:77–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

63.
Henriques AO, Moran CP. Structure and assembly of the bacterial endospore coat. Methods. 2000;20:95–110.
CAS  PubMed  Article  PubMed Central  Google Scholar 

64.
Julien B, Kaiser AD, Garza A. Spatial control of cell differentiation in Myxococcus xanthus. Proc Natl Acad Sci USA. 2000;97:9098–103.
CAS  PubMed  Article  PubMed Central  Google Scholar 

65.
Almatroudi A, Tahir S, Hu H, Chowdhury D, Gosbell IB, Jensen SO, et al. Staphylococcus aureus dry-surface biofilms are more resistant to heat treatment than traditional hydrated biofilms. J Hosp Infect. 2018;98:161–7.
CAS  PubMed  Article  PubMed Central  Google Scholar 

66.
Joyce R, Lavender H, Farrar J, Werth JT, Weber CF, D’Andrilli J, et al. Characterization and source identification of biological ice nucleating particles deposited year-round in subtropical precipitation. Appl Environ Microbiol. 2019;85:1–21.

67.
Rosenfeld D, Yu X, Liu G, Xu X, Zhu Y, Yue Z, et al. Glaciation temperatures of convective clouds ingesting desert dust, air pollution and smoke from forest fires. Geophys Res Lett. 2011;38:1–5.

68.
Koren I, Kaufman YJ, Remer LA, Martins JV. Measurement of the effect of Amazon Smoke on inhibition of cloud formation. Science (80-). 2004;303:1342–5.
CAS  Article  Google Scholar 

69.
Lohmann U. A glaciation indirect aerosol effect caused by soot aerosols. Geophys Res Lett. 2002;29:1052.
Article  Google Scholar 

70.
Sassen K, Khvorostyanov VI. Cloud effects from boreal forest fire smoke: evidence for ice nucleation from polarization lidar data and cloud model simulations. Environ Res Lett. 2008;3:025006.
Article  Google Scholar 

71.
Rosenfeld D, Rudich Y, Lahav R. Desert dust suppressing precipitation: a possible desertification feedback loop. Proc Natl Acad Sci USA. 2001;98:5975–80.
CAS  PubMed  Article  Google Scholar 

72.
Andreae MO, Rosenfeld D, Artaxo P, Costa AA, Frank GP, Longo KM, et al. Smoking rain clouds over the Amazon. Science (80-). 2004;303:1337–42.
CAS  Article  Google Scholar 

73.
Rosenfeld D. TRMM observed first direct evidence of smoke from forest fires inhibiting rainfall. Geophys Res Lett. 1999;26:3105–8.
Article  Google Scholar 

74.
Peterson DA, Campbell JR, Hyer EJ, Fromm MD, Kablick GP, Cossuth JH, et al. Wildfire-driven thunderstorms cause a volcano-like stratospheric injection of smoke. npj Clim Atmos Sci. 2018;1:1–8.

75.
Yu P, Toon OB, Bardeen CG, Zhu Y, Rosenlof KH, Portmann RW, et al. Black carbon lofts wildfire smoke high into the stratosphere to form a persistent plume. Science (80-). 2019;365:587–90.
CAS  Article  Google Scholar 

76.
Val Martin M, Kahn R, Tosca M. A Global analysis of wildfire smoke injection heights derived from space-based multi-angle imaging. Remote Sens. 2018;10:1609.
Article  Google Scholar 

77.
Shukla PR, Skea J, Buendia EC, Masson-Delmotte V, Pörtner H-O, Roberts DC, et al. Climate change and land: an IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems. The Intergovernmental Panel on Climate Change. 2019.

78.
Aylor DE. Spread of plant disease on a continental scale: role of aerial dispersal of pathogens. Ecology. 2003;84:1989–97.
Article  Google Scholar 

79.
Nagarajan S, Singh DV. Long-distance dispersion of rust pathogens. Annu Rev Phytopathol. 1990;28:139–53.
CAS  PubMed  Article  Google Scholar 

80.
Abatzoglou JT, Williams AP. Impact of anthropogenic climate change on wildfire across western US forests. Proc Natl Acad Sci USA. 2016;113:11770–5.
CAS  PubMed  Article  Google Scholar 

81.
Sharples JJ, Cary GJ, Fox-Hughes P, Mooney S, Evans JP, Fletcher MS, et al. Natural hazards in Australia: extreme bushfire. Clim Change. 2016;139:85–99.
Article  Google Scholar 

82.
Fernandes K, Verchot L, Baethgen W, Gutierrez-Velez V, Pinedo-Vasquez M, Martius C. Heightened fire probability in Indonesia in non-drought conditions: the effect of increasing temperatures. Environ Res Lett. 2017;12:054002.
Article  Google Scholar 

83.
Le Page Y, Morton D, Hartin C, Bond-Lamberty B, Pereira JMC, Hurtt G, et al. Synergy between land use and climate change increases future fire risk in Amazon forests. Earth Syst Dyn. 2017;8:1237–46.
Article  Google Scholar  More

1.
Veblen, T. T., Hadley, K. S., Reid, M. S. & Rebertus, A. J. The response of subalpine forests to spruce beetle outbreak in Colorado. Ecology 72, 213–231 (1991).
Article  Google Scholar 
2.
Edburg, S. L. et al. Cascading impacts of bark beetle-caused tree mortality on coupled biogeophysical and biogeochemical processes. Front. Ecol. Environ. 10, 416–424 (2012).
Article  Google Scholar 

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

4.
McFarlane, B. L., Stumpf-Allen, R. G. C. & Watson, D. O. Public perceptions of natural disturbance in Canada’s national parks: The case of the mountain pine beetle (Dendroctonus ponderosae Hopkins). Biol. Conserv. 130, 340–348 (2006).
Article  Google Scholar 

5.
Morris, J. L. et al. Bark beetles as agents of change in social-ecological systems. Front. Ecol. Environ. 16, S34–S43 (2018).
Article  Google Scholar 

6.
Beudert, B. et al. Bark beetles increase biodiversity while maintaining drinking water quality. Conserv. Lett. 8, 272–281 (2014).
Article  Google Scholar 

7.
Colorado State Forest Service. Report on the Health of Colorado Forests (Colorado State Forest Service Media, Fort Collins, 2014).
Google Scholar 

8.
Meddens, A. J. & Hicke, J. A. Spatial and temporal patterns of Landsat-based detection of tree mortality caused by a mountain pine beetle outbreak in Colorado, USA. For. Ecol. Manag. 322, 78–88 (2014).
Article  Google Scholar 

9.
Rhoades, P. R., Davis, T. S., Tinkham, W. T. & Hoffman, C. M. Effects of seasonality, forest structure, and understory plant richness on bee community assemblage in a southern Rocky Mountain mixed conifer forest. Ann. Entomol. Soc. Am. 111, 278–284 (2018).
Google Scholar 

10.
Potts, S. G. et al. Global pollinator declines: Trends, impacts and drivers. Trends Ecol. Evol. 25, 345–353 (2010).
Article  Google Scholar 

11.
Harrington, T. B. & Edwards, M. B. Understory vegetation, resource availability, and litterfall responses to pine thinning and woody vegetation control in longleaf pine plantations. Can. J. For. Res. 29, 1055–1064 (1999).
Article  Google Scholar 

12.
Takafumi, H. & Hiura, T. Effects of disturbance history and environmental factors on the diversity and productivity of understory vegetation in a cool-temperate forest in Japan. For. Ecol. Manag. 257, 843–857 (2009).
Article  Google Scholar 

13.
Coleman, T. W. et al. Accuracy of aerial detection surveys for mapping insect and disease disturbances in the United States. For. Ecol. Manag. 430, 321–326 (2018).
Article  Google Scholar 

14.
Holway, J. G. & Ward, R. T. Phenology of alpine plants in northern Colorado. Ecology 46, 73–83 (1965).
Article  Google Scholar 

15.
R Core Team. R: A Language and Environment for Statistical Programming. R Foundation for Statistical Computing, Vienna, Austria. www.R-project.org. (2020).

16.
Brown, M. B. & Forsythe, A. B. Robust tests for the equality of variances. J. Am. Stat. Assoc. 69, 346–367 (1974).
MATH  Google Scholar 

17.
Colwell, R. K. et al. Models and estimators linking individual-based and sample-based rarefaction, extrapolation and comparison of assemblages. J. Plant Ecol. 5, 3–21 (2012).
Article  Google Scholar 

18.
Chao, A. et al. Rarefaction and extrapolation with Hill numbers: A framework for sampling and estimation in species diversity studies. Ecol. Monogr. 84, 45–67 (2014).
Article  Google Scholar 

19.
Hsieh, T.C., Ma, K.H. & Chao, A. iNext: Interpolation and extrapolation for species diversity. R package V 2.0.20 (2020).

20.
Galbraith, S. M., Cane, J. H., Moldenke, A. R. & Rivers, J. W. Wild bee diversity increases with local fire severity in a fire-prone landscape. Ecosphere 10, e02668. https://doi.org/10.1002/ecs2.2668 (2019).
Article  Google Scholar 

21.
Anderson, M. J. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46 (2001).
Google Scholar 

22.
Oksanen, J., Guillaume-Blanchet, F., Friendly, M., Kindt, R., Legendre, P. & McGlinn, D., et al. Community ecology package ‘vegan’. R package V 2.5-6 (2019).

23.
McCabe, L. M., Cobb, N. S. & Butterfield, B. J. Environmental filtering of body size and darker coloration in pollinator communities indicate thermal restrictions on bees, but not flies, at high elevations. PeerJ 7, e7867. https://doi.org/10.7717/peerj.7867 (2019).
Article  PubMed  PubMed Central  Google Scholar 

24.
Oyen, K. J., Giri, S. & Dillon, M. E. Altitudinal variation in bumble bee (Bombus) critical thermal limits. J. Therm. Biol. 59, 52–57 (2016).
Article  Google Scholar 

25.
Woodard, S. H. Bumble bee ecophysiology: Integrating the changing environment and the organism. Curr. Opin. Insect Sci. 22, 101–108 (2017).
Article  Google Scholar 

26.
Carper, A. L. & Bowers, M. D. The Conservation Value of Woody Debris for Cavity-Nesting Bees on Boulder County Open Space (Boulder County Open Space Final Report, Boulder, 2017).
Google Scholar 

27.
Klutsch, J. G. et al. Stand characteristics and downed woody debris accumulations associated with a mountain pine beetle (Dendroctonus ponderosae Hopkins) outbreak in Colorado. For. Ecol. Manag. 258, 641–649 (2009).
Article  Google Scholar 

28.
Fayt, P., Machmer, M. M. & Steeger, C. Regulation of spruce bark beetles by woodpeckers—A literature review. For. Ecol. Manag. 206, 1–14 (2005).
Article  Google Scholar 

29.
Galbraith, S. M., Cane, J. H., Moldenke, A. R. & Rivers, J. W. Salvage logging reduces wild bee diversity, but not abundance, in severely burned mixed-conifer forest. For. Ecol. Manag. 453, 117622 (2019).
Article  Google Scholar 

30.
Angers, V. A., Drapeau, P. & Bergeron, Y. Mineralization rates and factors influencing snag decay in four North American boreal tree species. Can. J. For. Res. 42, 157–166 (2011).
Article  Google Scholar 

31.
Miller-Struttmann, N. E. & Galen, C. High-altitude multi-taskers: Bumble bee food plant use broadens along an altitudinal productivity gradient. Oecologia 176, 1033–1045 (2014).
ADS  Article  Google Scholar 

32.
Burkle, L. A., Simanonok, M. P., Durney, J. S., Myers, J. A. & Belote, R. T. Wildfires influence abundance, diversity, and intraspecific and interspecific trait variation of native bees and flowering plants across burned and unburned landscapes. Front. Ecol. Evol. 7, 252. https://doi.org/10.3389/fevo.2019.00252 (2019).
Article  Google Scholar 

33.
Owen, E. L., Bale, J. S. & Hayward, S. A. L. Can winter-active bumblebees survive the cold? Assessing the cold tolerance of Bombus terrestris audax and the effects of pollen feeding. PLoS ONE 8, e80061. https://doi.org/10.1371/journal.pone.0080061 (2013).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

34.
Rodriguez, A. & Kouki, J. Disturbance-mediated heterogeneity drives pollinator diversity in boreal managed forest ecosystems. Ecol. Appl. 27, 589–602 (2017).
Article  Google Scholar 

35.
Cane, J. H. & Neff, J. L. Predicted fates of ground-nesting bees in soil heated by wildfire: Thermal tolerances of life stages and a survey of nesting depths. Biol. Conserv. 144, 2631–2636 (2011).
Article  Google Scholar 

36.
Odanaka, K., Gibbs, J., Turley, N. E., Isaacs, R. & Brudvig, L. A. Canopy thinning, not agricultural history, determines early responses of wild bees to longleaf pine savanna restoration. Restor. Ecol. 28, 138–146 (2020).
Article  Google Scholar 

37.
Rubene, D., Schroeder, M. & Ranius, T. Diversity patterns of wild bees and wasps in managed boreal forests: Effects of spatial structure, local habitat and surrounding landscape. Biol. Conserv. 184, 201–208 (2015).
Article  Google Scholar 

38.
Mielke, J. L. Rate of deterioration of beetle-killed Engelmann spruce. J. For. 48, 882–888 (1950).
Google Scholar 

39.
Raphael, M. G. & Morrison, M. L. Decay and dynamics of snags in the Sierra Nevada, California. For. Sci. 33, 774–783 (1987).
Google Scholar 

40.
Rhoades, P. R. et al. Sampling technique affects detection of habitat factors influencing wild bee communities. J. Insect Conserv. 21, 703–714 (2017).
Article  Google Scholar 

41.
Westphal, C. et al. Measuring bee diversity in different European habitats and biogeographical regions. Ecol. Monogr. 78, 653–671 (2008).
Article  Google Scholar 

42.
Romme, W. H., Knight, D. H. & Yavitt, J. B. Mountain pine beetle outbreaks in the Rocky Mountains: Regulators of primary productivity?. Am. Nat. 127, 484–494 (1986).
Article  Google Scholar 

43.
Nelson, K. N., Rocca, M. E., Diskin, M., Aoki, C. F. & Romme, W. H. Predictors of bark beetle activity and scale-dependent spatial heterogeneity change during the course of an outbreak in a subalpine forest. Landsc. Ecol. 29, 97–109 (2014).
Article  Google Scholar 

44.
Lozier, J. D., Strange, J. P. & Koch, J. B. Landscape heterogeneity predicts gene flow in a widespread polymorphic bumble bee, Bombus bifarius (Hymenoptera: Apidae). Conserv. Genet. 14, 1099–1110 (2013).
Article  Google Scholar 

45.
Boscolo, D., Tokumoto, P. M., Ferreira, P. A., Ribeiro, J. W. & dos Santos, J. S. Positive responses of flower visiting bees to landscape heterogeneity depend on functional connectivity levels. Perspect. Ecol. Evol. 15, 18–24 (2017).
Google Scholar 

46.
Ründlof, M., Nilsson, H. & Smith, H. G. Interacting effects of farming practice and landscape context on bumble bees. Biol. Conserv. 141, 417–426 (2008).
Article  Google Scholar 

47.
Andersson, G. K., Ekroos, J., Stjernman, M., Ründlof, M. & Smith, H. G. Effects of farming intensity, crop rotation and landscape heterogeneity on field bean pollination. Agric. Ecosyst. Environ. 184, 145–148 (2014).
Article  Google Scholar 

48.
Moreira, E. F., Boscolo, D. & Viana, B. F. Spatial heterogeneity regulates plant–pollinator networks across multiple landscape scales. PLoS ONE 10, e0123628. https://doi.org/10.1371/journal.pone.0123628 (2015).
CAS  Article  PubMed  PubMed Central  Google Scholar  More

1.
Li, Y. H., Hallerman, E. M., Liu, Q. S., Wu, K. M. & Peng, Y. F. The development and status of Bt rice in China. Plant Biotechnol. J. 14, 839–848 (2016).
Article  Google Scholar 
2.
Li, X. et al. Comparison of nutritional quality between chinese indica rice with sck and crylAc genes and its nontransgenic counterpart. J. Food Sci. 72, 420–424 (2007).
Article  Google Scholar 

3.
Liu, Q. S., Hallerman, E., Peng, Y. F. & Li, Y. H. Development of Bt rice and Bt maize in china and their efficacy in target pest control. Int. J. Mol. Sci. 17, 1–15 (2016).
ADS  Google Scholar 

4.
Chen, M., Shelton, A. & Ye, G. Y. Insect-resistant genetically modified rice in china: from research to commercialization. Annu. Rev. Entomol. 56, 81–101 (2011).
CAS  Article  Google Scholar 

5.
Mannakkara, A., Niu, L., Ma, W. H. & Lei, C. L. Zero effect of Bt rice on expression of genes coding for digestion, detoxification and immune responses and developmental performances of brown planthopper Nilaparvata lugens (Stål). J. Insect Physiol. 59, 985–993 (2013).
CAS  Article  Google Scholar 

6.
Tian, J. C. et al. Assessing the effects of Cry1C rice and Cry2A rice to Pseudogonatopus flavifemur, a parasitoid of rice planthoppers. Sci. Rep. 7, 7838. https://doi.org/10.1038/s41598-017-08173-w (2017).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

7.
Tian, J. C. et al. The rice planthopper parasitoid Anagrus nilaparvatae is not at risk when feeding on honeydew derived from Bacillus thuringiensis (Bt) rice. Pest Manag. Sci. 74, 1854–1860 (2018).
CAS  Article  Google Scholar 

8.
Zhang, L., Guo, R., Fang, Z. & Liu, B. Genetically modified rice Bt-Shanyou63 expressing Cry1Ab/c protein does not harm Daphnia magna. Ecotox. Environ. Safe. 132, 196–201 (2016).
CAS  Article  Google Scholar 

9.
Chen, Y. et al. Bt rice expressing Cry1Ab does not stimulate an outbreak of its non-target herbivore, Nilaparvata lugens. Transgen. Res. 21, 279–291 (2012).
CAS  Article  Google Scholar 

10.
Akhtar, Z. R. et al. Impact of six transgenic Bacillus thuringiensis rice lines on four nontarget thrips species attacking rice panicles in the paddy field. Environ. Entomol. 42, 173–180 (2013).
CAS  Article  Google Scholar 

11.
Zhou, X., Cheng, J. A., Yang, H. U. & Lou, Y. G. Effects of transgenic Bt rice on the population development of Nephotettix cincticeps. Chin. J. Rice Sci. 19, 74–78 (2005).
Google Scholar 

12.
Lu, Z. B. et al. Impacts of Bt rice expressing Cry1C or Cry2A protein on the performance of nontarget leafhopper, Nephotettix cincticeps (Hemiptera: Cicadellidae), under laboratory and field conditions. Environ. Entomol. 43, 209–217 (2014).
CAS  Article  Google Scholar 

13.
Hagenbucher, S. et al. Pest trade-offs in technology: reduced damage by caterpillars in Bt cotton benefits aphids. Proc. Biol. Sci. 280, 20130042. https://doi.org/10.1098/rspb.2013.0042 (2013).
CAS  Article  PubMed  PubMed Central  Google Scholar 

14.
Hagenbucher, S., Wackers, F. L. & Romeis, J. Indirect multi-trophic interactions mediated by induced plant resistance: Impact of caterpillar feeding on aphid parasitoids. Biol. Lett. 10, 20130795. https://doi.org/10.1098/rsbl.2013.0795 (2014).
Article  PubMed  PubMed Central  Google Scholar 

15.
Wang, X. Y. et al. Bt rice could provide ecological resistance against nontarget planthoppers. Plant Biotechnol. J. 16, 1748–1755 (2018).
CAS  Article  Google Scholar 

16.
Chen, G. et al. Odor, not performance, dictates Bemisia tabaci’s selection between healthy and virus infected plants. Front. Physiol. 8, 146 (2017).
PubMed  PubMed Central  Google Scholar 

17.
Xu, H. X., He, X. C., Zheng, X. S., Yang, Y. J. & Lu, Z. X. Influence of rice black streaked dwarf virus on the ecological fitness of non-vector planthopper Nilaparvata lugens (Hemiptera: Delphacidae). Insect Sci. 21, 507–514 (2013).
Article  Google Scholar 

18.
Xu, H. X. et al. Effects of transgenic rice infected with SRBSDV on Bt expression and the eological fitness of non-vector brown planthopper Nilaparvata lugens. Sci. Rep. 7, 6328. https://doi.org/10.1038/s41598-017-02218-w (2017).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

19.
Wang, Q. J. et al. Combined influence of Bt rice and rice dwarf virus on biological parameters of a non-target herbivore, Nephotettix cincticeps (Uhler) (Hemiptera: Cicadellidae). PLoS ONE 12, e0181258 (2017).
Article  Google Scholar 

20.
Chen, H. Y. et al. Sequential infection of rice dwarf virus in the internal organs of its insect vector after ingestion of virus. Virus Res. 160, 389–394 (2011).
CAS  Article  Google Scholar 

21.
Honda, K. et al. Retention of rice dwarf virus by descendants of pairs of viruliferous vector insects after rearing for 6 years. Phytopathology 97, 712–716 (2007).
Article  Google Scholar 

22.
Wei, T. Y. & Li, Y. Rice reoviruses in insect vectors. Annu. Rev. Phytopathol. 54, 99–120 (2016).
CAS  Article  Google Scholar 

23.
Zhao, S. S. et al. A viral protein promotes host SAMS1 activity and ethylene production for the benefit of virus infection. elife. 6, e27529. https://doi.org/10.7554/eLife.27529 (2017).
Article  PubMed  PubMed Central  Google Scholar 

24.
Ling, K. C., Tiongco, E. R. & Aguiero, V. M. Rice ragged stunt a new virus disease. Plant. Dis. Rep. 62, 701–705 (1978).
Google Scholar 

25.
Hibino, H. et al. Rice grassy stunt virus: A planthopper-borne circular filament. Phytopathology 75, 894–899 (1985).
Article  Google Scholar 

26.
Shimizu, T. et al. Strong resistance against rice grassy stunt virus is induced in transgenic rice plants expressing double-stranded RNA of the viral genes for nucleocapsid or movement proteins as targets for RNA interference. Phytopathology 103, 513–519 (2013).
Article  Google Scholar 

27.
Dang, C. et al. Does Bt rice pose risks to non-target arthropods? Results of a meta-analysis in China. Plant Biotechnol. J. 15, 1047–1053 (2017).
CAS  Article  Google Scholar 

28.
Ge, L. Q., Wu, J. C., Sun, Y. C., Ouyang, F. & Ge, F. Effects of triazophos on biochemical substances of transgenic Bt rice and its nontarget pest Nilaparvata lugens Stål under elevated CO2. Pest. Biochem. Physiol. 107, 188–199 (2013).
CAS  Article  Google Scholar 

29.
Ge, L. Q., Sun, Y. C., Ouyang, F., Wu, J. C. & Ge, F. The effects of triazophos applied to transgenic Bt rice on the nutritional indexes, Nlvg expression, and population growth of Nilaparvata lugens Stal under elevated CO2. Pest. Biochem. Physiol. 118, 50–57 (2015).
CAS  Article  Google Scholar 

30.
Che, Q. Y., Liu, X. H., Liang, Y. Y., Yang, B. & Ge, F. Effects of transgenic Bt-rice and insecticides on the communitystructure of soil nematodes. J. Plant Protect. 42(5), 724–733 (2015).
Google Scholar 

31.
Yang, Y. J. et al. Impacts of nitrogen fertilizer on major insect pests and their predators in transgenic Bt rice lines T2A–1 and T1C–19. Entomol. Exp. Appl. 160, 281–291 (2016).
CAS  Article  Google Scholar 

32.
Hill, D. S. Agricultural insect pests of temperate regions and their control. Q. Rev. Biol. 63, 344 (1987).
Google Scholar 

33.
Wang, Q. J. et al. Rice dwarf virus infection alters green rice leafhopper host preference and feeding behavior. PLoS ONE 13, e0203364 (2018).
Article  Google Scholar 

34.
Sun, X., Yan, M. J., Zhang, A. J. & Wang, M. Q. Transgenic cry1C gene rough rice line T1C–19 does not change the host preferences of the non-target stored product pest, Rhyzopertha dominica (Fabricius) (Coleoptera: Bostrichidae), and its parasitoid wasp, Anisopteromalus calandrae (Howard) (Hymenoptera: Pteromalidae). Ecotox. Environ. Safe. 120, 449–456 (2015).
CAS  Article  Google Scholar 

35.
Lu, Z. B. et al. Transgenic cry1C or cry2A rice has no adverse impacts on the life-table parameters and population dynamics of the brown planthopper, Nilaparvata lugens (Hemiptera: Delphacidae). Pest Manag. Sci. 71, 937–945 (2015).
ADS  CAS  Article  Google Scholar 

36.
He, X. C. et al. Ecological fitness of non-vector planthopper Sogatella furcifera on rice plants infected with rice black streaked dwarf virus. Rice Sci. 19, 335–338 (2012).
Article  Google Scholar 

37.
Xu, H. X. et al. Effects of SRBSDV-infected rice plants on the fitness of vector and non-vector rice planthoppers. J. Asia-Pac. Entomol. 19, 707–710 (2016).
Article  Google Scholar 

38.
Gómez-Torres, M. L., Nava, D. E. & Parra, J. R. P. Life table of Tamarixia radiata (Hymenoptera: Eulophidae) on Diaphorina citri (Hemiptera: Psyllidae) at different temperatures. J. Econ. Entomol. 105, 338–343 (2012).
Article  Google Scholar 

39.
Tang, W. et al. Development of insect-resistant transgenic indica rice with a synthetic cry1C* gene. Mol. Breed. 18, 1–10 (2006).
CAS  Article  Google Scholar 

40.
Hulting, F. L., Orr, D. B. & Obrycki, J. J. A computer-program for calculation and statistical comparison of intrinsic rates of increase and associated life table parameters. Fla. Entomol. 73, 601–612 (1990).
Article  Google Scholar 

41.
Maia, A. D. N., Luiz, A. J. B. & Campanhola, C. Statistical inference on associated fertility life table parameters using Jackknife technique: computational aspects. J. Econ. Entomol. 93, 511–518 (2000).
Article  Google Scholar  More

Sampling site and regimes
Field work was carried out in the agricultural catchment of Grytelandsbekken (a creek of approx. 2.5 km in length) also known as Skuterud catchment, located in the municipality of Ås (20,000 people), 30 km southeast of Oslo (Fig. 1). The catchment area is of approx. 4.5 km2 and largely consists of farmlands (60%) and forest/marshlands (31%). Grytelandsbekken has previously been studied for spatial variations in the microbial communities of various lotic freshwater ecosystems in different regions of Norway21. Based on that study, it was characterised as a rural creek with the highest diversity and abundance of microbial communities. Thus, in this follow-up study, Grytelandsbekken, specifically, was selected for assessing the seasonal dynamics of lotic freshwater bacterial communities.
All field measurements and samplings were carried out over a 2-year period at the same site of the creek, i.e. at about 0.25 km. In total, there were 16 sampling events, split equally between cold and warm seasonal regimes. The cold season refers to December–March, while the warm season includes June–September. In detail, there were four events during the cold season in the first year, Cold_1 (Cold_1a/Dec.2014, Cold_1b/Jan.2015, Cold_1c/Feb.2015, and Cold_1d/Mar.2015) and four events during the cold season in the second year, Cold_2 (Cold_2a/Dec.2016, Cold_2b/Jan.2017, Cold_2c/Feb.2017, and Cold_2d/Mar.2017). A similar sample assembly was applied to the warm seasonal regimes, i.e. four in the first year, Warm_1 (Warm_1a/Jun.2015, Warm_1b/Jul.2015, Warm_1c/Aug.2015, and Warm_1d/Sep.2015) and four in the second year, Warm_2 (Warm_2a/Jun.2016, Warm_2b/Jul.2016, Warm_2c/Aug.2016, and Warm_2d/Sep.2016). The entire study duration and sampling sets were conceived based on similar settings reported in aquatic research worldwide, profiling microbial communities through seasonal and spatial variations20,35,36.
Environmental measures
Primary physico-chemical parameters that are routinely measured in standard catchment water quality control37,38 were selected for abiotic characteristics of the lotic environment. These included organic matter content (expressed as CODCr), TSS, nutrients (Ptot and Ntot), EC, pH, and Temp. The latter was measured in situ at one of the weather stations, administrated by the Norwegian Institute of Bioeconomy Research (NIBIO) and located at the field measuring/sampling site of Grytelandsbekken. The station is equipped with a number of automatic sensors, providing hourly measurements and registering various climatic data online, which are all open and available at the AgroMetBase hosted by NIBIO (https://lmt.nibio.no/agrometbase/getweatherdata.php). The former parameters were analysed post-sampling in an accredited laboratory of the ALS Laboratory Group Norway AS. These analyses were performed in accordance with ISO and national standards for the respective parameters: CODCr (ISO 15705), TSS (CSN EN 872, NS 4733), Ptot (ISO 6878, ISO 15681-1), Ntot (EN 12260), EC (EN 27 888, SM 2520B, EN 16192), and pH (ISO 10523, EPA 150.1, EN 16192).
Genomic DNA purification: 16S rRNA amplicon library preparation and MiSeq sequencing
The seasonal water samples underwent DNA extraction using the QIAGEN DNeasy PowerWater Kit (QIAGEN GmbH, Hilden, Germany). In practice, 0.5 L of water was subjected to ultrafiltration to obtain a solid mass on a membrane filter (0.45 µm). DNA was then extracted from the collected filters, following the manufacturer’s instruction. The concentration and quality of the purified DNA was analysed on the NanoDrop Spectrophotometer (Thermo Fisher Scientific, Wilmington, USA). Five nanograms of purified DNA was used in PCR to prepare an amplicon library using the NEXTflex 16S V4 Amplicon-Seq Kit 2.0 (Bioo Scientific Corporation, Austin, TX, USA), following the provided protocol, which has previously been described in detail21. The applied specific 16S V4 forward primer (5′-GACGCTCTTCCGATCTTATGGTAATTGTGTGCCAGCMGCCGCGGTAA-3′) and reverse primer (5′-TGTGCTCTTCCGATCTAGTCAGTCAGCCGGACTACHVGGGTWTCTAAT-3′) were provided by the manufacturer and included in the sequencing kit. The library was prepared in triplicate for each sample. The final concentration of each library was measured on a Qubit Fluorometer (Life Technologies, Eugene, OR, USA) using the Quant-IT dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA). Pooling of all prepared libraries was achieved after concentration normalization using the SequalPrep Normalization Plate Kit (Thermo Fisher Scientific, Wilmington, USA). The pooled library was analysed on the Agilent 2100 Bioanalyzer system using the Agilent High Sensitivity DNA Kit (Agilent, CA, USA). It indicated a single band/unique product at 451 bp. The multiplexed library was sequenced on the Illumina MiSeq system using the MiSeq Reagent Kit V3, 600 Cycles (Illumina Inc., San Diego, CA, USA), following the default standard procedures.
Sequence data processing
The output sequence datasets were analysed using Microbial Genomics Module 2.0, added onto the CLC Genomic Workbench 10.1.1 (CLC Bio, QIAGEN Company, Aarhus, Denmark, https://www.qiagenbioinformatics.com/products/clc-genomics-workbench). The processing workflow consisted of four key components: quality filtration, OTU clustering, and alpha and beta diversity measures. Adapter and primer sequences were trimmed. Unqualified reads were trashed when the quality score was less than 20 or a higher number of ambiguous nucleotides (more than two) were detected. The average length after trimming was between 220–230 bp. Chimeric sequences and singletons were detected and discarded. The remaining qualified reads were used to characterize OTUs based on a reference database (Greengenes v_13_5)39 at a 97% identity level. The bacterial alpha diversity of each sample was estimated in rarefaction analysis with a depth cutoff at 50,000 reads. The bacterial beta diversity applied the Euclidean distance criterion (EDC) to estimate the community similarities between the examined samples. All sequence data are available at NCBI Sequence Read Archive, under accession number SRR10835654-669, as part of BioProject PRJNA599104.
Statistical analyses
Alpha diversity differences were tested using a two-tailed Student’s t-test at 0.05 significance. This was performed in the XLSTAT-ECOLOGY statistical software package version 2019.1.1 (Addinsoft 2020, Boston, USA, https://www.xlstat.com). A hierarchical clustering heat map was created to elucidate the bacterial community similarities/relatedness (pairwise) among all examined samples. It was conducted on a subset of top 700 OTUs, based on the EDC using the trimmed mean of M-values and Z-score (standard deviation numbers from the population mean) normalizations. This was further supported by the PERMANOVA analysis, performed to ascertain statistical significance of the clusters (p-value = 0.00001). These tests were executed using a package of functional features included in the Microbial Genomics Module (CLC Bio, QIAGEN Company, Aarhus, Denmark, https://www.qiagenbioinformatics.com/products/clc-genomics-workbench). Furthermore, the LEfSe tool40 was applied to identify the responsible bacterial members, accounting for community discrepancy between cold and warm seasons. The formatted abundance table of bacterial classes was uploaded to the Galaxy/Hutlab application web-based platform (Biostatistics Department, Harvard School of Public Health, Boston, MA, USA, https://huttenhower.sph.harvard.edu/galaxy) for pairwise comparison. Statistical significance of the comparison was determined by the Wilcoxon rank-sum test at the alpha value of 0.05. The identified features characterising the microbial differences among samples were processed using LDA, with a threshold score set at 2.0. Beyond that, RDA was carried out to determine the key abiotic environmental variables driving seasonal changes of the bacterial community based on the Pearson correlation test, with a statistical significance level higher than 95% (p  More

Phenological and climatic data
We used data from the Pan European Phenology Project (PEP725)29, an open-access database with long-term plant phenological observations across 25 European countries (http://www.pep725.eu/). The regional/national network partners of PEP725 are following a consistent guideline for phenological observations30 and prepare the data for submission to the PEP725 database curators29. We selected 30 species for which sufficient observational data were available: 21 deciduous broadleaved trees or shrubs, 6 herbaceous perennials, 2 evergreen coniferous trees, and 1 deciduous coniferous tree (Supplementary Table 3). Particularly, our data set included one fruit tree (Prunus avium) and one nut tree (Corylus avellana) since some of the chilling models are specifically developed for fruit and nut trees. A total of 2,493,644 individual records from 15,533 phenological stations were used. The stations were mainly distributed in moderate climates in Central Europe (Supplementary Fig. 3). Four spring events based on the BBCH code were investigated: BBCH 10, 11, 60, and 69, representing first leaves separated, first leaves unfolded, first flowers open, and end of flowering, respectively31.
We used the E-OBS v19.0eHOM data set32 with a spatial resolution of 0.1 × 0.1° for 1950–2018 for calculating CA and HR of the in situ phenological records. This data set is provided by the European Climate Assessment & Data set project and includes homogenized series of daily mean, minimum, and maximum temperatures. We also use the daily maximum and minimum temperature data from the GHCN data set33 to assess the scale effect. The GHCN data set contains station-based measurements from over 90,000 land-based stations worldwide, but only parts of PEP725 stations match with the GHCN stations.
For future climatic data (2019–2099), we used daily minimum and maximum temperatures simulated by the HADGEM2-ES model (with a spatial resolution of 0.5 × 0.5°) under two climatic scenarios (RCP 4.5 and RCP 8.5). These data have been bias-corrected by applying the method used in the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP)34, which were available on the ISIMIP server (https://esg.pik-potsdam.de/projects/isimip2b/).
Chilling models
We used 12 chilling models to measure the amount of chilling. One type of chilling model is based on several specific temperature thresholds. The most commonly used model, developed in the 1930s and 1940s for peach35, calculates the number of hours or days with temperatures 0 °C for calculating CA11,38. Models C1–C6 were developed based on various combinations of the upper and lower temperature limits. The rate of chilling was 1 for daily temperatures 5,{mathrm{or}},{T} < - 10} hfill end{array}} right.,$$ (2) $${mathop{rm{CU}}nolimits} _3 = left{ {begin{array}{*{20}{l}} 1 hfill & {0 le {mathop{T}nolimits} le 5} hfill \ 0 hfill & {{T} > 5,{mathrm{or}},{T} < 0} hfill end{array}} right.,$$ (3) $${mathrm{CU}}_{mathrm{4}} = left{ {begin{array}{*{20}{c}} 1 & {{mathop{T}nolimits} le 7} \ 0 & {{T} > 7} end{array}} right.,$$
(4)

$${mathop{rm{CU}}nolimits} _5 = left{ {begin{array}{*{20}{l}} 1 hfill & { – 10 le {mathop{T}nolimits} le 7} hfill \ 0 hfill & {{T} > 7,{mathrm{or}},{T} < - 10} hfill end{array}} right.,$$ (5) $${mathop{rm{CU}}nolimits} _6 = left{ {begin{array}{*{20}{l}} 1 hfill & {0 le {mathop{T}nolimits} le 7} hfill \ 0 hfill & {{T} > 7,{mathrm{or}},{T} < 0} hfill end{array}} right.,$$ (6) where CUi is the rate of chilling for Model Ci, and T is the daily mean temperature (°C). Model C7 is also known as the Utah Model39, which assigned different weights to different ranges of temperatures and was first used to measure the chilling requirements of peach (Eq. (7)). The Utah Model was modified to produce Model C8 (Eq. (8)), which removed the negative contributions of warm temperatures to accumulated chilling20. $${mathop{rm{CU}}nolimits} _7 = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le 1.4} hfill \ {0.5} hfill & {1.4 < {mathop{T}nolimits} le 2.4} hfill \ 1 hfill & {2.4 < {mathop{T}nolimits} le 9.1} hfill \ {0.5} hfill & {9.1 < {mathop{T}nolimits} le 12.4} hfill \ 0 hfill & {12.4 < {mathop{T}nolimits} le 15.9} hfill \ { - 0.5} hfill & {15.9 < {mathop{T}nolimits} le 18} hfill \ { - 1} hfill & {{mathop{T}nolimits} > 18} hfill end{array}} right.,$$
(7)

$${mathop{rm{CU}}nolimits} _8 = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le 1.4} hfill \ {0.5} hfill & {1.4 < {mathop{T}nolimits} le 2.4} hfill \ 1 hfill & {2.4 < {mathop{T}nolimits} le 9.1} hfill \ {0.5} hfill & {9.1 < {mathop{T}nolimits} le 12.4} hfill \ 0 hfill & {{mathop{T}nolimits} > 12.4} hfill end{array}} right.,$$
(8)

where CUi is the rate of chilling for Model Ci, and T is the daily mean temperature (°C).
Model C9 is a dynamic model developed for peach in Israel and South Africa40,41 and now adopted for apricot cultivars42. The most important characteristic of Model C9 was that a previous intermediate product affected the rate of chilling in the current hour or day. We did not provide equations for Model C9 for simplicity (see the equation in Luedeling et al.20).
Harrington et al.43 summarized published results for chilling units and constructed a chilling function based on a three-parameter Weibull distribution, coded as Model C10 (Eq. (9)). Model C11 has a triangular form, which was fitted by Hänninen44 using previous experimental results for Finnish birch seedlings (Eq. (10)). Zhang et al.45 recently fitted observational data to the triangular model for 24 plant species and found that a mean optimal chilling temperature of 0.2 °C and an upper limit of the chilling temperature of 6.9 °C were most effective. Model C12, therefore, uses the triangular form with parameters of 0.2 and 6.9 °C (Eq. (11)).

$${mathop{rm{CU}}nolimits} _{10} = left{ {begin{array}{*{20}{l}} 1 hfill & {2.5 < {mathop{T}nolimits} < 7.4} hfill \ 0 hfill & {{mathop{T}nolimits} < - 4.7,{mathrm{or}},{mathop{T}nolimits} > 16} hfill \ {3.13left( {frac{{{mathop{T}nolimits} + 4.66}}{{10.93}}} right)^{2.10}{mathop{rm{e}}nolimits} ^{ – left( {frac{{{mathop{T}nolimits} + 4.66}}{{10.93}}} right)^{3.10}}} hfill & {{mathrm{else}}} hfill end{array}} right.,$$
(9)

$${mathrm{CU}}_{11} = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le – 3.4,{mathop{rm{or}}nolimits} ,{mathop{T}nolimits} ge 10.4} hfill \ {frac{{{mathop{T}nolimits} + 3.4}}{{5 + 3.4}}} hfill & { – 3.4 < {mathop{rm{T}}nolimits} le 5} hfill \ {frac{{{mathop{T}nolimits} - 10.4}}{{5 - 10.4}}} hfill & {5 < {mathop{T}nolimits} < 10.4} hfill end{array}} right.,$$ (10) $${mathop{rm{CU}}nolimits} _{12} = left{ {begin{array}{*{20}{l}} 0 hfill & {{mathop{T}nolimits} le - 6.5,{mathop{rm{or}}nolimits} ,{mathop{T}nolimits} ge 6.9} hfill \ {frac{{{T + }6.5}}{{6.9 - 0.2}}} hfill & { - 6.5 < {mathop{T}nolimits} le 0.2} hfill \ {frac{{6.9 - {T}}}{{6.9 - 0.2}}} hfill & {0.2 < {T} < 6.9} hfill end{array}} right.,$$ (11) where CUi is the rate of chilling for Model Ci, and T is the daily mean temperature in °C. Forcing models Forcing models were used to measure HR for the spring events of plants. The GDD model is the most commonly used forcing model, which assumes that the rate of forcing is linearly correlated with temperature if the temperature is above a particular threshold. We mainly used Model F1 (Eq. (12)), which adopts a temperature threshold of 0 °C46,47,48, for examining the relationship between CA and HR: $${mathrm{FU}}_{mathrm{1}} = {mathrm{max}}({T},0),$$ (12) where FU1 is the rate of forcing for Model F1, and T is the daily mean temperature (°C). We also validated the chilling models by correlating them with HR based on seven other forcing models to assess the impact of the choice of forcing model on the results. Model F2 (Eq. (13)) is also a GDD model but has a temperature threshold of 5 °C12,18: $${mathrm{FU}}_{mathrm{2}} = {mathrm{max}}(T - 5,0),$$ (13) where FU2 is the rate of forcing for Model F2, and T is the daily mean temperature (°C). Piao et al.48 found that leaf onset in the Northern Hemisphere was triggered more by daytime than nighttime temperature. They thus proposed a GDD model using maximum instead of mean temperature. Models F3 (Eq. (14)) and F4 (Eq. (15)) are based on maximum temperature with thresholds of 0 and 5 °C, respectively. $${mathop{rm{FU}}nolimits} _3 = max ({mathop{T}nolimits} _{rm{max}},0),$$ (14) $${mathop{rm{FU}}nolimits} _4 = max ({mathop{T}nolimits} _{rm{max}} - 5,0),$$ (15) where FUi is the rate of forcing for Model Fi, and Tmax is the daily maximum temperature (°C). A recent experiment demonstrated that the impact of daytime temperature on leaf unfolding for temperate trees was approximately threefold higher than the impact of nighttime temperature49. Model F5 (Eq. 16) thus uses two parameters (0.75 and 0.25) to weigh the impact of daytime and nighttime temperatures on HR. $${mathop{rm{FU}}nolimits} _5 = 0.75 times max ({mathop{T}nolimits} _{rm{max}} - 5,0){mathrm{ + }}0.25 times max ({mathop{T}nolimits} _{rm{min}} - 5,0),$$ (16) where FU5 is the rate of forcing for Model F5. Tmax and Tmin are the daily maximum and minimum temperatures (°C), respectively. Many studies have suggested that the rate of forcing followed a logistic function of temperature44,50. Model F6 (Eq. (17)) uses a logistic function proposed by Hänninen44, and Model F7 (Eq. (18)) uses another logistic function proposed by Harrington et al.43. $${mathop{rm{FU}}nolimits} _6 = left{ {begin{array}{*{20}{l}} {frac{{28.4}}{{1 + {mathop{rm{e}}nolimits} ^{ - 0.185({mathop{T}nolimits} - 18.5)}}}} hfill & {{mathop{T}nolimits} > 0} hfill \ 0 hfill & {{mathop{rm{else}}nolimits} } hfill end{array}} right.,$$
(17)

$${mathop{rm{FU}}nolimits} _7 = frac{1}{{1 + {mathop{rm{e}}nolimits} ^{ – 0.47{mathop{T}nolimits} {mathrm{ + 6}}{mathrm{.49}}}}},$$
(18)

where FUi is the rate of forcing for Model Fi, and T is the daily mean temperature (°C).
Model F8 is a growing degree hour (GDH) model, where species have an optimum temperature for growth and where temperatures above or below that optimum have a smaller impact51. Model F8 (Eq. (19)) was first designed for calculating HR at hourly intervals, but we applied it at a daily interval. The stress factor in the original GDH model was ignored, because we assumed that the plants were not under other stresses.

$${mathrm{FU}}_{mathrm{8}}{mathrm{ = }}left{ {begin{array}{*{20}{l}} {mathrm{0}} hfill & {{T < }T_{mathrm{L}},{mathrm{or}},{T > }T_{mathrm{c}}} hfill \ {frac{{T_{mathrm{u}} – T_{mathrm{L}}}}{{mathrm{2}}}left( {{mathrm{1 + cos}}left( {pi {mathrm{ + }}pi frac{{{T} – T_{mathrm{L}}}}{{T_{mathrm{u}} – T_{mathrm{L}}}}} right)} right)} hfill & {T_{mathrm{L}} ge {T} ge T_{mathrm{u}}} hfill \ {left( {T_{mathrm{u}} – T_{mathrm{L}}} right)left( {{mathrm{1 + cos}}left( {frac{{uppi }}{{mathrm{2}}}{mathrm{ + }}frac{{uppi }}{{mathrm{2}}}frac{{{T} – T_{mathrm{u}}}}{{T_{mathrm{c}} – T_{mathrm{u}}}}} right)} right)} hfill & {T_{mathrm{u}}{ < T} le T_{mathrm{c}}} hfill end{array}} right.,$$ (19) where FU8 is the rate of forcing for Model F8, T is the daily mean temperature (°C), TL = 4, Tu = 25, and Tc = 36. Analysis We assessed the ability of each chilling model to represent long-term trends in the chilling conditions by calculating CA using each chilling model for each station for 1951–2018. CA was calculated as the sum of CUi from 1 November in the previous year to 30 April. The trend of CA at each station was visualized as the slope of the linear regression of CA against year. We also calculated Pearson’s r between each pair of chilling models for each station to determine if the chilling models were interrelated. We calculated HR and CA of spring events for each species, station, and year to determine if the chilling models are consistent with the physiological assumption that the reduction in chilling would increase HR (Fig. 1). HR was calculated as the sum of FUi from 1 January to the date of onset of spring events using Model F1, and the performances of the other forcing models (F2–F8) were also tested. We also compared 1 January with the other two starting dates of temperature accumulation (15 January and 1 February) to test any potential difference causing by the date when temperature accumulation begins. CA was calculated as the sum of CUi from 1 November in the previous year to the date of onset of spring events. We chose 1 November as the start date for CA because the endodormancy of temperate trees began around 1 November52. We only tested the linear relationship because the data were better fitted by the linear regression than the exponential model (Fig. 3), even though CA was linearly or nonlinearly negatively correlated with HR17. Pearson’s r between CA and HR for all records was calculated for each species, with a significantly negative Pearson’s r (p  More

Acoustic data and analysis
Audio data were collected from online sound archives (Xeno-Canto—https://www.xeno-canto.org—and Macaulay libraries—https://www.macaulaylibrary.org), creating a pool of over 2000 audio tracks. We assessed the sound quality of these audio tracks by listening and through visual inspection of sound spectrograms. To capture intra-specific variation, we limited audio extraction to one drum per audio track (which also avoided pseudoreplication) and only included species for which at least 3 high-quality drums could be extracted. We retained 736 high-quality drums suitable for further analyses. These drums were distributed among 92 species (out of the 217 recognized species of woodpeckers50 and 22 genera, providing a representative sampling of the phylogenetic diversity found in this family (Fig.1b)). Background noise and other artifacts were reduced by wavelet continuous reconstruction (R ‘WaveletComp’ package72), following the methods and description outlined in previous work73. The full script is available on demand. Finally, 22 acoustic variables were extracted from these filtered sound samples using the R ‘Seewave’ package74. Given the pulsed-like nature of drumming, the chosen variables emphasized the temporal and amplitude-related (all normalized to the maximal amplitude within a given drumming signal) features of the sounds (Supplementary Table 1). These 22 variables were z-scored and then used in all subsequent analyses. Since these variables were partly correlated to varying degrees, we performed a principal component analysis (PCA) to reduce the number of descriptive variables quantifying drumming acoustic structure. This dimensionality reduction was useful for visualization and necessary for regularization (i.e. to prevent overfitting) in many of our analyses. This resulted in six principal components (PCs) with eigenvalues >1 which together explained 75% of the variance (Supplementary Table 11).
We used these variables to evaluate the similarity between species-specific drums, by performing a hierarchical cluster analysis (HCA)75,76 based on Euclidean distances and the ‘Ward.D2’ method (‘NbClust’ R package77). NbClust provides a clustering output resulting from the use of multiple indices (in the case of our analysis, 26 indices were used). The best number of clusters is chosen according to the majority rule, i.e. it is the one supported by the highest number of indices used. This entails creating a vector of acoustic features (22 raw acoustic measures or 22 PCs) for each of the 92 species in our dataset and calculating the Euclidean distance between these vectors to evaluate how close acoustically species were. Note that one can still use Euclidian distances in a non-orthonormal space to calculate the ‘distance’ between signals. The result is a distance metric that might give more weights to measures that co-vary. This could theoretically affect the clustering results. However, when we performed the same analysis using the 22 PCs, we obtained the same grouping (6 clusters) with very minor differences in species grouping and distances between clusters as shown in the relative length of branches (Supplementary Fig. 13). The output of this HCA established an optimal classification of woodpeckers’ drums into 6 main drumming types (Fig. 2a), described as follows:
– Acceleration (AC): Beak strikes decrease in amplitude as they are produced within successively shorter time intervals.
– Regular sequence (RS): Beak strikes are produced in bouts, each comprising a relatively fixed (stereotyped) number of strikes.
– Irregular sequence (IS): Beak strikes are produced in bouts, each comprising a variable number of strikes (as opposed to RS).
– Steady fast (SF): Beak strikes are produced with constant time intervals and at a similar amplitude, with a high pulse rate (on average >20 strikes/s).
– Steady slow (SS): Beak strikes are produced with constant time intervals and at a similar amplitude, with a low pulse rate (on average 2; see models Supplementary Tables 3–6), a likelihood ratio test (LRT) was conducted to test for the specific effect of predictor variables (Supplementary Table 4). Because both models (null and fitted) differ in their fixed effects, model comparison was performed on models fit by maximum likelihood (ML) with the phylogenetic correlation structure (Pagel’s λ) fixed to the estimates obtained from initial fit by REML. The statistics reported for model comparison are likelihood ratios.
PGLS models testing for a relationship between life-history variables and drumming structure included either of PC1 to PC6 as the dependent variable to investigate whether differences exist between these proxies for acoustic structure. Similarly, LDs were used to verify our results with these different loading combinations of drumming acoustic variables. No significant correlations were found between life-history traits and acoustic structure using LDs instead of PCs (Supplementary Table 6), indicating that the combination of structural variation captured by the LDs differed from that of the PCs, while not leading to fundamentally different conclusions. Similarly, no significant correlations were found between life-history traits and information content (no decrease in AICc >2; Supplementary Table 3), overall emphasizing that none of the variables investigated here (and which could have potentially affected drumming structure) seemed to have influenced species-specific information, or at least not directly.
Ancestral states reconstructions
We carried out two types of ancestral state reconstructions: discrete reconstruction of drumming status in Fig. 1b and of drumming types in Fig. 3a (using ‘ace’ from the R ‘ape’ package81), or continuous character reconstruction of drumming acoustic structure based on Brownian motion models (using ‘fastanc’ from the R ‘phytools’ package82 in Supplementary Figs. 4 and 5) and using relaxed Brownian motion model (using ‘rjmcmc.bm’ from the R ‘geiger’83 package in Fig. 4a and Supplementary Figs. 8 and 9).
While evaluating the likelihood that drumming was already present at an early stage of woodpecker’s phylogeny, we tried to represent the most complete tree of the family, based on very recent molecular data50. Note that strictly speaking, we evaluate the state at the root but at the next internal node, i.e. at the node including Picumninae and Picinae (the largest pie-chart in our Fig. 1b), as Wrynecks do no drum, and neither do honeyguides or barbets). To include species with an unknown drumming status in this discrete reconstruction, we attributed equal probability distribution between the 3 states (i.e. when the ‘drummer state’ of a species is unknown, the species is given, prior to ancestral state reconstruction, a 1/3 probability of belonging to each of the three categories ‘drummer’, ‘occasional drummer’ and ‘non-drummer’). Stochastic mapping was performed under an MCMC model, sampling the rate matrix from its posterior distribution for Q (‘Q = mcmc’ in make.simmap function from the R ‘phytools’ package), with an equiprobable default prior at the root, and 200 simulations. Under a symmetrical model for the probability to change among the three states, scaled likelihood on woodpeckers’ ancestral node indicated 56.4%, 38.3% and 5.3% probabilities of being a drummer, an occasional drummer and a non-drummer, respectively. This is in line with the fact that morphological adaptations for drilling (including reinforced rhamphotheca, frontal overhang and processus dorsalis pterygoidei) evolved in the ancestral lineage of Picumninae and Picinae64.
To prevent overfitting, the discrete reconstructions for drumming types were estimated for six different rate models: equal rate model (ER), symmetric rate model (SYM), all rates difference model (ARD) and three sequential transition models based on the normalized MIL as measures of complexity as shown in Supplementary Fig. 7 ((SF leftrightarrow SS leftrightarrow DK leftrightarrow AC leftrightarrow RS leftrightarrow IS)). These three models assumed (1) sequential and equal, (2) sequential and incremental and (3) sequential and reversed transition rates, respectively. The number of parameters for these 6 rate models were 25, 1, 15, 1, 2 and 10. The final regularized likelihoods of each ancestral states were then obtained by model averaging using Akaike weights.
Calculation of information at different evolutionary steps was carried out as an extension of the drumming types reconstruction described above. From the discrete ancestral reconstruction procedure, probability distributions of drumming types were obtained for each node of the phylogenetic tree. We then obtained probability distributions at 20 fixed time intervals (dt = 1 myr) by linear interpolation. Using these probability distributions, we sampled drumming types proportionally from extant species descending the node closest to the time interval to estimate ancestral information values. This bootstrap procedure was repeated 30 times in order to obtain reliable estimates of mean and standard error. In this manner, we obtained information-through-time plots. These plots quantify a putative diversity of drumming signals in the clade at a particular point in time. They are similar in spirit to the disparity-through-time plots that have been used to measure specific morphological diversity in a clade through time using phylogenetic trees based on molecular data in combination with morphological measures in extant species84.
Continuous ancestral character trait reconstruction of drumming acoustic structure was carried out using either the six PCs that explain variation among the 22 drumming acoustic variables, or the six LDs that explain the variation in discriminating potential among the same variables (see above, ‘Acoustic data and analysis’ and ‘Calculation of information’ sections; Supplementary Figs. 4 and 5). The results and conclusions were similar for all PC’s and since the PC1 component has strong loading of multiple acoustic variables and the highest acoustic structure variance explained (Supplementary Table 11) it serves well as an illustrative example. The measure of phylogenetic signal on continuous traits (i.e. the historical contingency between species-specific drums that renders a trait non-randomly distributed along the phylogenetic tree) was made using Pagel’s lamba (Supplementary Table 2).
Reconstructing information content from raw MIL values would not have been biologically relevant since information calculation is based on the number of species involved, a factor that changes as branches merge going backward along the phylogenetic tree. We thus reconstructed MIL based on the normalized MIL values to avoid this pitfall. We used a Bayesian model implemented in the R package ‘Geiger’83 (model ‘rbm’ in the function ‘rjmcmc.bm’) to estimate branch-specific rates of trait evolution (i.e. changes in rates through time and across lineages). In this method, a reversible jump Markov Chain Monte Carlo (MCMC) sampling algorithm is used to detect shifts in rates of continuous traits evolution under a relaxed Brownian motion model85. The results of the model fit were summarized by the branch-specific average rate, estimated from the posterior samples. To obtain relative variations in posterior average rates, drumming structure (PC1-PC6) and MIL were standardized, i.e. these traits were divided by their standard deviation prior to running the ‘rbm’ models.
Analytical simulations of selection for information
In Fig. 3c, we compared that reconstructed evolution of information to what might be expected in different scenarios to further support those conclusions. More specifically, we estimated the ancestral MI for two simulated scenarios using an analytical model that describes species-specific information based on the probability of correct detection and the number of species (see ‘Calculation of information’ section). In the ‘No Diversifying Selection’ scenario (dark brown), the probability of correct detection for the initial pair of species, p2, is first estimated from the data using the approach described in the main text. It is then assumed that that additional species are randomly just as different/similar than these original species pair, yielding a probability of correct detection through time given by (p_c(t) = p_2^{n_s(t) – 1}), where ns(t) is the number of species at a given time. In the ‘Strong Diversifying Selection’ scenario (light brown), the probability of correct detection estimated at the first time point in our reconstruction (−20 M years ago, 3 species) is kept constant, (p_c(t) = p_2). In other words, the only species that survive would be species that can discriminate themselves from all other species equally well than the currently existing species. The reconstructed (actual) scenario is found between these two extreme values, showing that the drumming types are clearly not random but were also not under high evolutionary pressure to increase species-specific information. New drumming types evolved and species within types used signals that were distinct enough to result in the maintenance of normalized MI.
In Supplementary Fig. 6, we showed that the non-normalized reconstructed MI increased more rapidly when new drumming types appeared but that the normalized MI was relatively constant, reflecting the fact that the appearance of novel drumming types could co-occur with rapid radiation and increase in species numbers.
Playback experiments
Initial preparation involved identifying and mapping the areas prone to high densities of great-spotted woodpeckers Dendrocopos major, the study species of this experimental phase, using GIS maps provided by the LPO (French Bird Protection Organization). D. major is commonly found in European forests, ranging from open coniferous to mature deciduous forests. Playback experiments were carried out on wild individuals around Saint-Etienne, France, during this species’ breeding season (February–April 2017). All experiments were performed in accordance with relevant guidelines and regulations including French national guidelines, permits and regulations regarding animal care and experimental use (approval no. D42-218-0901, ENES lab agreement, Direction Départementale de la Protection des Populations, Préfecture du Rhône).
Two sets of experiments were conducted over the course of the breeding season, although we implemented the same general design which consisted in simulating a territorial intrusion. Playback stimuli tracks consisted of eight drums spread unevenly over about 60 s, aiming at representing the variation encountered in natural sequences (ref. 44 and personal observations). The first experiment (Exp. 1) aimed at investigating D. major’s response to conspecific vs. heterospecific drums. The other experiment (Exp. 2) aimed at investigating D. major’s response to drums from conspecifics vs. drums modified through acoustic manipulation (i.e. signal re-synthesis). D. major typically drums with an ‘acceleration’ pattern, which is mainly characterized by a shortening of the inter-strike time interval, a progressive decrease in strikes’ amplitude, and a gradual change in spectral properties as strikes get faster and weaker.
In Exp. 1, we used a paired and randomized order design, presenting each focal individual with one D. major drum and one drum from one out of 4 different species: 2 of which have very different drumming patterns (Picus canus and Dryobates minor, both producing ‘steady fast’ drums), and 2 others which have similar (accelerating) drumming patterns (Dendrocopos syriacus and Dendrocopos hyperythrus). A potentially confounding factor (which is nevertheless in line with our phylogenetic analyses) lies in that the allopatric species producing a drum similar to that of D. major also happened to be closely related to our model species. We carried out 48 playback experiments (testing 24 individuals with one of 4 categories of paired signals).
In Exp. 2, we altered one of the 3 acoustic features described above or all of them together (thus having 4 categories of modified signals), using Praat sound analysis software86. The design was paired so that each focal individual was exposed to one conspecific drum and one modified drum, following a randomized presentation order. This led to 48 playback experiments (24 individuals, each tested with one of 4 categories of paired signals).
Within each of Exp. 1 and Exp. 2, tested individuals were all separated by at least 500 m, ensuring different identities since their territory sizes vary between 200 and 400 m87,88. Upon visual or aural detection of (an) active individual(s), the experimenter set up an Anchor Megavox loudspeaker at about 1–1.5 m from ground level. The speaker was connected to an Edirol R-09 recorder (stimuli tracks were created and stored as WAV files, 44.1 kHz sampling frequency). Playbacks started at about 50 m from where the experimenter last saw or heard the focal individual. Following the work from Schuppe and colleagues89, playback intensity was calibrated and kept at about 80 dB measured 1 m away from the speaker. Behavioural data collection started when the first drum of the stimuli track was broadcasted and lasted 10 min from that moment on. To document focal individuals’ responses, notes were taken manually and continuously, while audio was recorded with a Sennheiser ME67 microphone mounted on a tripod and connected to a digital recorder (Zoom H4N, 44.1 kHz, 16 bit). If a response was elicited from multiple individuals in the area, only the one from a particular individual (ideally the one seen or heard before setting up the experiment) was monitored and used in further analyses. Six behavioural variables were reported, namely the number of screams, the number of drums, the approach (which was divided into three categories: ‘within 25 m’, ‘25–50 m’ and ‘further than 50 m’) as well as the latencies to first scream and drum and the latency to closest approach. When no occurrence was observed for the first three behaviours, latencies were set by default to the maximum value, i.e. the duration of the full experiment (10 min = 600 s). To characterize D. major’s behaviour, a PCA was then performed on scaled/centred data, where we retained the first principal component (‘Playback-PC1’) as an indicator of the behavioural response’s strength. A higher Playback-PC1 score indicates a stronger territorial response, i.e. more screams, a closer approach to the speaker and shorter latencies to these 2 behaviours. A second significant component resulted from this PCA, which represented the drumming’s response (inversely related: a higher Playback-PC2 score indicates fewer drums and a longer latency to drum; see Supplementary Table 13). None of the pairwise comparisons were statistically significant for PC2, besides a stronger drumming response to drums resynthesized without temporal variation than to D. major drums (Supplementary Fig. 15). This can be explained by the fact that birds were tested during their breeding season. At this time, drumming behaviour is likely to occur more consistently and commonly across experiments, independently from the stimulus played back, while screams and approach do not occur unless threat of intrusion is clear. Therefore, we used Playback-PC1 to represent birds’ behavioural response in our analysis (as it is in addition explaining much more variance in the behavioural data than Playback-PC2). Note that, as two playback sets were involved in this study, while we considered them independently in our statistical analysis, for standardization of the behavioural scale, we used the same polynomial equation. More specifically, the linear equation obtained from the loading scores of Exp. 1 was applied to the behavioural data of Exp. 2 for computation of Playback-PC2 scores.
Finally, distances were approximated during continuous note-taking and confirmed post-experimentally using a National Geographic 4*21 rangefinder (measurement accuracy: ±1 m up to 200 m). Sex was not documented as sometimes birds were not seen (but just heard drumming or calling back at our playback), which we nevertheless believe to be negligible since both sexes drum and are territorial in this monogamous species44,90.
Statistical analyses tested for differential responses of focal birds to drums of their own species versus either another species or a modified resynthesized condition. A paired comparison design was used by means of LMMs and contrasts using R software (‘lme4’ and ‘lsmeans’ packages)91,92. LMMs included study day and time, order of presentation and focal bird identity as random factors, and tested for a fixed effect of the interaction between treatment and group of paired condition. Contrasts were then computed between treatments (i.e. conspecific versus non-conspecific drums) for each group (i.e. each paired testing condition, such as D. major versus D. minor for which n = 6 birds were exposed to paired playback presentations—see Fig. 5a, b). Before contrasts and using the ‘lsmeans’ function, a Tukey adjustment for multiple testing was used; two-sided statistics are reported.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

Publication collection and data compilation
The detailed methods of publication search and data collection were described in our related work7. In brief, 10 databases in Web of Science (WoS; 1 January 1900 to 13 December 2016) including BIOSIS Previews, Chinese Science Citation Database, Data Citation Index, Derwent Innovations Index, Inspec, KCI-Korean Journal Database, MEDLINE, Russian Science Citation Index, SciELO Citation Index, and WoS Core Collection were used for searching peer-reviewed publications that reported GCMEs. The 18 keywords for WoS title search were: global change, climate change, free-air carbon dioxide enrichment, free-air CO2 enrichment, elevated carbon dioxide, elevated CO2, elevated atmospheric CO2, CO2 enrichment, eCO2, [CO2], warming, elevated temperature, changing precipitation, increased precipitation, decreased precipitation, nitrogen deposition, nitrogen addition, and nitrogen application. Through these search, 310,177 publication records that might be relevant to our topic were found.
First, we identified all the 310,177 records via reading each title. Second, we read the abstracts of all the records collected in the first step to further screen publications. During the two steps, we excluded 291,436 records because these studies were reviews/meta-analyses or conducted in non-terrestrial ecosystems such as oceans. Third, we read the methods of the remaining 18,741 publications to identify which of them met the following three inclusion criteria:
1.
Publications reported results of outdoor GCMEs which had at least three control and global change treatment plots ( > = 1 m2).

2.
The GCMEs were conducted in terrestrial ecosystems except for croplands and lab incubation studies.

3.
The GCMEs aimed to examine effects of simulated global change drivers on carbon, nitrogen, and water-cycle variables as well as plant and microbial parameters.

During the screening in the third step, 1,290 publications met these defined criteria.
We subsequently cross-checked the list of the 1,290 publications with references cited by the previous review/meta-analysis articles in global change research as well as the 1,290 publications, and collected 756 publications. In addition, 184 studies were collected by searching the websites of ecology laboratories and experiment networks and checking the references of the papers downloaded from these websites. In total, 2,230 publications were collected in the original version of the database7. Moreover, another 12 publications were found when we checked and reorganized all the data extracted from the 2,230 publications8,9,10,11,12,13,14,15,16,17,18,19. This database compiled 11 plant production and ecosystem carbon exchange variables including net primary productivity (NPP), above- and below-ground NPP (ANPP and BNPP), total biomass, aboveground biomass (AGB), root biomass, litter mass, gross and net ecosystem productivity (GEP and NEP), and ecosystem and soil respiration (ER and SR). Data of mean values, standard deviations or standard errors, and sample sizes (number of plot replications) of these variables in the control and treatment (e.g., elevated CO2, nitrogen addition, warming, increased/decreased precipitation, or their combinations) groups were extracted from each publication when possible. The figures were digitized using SigmaScan Pro 5.0 (SPSS, Inc.) and the numerical values were extracted when a publication presented experimental data graphically. Data of the experiments that were conducted over less than one year/growing season were excluded in this database. However, we included short-term data from tundra studies because most of measurements in this ecosystem were performed during July-August. Overall, 5,213 pairs (the control versus global change treatment) of plant production and ecosystem carbon exchange samples were collected in this database, having 2,247, 2,120, 81, and 765 pairs from single-, two-, three-, and four-factor manipulative experiments, respectively (Fig. 1).
Fig. 1

Number of samplings. Number of sample pairs of ecosystem carbon-cycling variables including net primary productivity (NPP), above- and below-ground NPP (ANPP and BNPP), total biomass, aboveground biomass (AGB), root biomass, litter mass, gross and net ecosystem productivity (GEP and NEP), and ecosystem and soil respiration (ER and SR) extracted from publications reporting single-, two-, three-, and four-factor global change manipulative experiments.

Full size image

Environmental metadata: Climate and vegetation
Information on the locations and altitudes of each experimental site, site climate including mean annual temperature (MAT) and precipitation (MAP) as well as wetness index ((frac{{rm{MAP}}}{{rm{MAT+10}}})), ref. 20, and vegetation types were extracted from each of the 2,242 publications. If a study did not report climate characteristics for its experimental site, data of MAT and MAP were downloaded from Climate Model Intercomparison Project phase 5 (CMIP5; https://esgf-node.llnl.gov/projects/cmip5/) based on the site coordinate. The dataset selection in CMIP5 was “historical (simulation of recent past 1850–2005)” and the climate data averaged from 20 (i.e. BCC_CSM1_1, BCC_CSM1_1_M, CANESM2, CCSM4, CMCC_CM, CMCC-CMS, CNRM-CM5, CSIRO_MK3_6_0, GFDL_CM3, GISS_E2_H, HADGEM2_AO, HADGEM2_ES, INMCM4, MIROC_ESM, MIROC_ESM_CHEM, MIROC5, MPI_ESM_LR, MPI_ESM_MR, MRI_ESM1, and NORESM1_M)21, that contained historical climate data, out of the 35 global climate models available in CMIP5 were used in this study. In addition, we downloaded data of climate means at global 1 × 1° land grid cells from Princeton University (http://hydrology.princeton.edu/data/pgf/v3/) to construct global climate space. Moreover, we classified ecosystems subjected to ecosystem manipulative experiments into five typical types: forests (mature forests and tree seedlings), grasslands (grasslands, meadows, short- and tall-grass prairies, temperate/semi-arid steppes, shrublands, savannas, pastures, and old-fields), tundra, wetlands (peatlands, bogs, marshes, and fens), and deserts.
Metadata of experimental facilities and performance
Information on CO2 enrichment and warming facilities were also extracted from the related publications reporting CO2 and warming effects on plant production and ecosystem carbon exchange. Facilities used in elevated CO2 experiments included greenhouse, open-top chamber, free-air CO2 enrichment, and tunnels. Warming experiments primarily used greenhouse, open-top chamber, soil heating cables, infrared radiator, and infrared reflector to elevate vegetation canopy and soil temperature. In addition, the manipulation magnitudes of global change drivers imposed by manipulative experiments, such as the increases in CO2 concentrations (ppm) and temperature (°C), the changes in precipitation amount (mm), and the rates of nitrogen input (g N m−2 yr−1), were also collected and added into this updated database. More

For sustainable use of marine-animal resources, preservation of endangered species, and conservation of ecosystems, it is very important to understand the biology of individual marine animal. From the viewpoints of physiology, ethology, and environmentology, marine animals have been studied by bioelectric measurement1,2,3,4, bio-logging5,6,7,8,9, and DNA (genome) analysis10,11,12,13,14,15, respectively. Recent technological innovations helped studies on bio-logging and DNA analysis advance rapidly, but advancement of bioelectric-measurement technology, which has existed for a long time, lags behind those of bio-logging and DNA analysis.
Now, aiming to obtain good harvests, the aquaculture industry requires bioelectric measurements to grasp the health condition of marine animals from pathophysiological viewpoints. Moreover, the electrocardiogram (ECG), which is a kind of bioelectric measurement, carries high expectations because it can evaluate psychological stress of marine animals just as it can evaluate that of humans16,17,18,19. Moreover, ECG can be used in fish ethological- and physiological studies2,4, so innovating techniques and devices for ECG measurement will contribute to developing these studies.
In regards to bioelectric measurement targeting marine animals, to prevent electric short-circuiting between the pair of bioelectrodes via seawater (which is conductive), one or multiple pairs of bioelectrodes are embedded inside the living body by incision surgery20,21, which can impose a heavy workload on inexperienced experimenters. Moreover, the animal can often become agitated without anesthesia and consume much physical energy when the electrodes are implanted into its body. To reduce these burdens, we propose a novel method of measuring bioelectric signals—which utilizes the conductivity of seawater surrounding the animal—by using only one bioelectrode attached at each measurement point (in contrast to the conventional method, which requires a pair of bioelectrodes). To the best of our knowledge, a similar method has not been reported.
In this paper, the proposed method of bioelectric measurement for marine animals under the seawater is first overviewed. Next, the bioelectric measurement system for the chosen experimental subjects, namely, fish, is described, and the availability of the proposed method is verified. Then, the experimental procedures and results of bioelectric measurements are presented. Finally, possible applications of the proposed method are discussed. More