Ecology

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

The decision by an individual to move from one area to another is mediated by a number of factors, such as resource quality and availability, predation risk and local environmental conditions, all of which will influence its survival and reproductive output1,4. The challenge for conservationists is understanding how these individual decisions can affect population dynamics, home ranges and ultimately species’ survival1.
Home ranges of carnivores should overlap and in some cases envelop those of their prey species. Southern Ground-hornbills feed on a variety of prey, ranging from snakes, rabbits and birds to invertebrates12,18. Through tracking Southern Ground-hornbill movements, we were able to show that group home ranges during the early and late dry seasons were larger than in the wet season. As the Southern Ground-hornbill breeding season in South Africa coincides with the warm, wet summer months, prey availability, especially that of invertebrates, is expected to be higher20,21, suggesting that individuals would not need to travel as extensively to find sufficient food. Furthermore, in the late dry season, groups used between 76 and 115% of their home ranges. This was likely a result of having to increase their search for food and relaxation of the central place foraging required around the nest during the breeding season.
Previous research on Southern Ground-hornbill home ranges has recorded group densities ranging from one group per 4000 ha (communal areas in Zimbabwe22), to one group every 10,000 ha (KNP14), with one group in the Limpopo Valley having a home range close to 20,000 ha21. These results were obtained by direct observations of active nest sites or using VHF radio transmitters. In our study using GPS data, we showed that home range sizes of Southern Ground-hornbills within KNP vary considerably. Despite this, our results confirmed the findings of Theron et al.21 and Zoghby et al.20, demonstrating a restricted and contracted home range during the breeding season, when group movements are concentrated around the nest site (central place foraging). Presumably, breeding success would influence the extent of wet seasonal home range for Southern Ground-hornbills, with groups abandoning their central place foraging behaviour when nests fail. Wyness19 reported that of four Southern Ground-hornbill groups studied in the Association of Private Nature Reserves (APNR) adjacent to the KNP, the three that bred successfully in the year of their study showed a breeding season range reduction to between 24–36% of their non-breeding home range. The unsuccessful group used 70% of their home range during this time19. Surprisingly, the groups within the KNP did not show such a definitive pattern in home range size reduction associated with breeding success, although all groups that attempted breeding did show a wet seasonal home range reduction. Of the six Southern Ground-hornbill groups monitored in our study, four groups bred successfully, one group’s attempt failed (Ngotso Camp), and the breeding status for the third group (Shingwedzi) was unknown. The groups that bred successfully used 21–97% of their respective home ranges, with the unsuccessful group using 85% of their home range (See Table 1).
Southern Ground-hornbills are known to favour more open habitats for foraging20,23. Our results supported this, with groups selecting the open woodland and grassland habitat types year-round, following their availability within the landscape.
Although Southern Ground-hornbill seasonal territory size differed significantly amongst the groups, they all showed a decrease in the amount of low shrubland and an increase in the amount of grassland habitat used with increased territory size. Similarly, as seasonal territory sizes increased, the amount of low-medium woody cover (25–50%) decreased. Thus, when selecting an area for a reintroduction of Southern Ground-hornbill groups, the ratio of low-medium woody cover (low shrubland) to grassland, calculated based on the national land cover datasets available, should be taken into account, as this will likely influence the home range size and the number of groups that could be supported in an area.
Although an understanding of the changes and restrictions in territory size is important for the management of a species, the types of movements adopted within a population will influence the management actions needed for their conservation, such as ensuring connectivity or access to certain resources1. Conservation policy and management actions are less effective when interventions do not integrate both the spatial and temporal changes in habitat use and the scale of species movements1,3. The results from the first-passage time analysis of Southern Ground-hornbill movements showed that the different groups did not consistently demonstrate seasonal patterns in the scale at which they concentrated their foraging efforts. The mean distances travelled for all trajectory paths, classified as active foraging behaviour, were similar and lower in the late wet and early dry seasons compared with the late dry and early wet seasons. Movement between foraging resource patches or mean relocation distances were highest in the wet season months, with the maximum mean distances travelled during the early wet season and the start of the breeding period. Overall prey abundance for Southern Ground-hornbills is generally higher in the wetter months, resulting in a decrease in relocation distances. Our results support the theory that Southern Ground-hornbill wet season movements are most likely influenced by the need to travel to and from the nest site to provision prey to the incubating female and growing nestling. Once resources closer to the nest are depleted, the distances travelled to access additional habitats and prey would likely increase.
Southern Ground-hornbills seemingly prefer nest sites surrounded by more open woodland habitat24,25. Habitat structure and the diversity of habitat types within a 3 km radius around the nest site positively influenced Southern Ground-hornbill nesting success. An increase in the density of woody habitat surrounding the nest site, however, had a negative impact on Southern Ground-hornbill breeding success24, possibly owing to decreased foraging opportunities, an increased risk of predation or an increase in foraging effort beyond a value which is beneficial.
Habitat structure will likely promote or inhibit the types of movement that can occur in an area. The results from the multinomial regression (Table 5) indicate that the likelihood of a movement behaviour being classified as “foraging” within the open woodland, grassland and dense thicket habitat types was higher than the behaviour being attributed to “relocating”. This is to be expected for open woodland, and grassland habitats as these are both ideal open foraging habitats for Southern Ground-hornbills20,23 and are used year-round in proportion to their availability. Southern Ground-hornbills spend around 70% of their day walking12 and have been shown to travel distances of up to 10.6 km in a day20. Having to navigate through dense thicket vegetation in an area may increase the amount of time spent there, possibly accounting for why this habitat type is predicted to be used more for “foraging”-type behaviour as opposed to “relocating” behaviour. Travel through areas of low shrubland habitat was considered “relocating” behaviour, suggesting that within this habitat type, it is more profitable for Southern Ground-hornbills to move further, and the corresponding chance of finding food greater, than conducting area-restricted searches and spending longer periods concentrated in one patch.
When comparing movements between habitats allocated to “resting” as opposed to “foraging”, the time spent in all habitats was most likely as a result of “foraging”. As GPS locations were only recorded during the day, switching off at dusk (~ 18h00) when Southern Ground-hornbills would roost for the night, habitat preferences for “resting” movements may not have been recorded. Moreover, during the day, Southern Ground-hornbills may not be actively selecting for specific habitat types in which to roost or rest. They may simply be roosting or resting at a chosen site to escape the midday heat within the habitat type in which they were “foraging” or “relocating”.
We were unable to explore differences in movement relating to specific characteristics of the tagged bird (age, sex, helper versus breeder status, etc.) in our study. However, future research should consider study designs able to account for these potential differences, as García-Jiménez et al.26 showed that both the breeding season and sex of the individual influence displacement and distance travelled in Pyrenean Bearded Vultures (Gypaetus barbatus). They found that all individuals travelled more in the breeding season, with females having greater cumulative and maximum distances regardless of the season. More

1.
Lomolino, M. V., Riddle, B. R. & Brown, J. H. Biogeography (Sinauer Associates Inc, Sunderland, 2006).
Google Scholar 
2.
Chamorro, D., Olivero, J., Real, R. & Muñoz, A.-R. Environmental factors determining the establishment of the African Long-legged Buzzard Buteo rufinus cirtensis in Western Europe. Ibis (Lond. 1859). 159, 331–342 (2017).
Google Scholar 

3.
Pearson, R. G. & Dawson, T. P. Predicting the impacts of climate change on the distribution of species: Are bioclimate envelope models useful ?. Glob. Ecol. Biogeogr. 12, 361–371 (2003).
Google Scholar 

4.
Sun, J. et al. Modeling the potential distribution of Zelkova schneideriana under different human activity intensities and climate change patterns in China. Glob. Ecol. Conserv. 21, e00840 (2020).
Google Scholar 

5.
Mackey, B. & Lindemayer, D. Towards a hierarchical the spatial distribution of animals. J. Biogeogr. 28, 1147–1166 (2011).
Google Scholar 

6.
Márquez, A. L., Real, R., Olivero, J. & Estrada, A. Combining climate with other influential factors for modelling the impact of climate change on species distribution. Clim. Change 108, 135–157 (2011).
ADS  Google Scholar 

7.
IPCC. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (2014) https://doi.org/10.1017/CBO9781107415324.004.

8.
Root, T. L. et al. Fingerprints of global warming on wild animals and plants. Nature 421, 57–60 (2003).
ADS  CAS  PubMed  Google Scholar 

9.
Karl, T. R. & Trenberth, K. E. Modern global climate change. Science (80-). 302, 1719–1723 (2003).
ADS  CAS  Google Scholar 

10.
Walther, G.-R. et al. Ecological responses to recent climate change. Nature 416, 389–395 (2002).
ADS  CAS  PubMed  Google Scholar 

11.
Lazo-Cancino, D. et al. The impacts of climate change on the habitat distribution of the vulnerable Patagonian-Fueguian species Ctenomys magellanicus (Rodentia, Ctenomyidae). J. Arid Environ. 173, 104016 (2019).
Google Scholar 

12.
Borzée, A. et al. Climate change-based models predict range shifts in the distribution of the only Asian plethodontid salamander: Karsenia koreana. Sci. Rep. 9, 11838 (2019).
ADS  PubMed  PubMed Central  Google Scholar 

13.
Thomas, C. D. & Lennon, J. J. Birds extend their ranges northwards. Nature 399, 213 (1999).
ADS  CAS  Google Scholar 

14.
Wallingford, P. D. et al. Adjusting the lens of invasion biology to focus on the impacts of climate-driven range shifts. Nat. Clim. Change. 10, 398–405 (2020).
ADS  Google Scholar 

15.
Castro, A., Muñoz, A.-R. & Real, R. Modelling the spatial distribution of the Tengmalm’s owl Aegolius funereus in its Southwestern Palaeartic limit (NE Spain). Ardeola 55, 71–85 (2008).
Google Scholar 

16.
Thomas, C. D. et al. Ecological and evolutionary processes at expanding range margins. Nature 411, 577–581 (2001).
ADS  CAS  PubMed  Google Scholar 

17.
Huntley, B. et al. Potential impacts of climatic change upon geographical distribution of birds. Ibis (Lond. 1859). 148, 8–28 (2006).
Google Scholar 

18.
Maclean, I. M. D. et al. Climate change causes rapid changes in the distribution and site abundance of birds in winter. Glob. Change. Biol. 14, 2489–2500 (2008).
Google Scholar 

19.
Elorriaga, J. & Muñoz, A.-R. First breeding record of North African Long-legged Buzzard Buteo rufinus cirtensis in continental Europe. Br. Birds 103, 399–401 (2010).
Google Scholar 

20.
Nieto, I., Chamorro, D., Palomo, L. J., Real, R. & Muñoz, A.-R. Is the Eurasian Reed Warbler a regular wintering species in the Iberian Peninsula? Ringing data say yes. Acta Ornithol. 53, 61–68 (2018).
Google Scholar 

21.
Ambrosini, R. et al. Climate change and the long-term northward shift in the African wintering range of the barn swallow Hirundo rustica. Clim. Res. 49, 131–141 (2011).
Google Scholar 

22.
Møller, A. P., Rubolini, D. & Lehikoinen, E. Populations of migratory bird species that did not show a phenological response to climate change are declining. Proc. Natl. Acad. Sci. U. S. A. 105, 16195–16200 (2008).
ADS  PubMed  PubMed Central  Google Scholar 

23.
Real, R., Márquez, A. L., Olivero, J. & Estrada, A. Species distribution models in climate change scenarios are still not useful for informing policy planning: An uncertainty assessment using fuzzy logic. Ecography (Cop.) 33, 304–314 (2010).
Google Scholar 

24.
Root, T. L. & Schneider, S. H. Conservation and climate change: The challenges ahead. Conserv. Biol. 20, 706–708 (2006).
PubMed  Google Scholar 

25.
Hahn, S., Bauer, S. & Liechti, F. The natural link between Europe and Africa—2.1 billion birds on migration. Oikos 118, 624–626 (2009).
Google Scholar 

26.
Evans, P. R. & Lathbury, G. W. Raptor migration across the Strait of Gibraltar. Ibis (Lond. 1859). 115, 572–585 (1972).
Google Scholar 

27.
Bijlsma, R. G. Bottleneck Areas for Migratory Birds in the Mediterranean Region: An Assessment of the Problems and Recommendations for Action. (ICBP Study report no. 18, 1987).

28.
Gantenbein, B. & Largiadèr, C. R. The phylogeographic importance of the Strait of Gibraltar as a gene flow barrier in terrestrial arthropods: A case study with the scorpion Buthus occitanus as model organism. Mol. Phylogenet. Evol. 28, 119–130 (2003).
CAS  PubMed  Google Scholar 

29.
García-Mudarra, J. L., Ibáñez, C. & Juste, J. The Straits of Gibraltar: Barrier or bridge to Ibero-Moroccan bat diversity?. Biol. J. Linn. Soc. 96, 434–450 (2009).
Google Scholar 

30.
Carranza, S., Harris, D. J., Arnold, E. N., Batista, V. & de la Gonzalez, V. Phylogeography of the lacertid lizard, Psammodromus algirus, in Iberia and across the Strait of Gibraltar. J. Biogeogr. 33, 1279–1288 (2006).
Google Scholar 

31.
de Juana, E. & Comite de Rarezas de la Sociedad Española de Ornitología, (SEO). Observaciones de Aves raras en España, año 1995. Ardeola44, 119–141 (1997).

32.
de Juana, E. & Comité Ibérico de Rarezas de la Sociedad Española de Ornitología, (SEO). Observaciones homologadas de aves raras en España y Portugal. Informe de 1992. Ardeola41, 103–117 (1994).

33.
Copete, J. L. et al. Observaciones de Aves Raras en España, 2012 y 2013. Ardeola 62, 453–508 (2015).
Google Scholar 

34.
Dies, J. I. et al. Observaciones de Aves Raras en España, 2008. Ardeola 57, 481–516 (2010).
Google Scholar 

35.
Elorriaga, J. & Muñoz, A.-R. Hybridisation between the Common Buzzard Buteo buteo buteo and the North African race of Long-legged Buzzard Buteo rufinus cirtensis in the Strait of Gibraltar: Prelude or preclude to colonisation?. Ostrich J. Afr. Ornithol. 84, 41–45 (2013).
Google Scholar 

36.
Chamorro, D., Olivero, J., Real, R. & Muñoz, A.-R. El cambio del clima y la Barrera Biogeográfica del Estrecho de Gibraltar para las aves africanas. In Avances en biogeografía: Áreas de distribución: entre puentes y barreras (eds. Zotano, J. G., et al.) 172–180 (Tundra Ediciones, 2016).

37.
Ramírez, J. et al. Spring movements of Rüppell’s Vulture Gyps rueppellii across the Strait of Gibraltar. Ostrich 82, 71–73 (2011).
Google Scholar 

38.
del Hoyo, J., Elliott, A. & Sargatal, J. Handbook of the Birds of the World. Vol. 2. New World Vultures to Guineafowl. (Lynx Edicions, Barcelona, 1994).

39.
Ferguson-Lees, J. & Christie, D. A. Raptors of the World (Houghton Mifflin Harcourt, Boston, 2001).
Google Scholar 

40.
Irby, L. H. The Ornithology of the Straits of Gibraltar (Taylor & Francis, Abingdon, 1895).
Google Scholar 

41.
Ávila, A., Guirado, M. A. & Navarrete, J. Observaciones de Aves Raras. Busardo Moro. Ardeola 51, 548 (2004).
Google Scholar 

42.
Cramp, S. & Simmons, K. E. L. Handbook of the Birds of Europe, the Middle East and North Africa: The Birds of the Western Palearctic. Vol. 2: Hawks to Bustards. (Oxford University Press., Oxford, 1980).

43.
Zadeh, L. A. Fuzzy sets. Inf. Control 8, 338–353 (1965).
MATH  Google Scholar 

44.
Real, R. et al. Conservation biogeography of ecologically interacting species: The case of the Iberian lynx and the European rabbit. Divers. Distrib. 15, 390–400 (2009).
Google Scholar 

45.
Olivero, J., Real, R. & Márquez, A. L. Fuzzy chorotypes as a conceptual tool to improve insight into biogeographic patterns. Syst. Biol. 60, 645–660 (2011).
PubMed  Google Scholar 

46.
Olivero, J. et al. Recent loss of closed forests is associated with Ebola virus disease outbreaks. Sci. Rep. 7, 14291 (2017).
ADS  MathSciNet  PubMed  PubMed Central  Google Scholar 

47.
Robertson, M. P., Villet, M. H. & Palmer, A. R. A fuzzy classification technique for predicting species’ distributions: Applications using invasive alien plants and indigenous insects. Divers. Distrib. 10, 461–474 (2004).
Google Scholar 

48.
Salski, A. Ecological applications of fuzzy logic. In Ecological Informatics (ed. Recknagel, F.) 3–14 (Springer, New York, 2006).

49.
Real, R., Barbosa, A. M. & Bull, J. W. Species distributions, quantum theory, and the enhancement of biodiversity measures. Syst. Biol. 66, 453–462 (2017).
PubMed  Google Scholar 

50.
Acevedo, P. & Real, R. Favourability: Concept, distinctive characteristics and potential usefulness. Naturwissenschaften 99, 515–522 (2012).
ADS  CAS  PubMed  Google Scholar 

51.
Hosmer, D. W. & Lemeshow, S. Goodness of fit tests for the multiple logistic regression model. Commun. Stat. Theory Methods 9, 1043–1069 (1980).
MATH  Google Scholar 

52.
Thévenot, M., Vernon, R. & Bergier, P. The Birds of Morocco: An Annotated Checklist. BOU Checklist No. 20. (British Ornithologist’s Union, 2003).

53.
Muñoz, A.-R., Márquez, A. L. & Real, R. An approach to consider behavioral plasticity as a source of uncertainty when forecasting species’ response to climate change. Ecol. Evol. 5, 2359–2373 (2015).
PubMed  PubMed Central  Google Scholar 

54.
Muñoz, A.-R., Real, R., Barbosa, A. M. & Vargas, J. M. Modelling the distribution of Bonelli’s eagle in Spain: Implications for conservation planning. Divers. Distrib. 11, 477–486 (2005).
Google Scholar 

55.
Guisan, A. & Zimmermann, N. E. Predictive habitat distribution models in ecology. Ecol. Modell. 135, 147–186 (2000).
Google Scholar 

56.
Guisan, A. & Thuiller, W. Predicting species distribution: Offering more than simple habitat models. Ecol. Lett. 8, 993–1009 (2005).
Google Scholar 

57.
Gouveia, S. F. et al. Ecophysics reload—exploring applications of theoretical physics in macroecology. Ecol. Modell. 424, 109032 (2020).
Google Scholar 

58.
Cramp, S. Handbook of the Birds of Europe, the Middle East and North Africa: The birds of the Western Paleartic, Vol. VI. (Oxford University Press, Oxford, 1992).

59.
Hardin, G. The competitive exclusion principle. Science (80-). 131, 1292–1297 (1960).
ADS  CAS  Google Scholar 

60.
Jowers, M. J. et al. Unravelling population processes over the Late Pleistocene driving contemporary genetic divergence in Palearctic buzzards. Mol. Phylogenet. Evol. 134, 269–281 (2019).
PubMed  Google Scholar 

61.
Waters, J. M., Fraser, C. I. & Hewitt, G. M. Founder takes all: Density-dependent processes structure biodiversity. Trends Ecol. Evol. 28, 78–85 (2013).
PubMed  Google Scholar 

62.
Gil-Velasco, M. et al. Observaciones de Aves Raras en España, 2017. Ardeola 66, 169–204 (2019).
Google Scholar 

63.
Gutiérrez, R. et al. Observaciones De Aves Raras En España, 2011. Ardeola 60, 437–506 (2013).
Google Scholar 

64.
Gutiérrez, R. et al. Observaciones aves raras en España, 2010. Ardeola 59, 353–411 (2012).
Google Scholar 

65.
Jara, J. et al. Comité Português de Raridades. Anuário Ornitológico da SPEA 7, 140 (2010).
Google Scholar 

66.
Corso, A. Successful mixed breeding of Atlas Long-legged Buzzard and Common Buzzard on Pantelleria, Italy, in 2008. Dutch Bird. 31, 224–226 (2009).
Google Scholar 

67.
Väli, Ü et al. Widespread hybridization between the Greater Spotted Eagle Aquila clanga and the Lesser Spotted Eagle Aquila pomarina (Aves: Accipitriformes) in Europe. Biol. J. Linn. Soc. 100, 725–736 (2010).
Google Scholar 

68.
del Junco, O. & González, B. L. nueva especie de Vencejo en el Paleártico: Apus caffer. Ardeola 13, 115–127 (1969).
Google Scholar 

69.
Ferrero, J. J. Situación del Elanio azul (Elanus caeruleus) en el Mediterráneo. In Biología y Conservación de las Rapaces Mediterráneas,1994 (eds. Muntaner, J. & Mayol, J.) 101–115 (SEO, Sociedad Española de Ornitología, 1996).

70.
Ramírez, J., Simón, M., Solís, S., Pérez, C. & García, E. Vencejo Moro Apus affinis Observaciones de Aves Raras en España. Ardeola 49, 161 (2002).
Google Scholar 

71.
Molina, B., Prieta, J., Lorenzo, J. A. & López-Jurado, C. Noticiario Ornitológico. Ardeola 66, 205–255 (2019).
Google Scholar 

72.
Logeais, J.-M. Première nidification de l’Élanion blanc Elanus caeruleus en Maine-et-Loire. Crex 13, 45–50 (2015).
Google Scholar 

73.
Balbontín, J., Negro, J. J., Sarasola, J. H., Ferrero, J. J. & Rivera, D. Land-use changes may explain the recent range expansion of the Black-shouldered Kite Elanus caeruleus in southern Europe. Ibis (Lond. 1859). 150, 707–716 (2008).
Google Scholar 

74.
Ferrero, J. J. & Onrubia, A. Elanio azul Elanus caeruleus. In Libro Rojo de las Aves de España (eds Madroño, A. et al.) 113–116 (SEO/BirdLife, Madrid, 2004).
Google Scholar 

75.
Knaus, P. et al. Swiss Breeding Bird Atlas 2013–2016. Distribution and Population Trends of Birds in Switzerland and Liechtenstein (Swiss Ornithological Institute, Sempach, 2018).
Google Scholar 

76.
Pulido-Pastor, A., Márquez, A. L., García-Barros, E. & Real, R. Identification of potential source and sink areas for butterflies on the Iberian Peninsula. Insect Conserv. Divers. 11, 479–492 (2018).
Google Scholar 

77.
Muñoz, A.-R. & Real, R. Assessing the potential range expansion of the exotic monk parakeet in Spain. Divers. Distrib. 12, 656–665 (2006).
Google Scholar 

78.
Niamir, A., Skidmore, A. K., Muñoz, A. R., Toxopeus, A. G. & Real, R. Incorporating knowledge uncertainty into species distribution modelling. Biodivers. Conserv. 28, 571–588 (2019).
Google Scholar 

79.
Muñoz, A.-R. & Real, R. Distribution of Bonelli’s Eagle Aquila fasciata in southern Spain: Scale may matter. Acta Ornithol. 48, 93–101 (2013).
Google Scholar 

80.
Barbosa, A. M. & Real, R. Applying fuzzy logic to comparative distribution modelling: A case study with two sympatric amphibians. Sci. World J. 2012, 1–10 (2012).
Google Scholar 

81.
Maguire, B. J. Niche response structure and the analytical potentials of its relationship to the habitat. Am. Nat. 107, 213–246 (1973).
Google Scholar 

82.
Mill, H. The International Geography (D. Appleton and Company, Boston, 1902).
Google Scholar 

83.
Udvardy, M. D. F. A classification of the Biogeographical Provinces of the World. (1975).

84.
Font, I. Climatología de España y Portugal (Ediciones Universidad de Salamanca, Salamanca, 2000).
Google Scholar 

85.
Crovello, T. J. Quantitative Biogeography: An Overview. Taxon 30, 563–575 (1981).
Google Scholar 

86.
Márquez, A. L., Real, R. & Vargas, J. M. Methods for comparison of biotic regionalizations: The case of pteridophytes in the Iberian Peninsula. Ecography (Cop.) 24, 659–670 (2001).
Google Scholar 

87.
Barreau, D., Bergier, P. & Lesne, L. Lávifaune de l’Oukaïmeden, 2200–3600m (Haut Atlas, Maroc). L’Oiseau la R.F.O.57, 307–367 (1987).

88.
Rodriguez, G., Elorriaga, J. & Ramirez, J. Identification of Atlas Long-Legged Buzzard (Buteo rufinus cirtensis) and its status in Europe. Bird. World 26, 147–173 (2013).
Google Scholar 

89.
IUCN. Buteo rufinus. The IUCN Red List of Threatened Species. Version 2018-2. https://www.iucnredlist.org (2017).

90.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).
Google Scholar 

91.
Kirkevåg, A. et al. Aerosol-cloud-climate interactions in the climate model CAM-Oslo. Tellus Ser. A Dyn. Meteorol. Oceanogr. 60A, 492–512 (2008).
Google Scholar 

92.
Collins, W. J. et al. Development and evaluation of an Earth-System model—HadGEM2. Geosci. Model Dev. 4, 1051–1075 (2011).
ADS  Google Scholar 

93.
Fa, J. E. et al. Integrating sustainable hunting in biodiversity protection in central Africa: Hot spots, weak spots, and strong spots. PLoS ONE 9, e112367 (2014).
ADS  PubMed  PubMed Central  Google Scholar 

94.
Benjamini, Y. & Hochberg, Y. Controlling the False Discovery Rate: A practical and powerful approach to Multiple Testing. J. R. Stat. Soc. Ser. B. 57, 289–300 (1995).
MathSciNet  MATH  Google Scholar 

95.
García, L. V. Controlling the false discovery rate in ecological research. Trends Ecol. Evol. 18, 553–554 (2003).
Google Scholar 

96.
Benjamini, Y. & Yekutieli, D. The control of the false discovery rate in multiple testing under depencency. Ann. Stat. 29, 1165–1188 (2001).
MATH  Google Scholar 

97.
Hosmer, D. W. & Lemeshow, S. Applied Logistic Regression (Wiley, Hoboken, 2000).
Google Scholar 

98.
Real, R., Barbosa, A. M. & Vargas, J. M. Obtaining environmental favourability functions from logistic regression. Environ. Ecol. Stat. 13, 237–245 (2006).
MathSciNet  Google Scholar 

99.
Muñoz, A.-R., Jiménez-Valverde, A., Márquez, A. L., Moleón, M. & Real, R. Environmental favourability as a cost-efficient tool to estimate carrying capacity. Divers. Distrib. 21, 1388–1400 (2015).
Google Scholar 

100.
Hosmer, D. W. & Lemeshow, S. Assessing the fit of the model. In Applied Logistic Regression 157–158 (Wiley, Hoboken, 2000).

101.
Chamorro, D., Nieto, I., Real, R. & Muñoz, A. Wintering areas on the move in the face of warmer winters. Ornis Fenn. 96, 1–14 (2019).
Google Scholar 

102.
Wald, A. Tests of statistical hypotheses concerning several parameters when the number of observations is large. Trans. Am. Math. Soc. 54, 426–482 (1943).
MathSciNet  MATH  Google Scholar 

103.
Lobo, J. M., Jiménez-valverde, A. & Real, R. AUC: A misleading measure of the performance of predictive distribution models. Glob. Ecol. Biogeogr. 17, 145–151 (2008).
Google Scholar 

104.
Romero, D., Olivero, J. & Real, R. Comparative assessment of different methods for using land-cover variables for distribution modelling of Salamandra salamandra longirotris. Environ. Conserv. 40, 48–59 (2012).
Google Scholar 

105.
Cohen, J. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 20, 37–46 (1960).
Google Scholar 

106.
Fielding, A. H. & Bell, J. F. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 24, 38–49 (1997).
Google Scholar 

107.
Barbosa, A. M., Real, R., Muñoz, A.-R. & Brown, J. A. New measures for assessing model equilibrium and prediction mismatch in species distribution models. Divers. Distrib. 19, 1333–1338 (2013).
Google Scholar 

108.
Romo, H., García-Barros, E., Márquez, A. L., Moreno, J. C. & Real, R. Effects of climate change on the distribution of ecologically interacting species: Butterflies and their main food plants in Spain. Ecography (Cop.) 37, 1063–1072 (2014).
Google Scholar 

109.
Muñoz, A.-R., Márquez, A. L. & Real, R. Updating known distribution models for forecasting climate change impact on endangered species. PLoS ONE 8, e65462 (2013).
ADS  PubMed  PubMed Central  Google Scholar 

110.
Romero, D., Olivero, J. & Real, R. Accounting for uncertainty in assessing the impact of climate change on biodiversity hotspots in Spain. Anim. Biodivers. Conserv. 42, 355–367 (2019).
Google Scholar 

111.
Daget, P. Ordenation des profils ecologiques. Nat. Monspel. Ser. Bot. 26, 109–128 (1977).
Google Scholar 

112.
Gauch, H. G. Multivariate Analysis in Community Ecology. (Cambridge University Press, Cambridge, 1982). https://doi.org/10.1017/CBO9780511623332.

113.
Real, R., Guerrero, J. C., Antúnez, A., Olivero, J. & Vargas, J. M. Geographic responses of amphibian species to environmental gradients in Southern Spain. I. Individualistic patterns. Boletín la Real Soc. Española Hist. Nat. (Sec. Biol.) 96, 251–261 (2001).
Google Scholar 

114.
Antúnez, A. & Mendonza, M. Factores que determinan el área de distribución geográfica de las especies: Conceptos, modelos y métodos de análisis. In Objetivos y Métodos Biogeográficos. Aplicaciones en Herpetología. Monografías de Herpetología (eds. Vargas, J. M., et al.) vol. 2, 51–72 (Asociación Herpetológica Española, 1992).

115.
Kou, X., Li, Q. & Liu, S. Quantifying species’ range shifts in relation to climate change: A case study of Abies spp. in China. PLoS ONE 6, e23115 (2011).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

116.
US Geological Survey. GTOPO30. (Land processes distributed active archive center (LPDAAC), EROS data center, 1996). More