Philippa Kaur

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

When assessing the ecological status and the development of populations, one important factor to consider is the mortality rate and its underlying causes15. If the status of a population is deemed unsustainable due to high mortality rates, this information can then be used to develop and implement specific management measures, for example addressing the major causes of unnatural mortality16. With regard to the potential ecological relevance of the phenomenon of grey seal predation, it is therefore also important to have reliable estimates of harbour porpoise mortalities resulting from grey seal predation as one natural cause. To allow for a standardised assessment of lesions found in suspected grey seal predation cases, this study aims to summarise the knowledge that has been gathered to date.
The parameters described resemble the most commonly detected lesions in “definite”, “likely” and “fox” related cases of the 246 stranding records categorised as “suspicious” in terms of grey seal predation from the coasts of Schleswig–Holstein. With regard to grey seal predation, parameters 1–9 represent typical lesions, whereas the presence of parameters 10 and 11 is consistent with an interaction with a red fox.
Similar to lesions detected in seals, lesions in porpoises most often resemble puncture lesions in the skin and blubber (parameter 1). Yet, visually most striking is the commonly detected large tissue defect with straight, cut-like wound margins with often flaps of skin and blubber remaining only partly attached to the body (parameters 2, 5, 7). Missing blubber (parameter 4) is also recorded, either as reduced blubber thickness on the flaps of skin or as fully removed parts of blubber and skin. As has been described for seals14, the lesion most often originates in the cervical area (parameter 3). A difference that has been detected between the lesions in seals and porpoises is the rate of clear parallel running bite and / or scratch marks in the skin of the animals. Whereas this is rarely detected in seals14, most porpoises show respective marks (parameter 6). One probable explanation for this dissimilarity could be the different physical and morphological properties of the two types of skin. Seal skin is very dense and elastic; tearing and rupturing the skin requires a considerable amount of force17. Porpoise skin, however, is rather susceptible to applied mechanic force and puncturing or tearing it requires comparably little force18. These different mechanical properties might also be the reason why rake marks are found in the blubber (parameter 9) more often in seals (91% of likely cases14) than in porpoises (62% of likely cases). For seals, in the majority of cases, little to no skin is missing (skin missing in 49% of likely cases14), whereas in porpoises, a considerable number of cases (81% of likely cases) have been found where skin is missing (parameter 7). Grey seals have been described to mainly target the energy rich blubber tissue of their prey11,14. For the elastic and robust seal skin, this is done by scraping off the blubber with the teeth. As porpoise skin is fragile, we suggest that the blubber, including the skin, is more often fully removed by the grey seal and swallowed whole. If true, this may also influence the net energetic gain, which is acquired by the predator. Scraping of blubber tissue from seal skin will likely yield less tissue and cost more energy than tearing off whole pieces of blubber (including skin). Thus, it may result in a lower energetic gain. However, it is still unclear if the process of catching a porpoise in comparison to younger seals might also cost a considerably larger amount of energy, negatively influencing the net gain.
Similar to what has been described for seals, the avulsion of one or both scapulae (parameter 8) can be found and is also likely the result of the force applied when detaching the epidermis and blubber from the body of the prey14.
For porpoises, all nine suggested parameters were found in the definite case of grey seal predation. Parameters 1–5 showed a very high (100%) and parameter 6 a high rate (95%) of occurrence in likely cases. Parameters 7–9 occurred less frequently but were still found in > 60% of all likely predation cases. These high rates of occurrence throughout all parameters suggest that wound patterns found in porpoises are less variable than the patterns found in seals14. Whether this difference is a result of the different mechanical skin properties or if other factors are responsible, is beyond the scope of this study.
While for seals a skeletal trauma is used as an indicator of grey seal predation, for respective harbour porpoise cases, this is hardly ever (19% of likely cases) documented. In contrast, for porpoises, a skeletal trauma (parameter 10) is quite frequently detected in cases related to scavenging by foxes (46% of fox cases) where for example extremities can be manipulated19. As has been reported in seals14, the most often detected parameter in fox related cases is the ragged wound margin (parameter 11, 94% of the cases). Therefore, this can be seen as a good indicator of an interaction with a fox in porpoises. This is also supported by a definite case of fox scavenging, which was confirmed using genetic methods20. It needs to be stated though, that scavenging by birds can result in similar looking lesions, increasing the chance of misinterpretation. Scavenging by birds usually also leaves an irregular wound margin with extensive tissue loss. If parallel running lesions are present, the origin of the lesion can additionally be assessed by measuring the distances in between adjacent lesions and comparing them to published values of grey seal, fox and cetacean inter-teeth distances e.g.1. This is especially important when differentiating between for example rake marks by dolphins, which have been documented in porpoises21,22 and marks induced by grey seals. Here, it can be useful to assess the pattern of inter-teeth distances with those of dolphins expected to be consistent in length, whereas for grey seals, variable distances are expected as the result of the polydont dental morphology23. Despite a lack of available data, a differentiation between grey seal claw-induced marks and dolphin rake marks should be possible, as distances between claws of a subadult / adult grey seal male are expected to be considerably larger than for dolphin inter-teeth distances.
Single puncture lesions, in turn, are not considered as a very good indicator despite being present in the definite and all of the likely cases. Mainly due to the susceptibility of the porpoise skin, such lesions can have many different causes (e.g. feeding by birds).
Whether a loss of muscle tissue can be attributed to grey seal predation or is largely caused by scavengers like gulls as has been suggested for seals14, is still not entirely clear. In German as well as bordering waters, no clear pattern prevails. Carcasses with mainly intact as well as fully removed muscle tissue have been documented c.f.13. However, the reports by Stringell et al.4 suggest that not only the blubber tissue is targeted, but that there may also be some individual behavioural variation.
The findings and the resulting parameters described here are in line with wound patterns reported in earlier publications from other areas1,6,7,13. This shows that the documented wound patterns make a reliable set of parameters when assessing harbour porpoises carcasses potentially predated by a grey seal and should be used in future assessments.
As a complementary tool to the suggested parameters, corresponding to porpoises, we developed a decision tree with the aim of supporting a standardised and information-based decision-making process. Despite an accuracy in decision-making of 96% when using our data set, the example in Fig. 5 illustrates the limitations of such static tools when it comes to judging more complex cases. Furthermore, when comparing the suggestion given by the tree with the one made by the experts, in only 50% of unmatched cases, a rather precautionary judgement was made, bearing the risk of an overestimation of case numbers. Therefore, we recommend using the suggested tree only as an informational tool in supporting decision-making and final judgments should always be made by the responsible expert based on all available information.
In addition to cases for which the attack of a grey seal directly led to the death of the animal, interestingly, it seems not unusual that porpoises escape this predator. Several observations have been described in the literature5,6,13,24 and nine cases were documented in German waters (Figs. 1, 2). In order to be able to verify the origin of recorded teeth marks in porpoise skin, it is crucial to record marks in detail including their pattern, location and inter-teeth distances. Using the latter, for example, interactions with dolphins can potentially be excluded. Although there has been the odd case of a severely injured seal showing comparable lesions to what is associated with grey seal predation14, such high rates of escape cases as described for porpoises have not been reported.
Despite the co-occurrence of porpoises and grey seals in the Baltic, no case of grey seal predation on a porpoise can be confirmed by the presented results. It remains unclear whether grey seals in this area of the Baltic just don’t prey on porpoises or whether other factors like differences in behaviour (e.g. primary area of predation further offshore) are involved.
Some of the observed behaviour of grey seals when catching a porpoise can be directly linked to the detected lesions. For example, Stringell et al.4 as well as Bouveroux et al.7 described the grey seal acting as an ambush predator and attacking the porpoise from below using its jaws to catch and retain the prey. Lesions starting in the throat area (parameter 3) combined with parallel multifocal puncture lesions (parameter 6) resemble what would be expected as the result of such an attack.
Despite a lower rate of variability in detected wound patterns in porpoise carcasses, care should be applied when assessing lesions, as there is always the chance of other factors being involved. Therefore, if possible, a combination of data sources (necropsy results, genetic detection of predator DNA, indicators at the stranding site, eye witness reports, etc.) should be used in a systematic evaluation.
Future research should focus on continuing thorough investigations of stranded marine mammal carcasses in order to further update and refine the suggested set of parameters. Additionally, results of current as well as retrospective analysis of stranding data should be used to support an evaluation of the ecological relevance of this behaviour. More

1.
Mucina, L. & Rutherford, M. C. The Vegetation of South Africa, Lesotho and Swaziland (South African National Biodiversity Institute, Pretoria, 2006).
Google Scholar 
2.
van Zinderen Bakker, E. M. The evolution of late Quaternary paleoclimates of Southern Africa. Palaeoecol. Afr. 9, 160–202 (1976).
Google Scholar 

3.
Cockcroft, M. J., Wilkinson, M. J. & Tyson, P. D. The application of a present-day climatic model to the late Quaternary in southern Africa. Clim. Change 10, 161–181 (1987).
ADS  Google Scholar 

4.
Chase, B. M. & Meadows, M. E. Late Quaternary dynamics of southern Africa’s winter rainfall zone. Earth Sci. Rev. 84(3), 103–138 (2007).
ADS  Google Scholar 

5.
Bistinas, I., Harrison, S. P., Prentice, I. C. & Pereira, J. M. C. Causal relationships vs. emergent patterns in the global controls of fire frequency. Biogeosciences 11, 5087–5101 (2014).
ADS  Google Scholar 

6.
Hoetzel, S., Dupont, L., Schefuß, E., Rommerskirchen, F. & Wefer, G. The role of fire in Miocene to Pliocene C 4 grassland and ecosystem evolution. Nat. Geosci. 6(12), 1027–1030 (2013).
ADS  CAS  Google Scholar 

7.
Bond, W. J., Woodward, F. I. & Midgley, G. F. The global distribution of ecosystems in a world without fire. New Phytol. 165(2), 525–538 (2005).
CAS  PubMed  Google Scholar 

8.
Ripley, B. et al. Fire ecology of C3 and C4 grasses depends on evolutionary history and frequency of burning but not photosynthetic type. Ecology 96(10), 2679–2691 (2015).
PubMed  Google Scholar 

9.
Pinto, H., Sharwood, R. E., Tissue, D. T. & Ghannoum, O. Photosynthesis of C3, C3–C4, and C4 grasses at glacial CO2. J. Exp. Bot. 65(13), 3669–3681 (2014).
PubMed  PubMed Central  Google Scholar 

10.
Roth-Nebelsick, A. & Konrad, W. Habitat responses of fossil plant species to palaeoclimate—possible interference with CO2?. Palaeogeogr. Palaeoclimatol. Palaeoecol. 467, 277–286 (2017).
Google Scholar 

11.
Ehleringer, J. R., Cerling, T. E. & Helliker, B. R. C4 photosynthesis, atmospheric CO2, and climate. Oecologia 112(3), 285–299 (1997).
ADS  PubMed  Google Scholar 

12.
Edwards, E. J., Osborne, C. P., Strömberg, C. A., Smith, S. A. & C4 Grasses Consortium. The origins of C4 grasslands: integrating evolutionary and ecosystem science. Science 328(5978), 587–591 (2010).
CAS  PubMed  Google Scholar 

13.
Hönisch, B., Hemming, N. G., Archer, D., Siddall, M. & McManus, J. F. Atmospheric carbondioxide concentration across the mid-Pleistocene transition. Science 324(5934), 1551–1554 (2009).
ADS  PubMed  Google Scholar 

14.
Yan, Y. et al. Two-million-year-old snapshots of atmospheric gases from Antarctic ice. Nature 574(7780), 663–666 (2019).
ADS  CAS  PubMed  Google Scholar 

15.
Faith, J. T., Rowan, J. & Du, A. Early hominins evolved within non-analog ecosystems. Proc. Natl. Acad. Sci. 116(43), 21478–21483 (2019).
ADS  CAS  PubMed  Google Scholar 

16.
Sealy, J., Naidoo, N., Hare, V. J., Brunton, S. & Faith, J. T. Climate and ecology of the palaeo-Agulhas Plain from stable carbon and oxygen isotopes in bovid tooth enamel from Nelson Bay Cave, South Africa. Quat. Sci. Rev. 235, 105974 (2019).
Google Scholar 

17.
Horwitz, L. K. & Chazan, M. Past and present at Wonderwerk Cave (Northern Cape Province, South Africa). Afr. Archaeol. Rev. 32(4), 595–612 (2015).
Google Scholar 

18.
Ecker, M. et al. The palaeoecological context of the Oldowan-Acheulean in southern Africa. Nat. Ecol. Evol. 2(7), 1080–1086 (2018).
PubMed  Google Scholar 

19.
Matmon, A. et al. New chronology for the southern Kalahari Group sediments with implications for sediment-cycle dynamics and early hominin occupation. Quat. Res. 84(1), 118–132 (2015).
Google Scholar 

20.
Vainer, S., Erel, Y. & Matmon, A. Provenance and depositional environments of Quaternary sediments in the southern Kalahari Basin. Chem. Geol. 476, 352–369 (2018).
ADS  CAS  Google Scholar 

21.
Prentice, I. C. et al. Modeling fire and the terrestrial carbon balance. Glob. Biogeochem. Cycles 25(3), 2–13 (2011).
Google Scholar 

22.
Braconnot, P. et al. Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum-Part 1: experiments and large-scale features. Clim. Past 3(2), 261–277 (2007).
Google Scholar 

23.
Kelley, D. I. et al. A comprehensive benchmarking system for evaluating global vegetation models. Biogeosciences 10(5), 3313–3340 (2013).
ADS  Google Scholar 

24.
Chazan, M. et al. Archaeology, paleoenvironment and chronology of the early middle stone age component of Wonderwerk cave in the interior of southern Africa. J. Palaeolithic Archaeol. https://doi.org/10.1007/s41982-020-00051-8 (2020).
Article  Google Scholar 

25.
Lee-Thorp, J. A. & Beaumont, P. B. Vegetation and seasonality shifts during the late Quaternary deduced from 13C/12C ratios of grazers at Equus Cave, South Africa. Quat. Res. 43, 426–432 (1995).
Google Scholar 

26.
Vogel, J. C. The geographical distribution of Kranz species in southern Africa. South Afr. J. Sci. 75, 209–215 (1978).
Google Scholar 

27.
Zhou, H., Helliker, B. R., Huber, M., Dicks, A. & Akçay, E. C4 photosynthesis and climate through the lens of optimality. Proc. Natl. Acad. Sci. 115(47), 12057–12062 (2018).
CAS  PubMed  Google Scholar 

28.
Rubin, F., Palmer, A. R. & Tyson, C. Patterns of endemism within the Karoo National Park, South Africa. Bothalia 31(1), 117–133 (2001).
Google Scholar 

29.
Walker, S. J., Lukich, V. & Chazan, M. Kathu townlands: a high density earlier stone age locality in the interior of South Africa. PLoS ONE 9(7), e103436 (2014).
ADS  PubMed  PubMed Central  Google Scholar 

30.
Lee-Thorp, J. A., Sponheimer, M. & Luyt, J. Tracking changing environments using stable carbon isotopes in fossil tooth enamel: an example from the South African hominin sites. J. Hum. Evol. 53(5), 595–601 (2007).
PubMed  Google Scholar 

31.
Codron, D., Brink, J. S., Rossouw, L. & Clauss, M. The evolution of ecological specialization in southern African ungulates: competition- or physical environmental turnover?. Oikos 117, 344–353 (2008).
Google Scholar 

32.
Plummer, T. W. et al. The environmental context of Oldowan hominin activities at Kanjera South, Kenya. In Interdisciplinary approaches to the Oldowan (eds Hovers, E. & Braun, D. R.) 149–160 (Springer, Berlin, 2009).
Google Scholar 

33.
Cerling, T. E. et al. Dietary changes of large herbivores in the Turkana Basin, Kenya from 4 to 1 Ma. Proc. Natl. Acad. Sci. 112(37), 11467–11472 (2015).
ADS  CAS  PubMed  Google Scholar 

34.
Harris, I., Osborn, T. J., Jones, P. & Lister, D. Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. Sci. Data. https://doi.org/10.1038/s41597-020-0453-3 (2020).
Article  PubMed  PubMed Central  Google Scholar 

35.
Sitch, S. et al. Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model. Glob. Change Biol. 9(2), 161–185 (2003).
ADS  Google Scholar 

36.
Thonicke, K. et al. The influence of vegetation, fire spread and fire behaviour on biomass burning and trace gas emissions: results from a process-based model. Biogeosciences 7(6), 1991–2011 (2010).
ADS  CAS  Google Scholar 

37.
Haxeltine, A. & Prentice, I. C. BIOME3: an equilibrium terrestrial biosphere model based on ecophysiological constraints, resource availability, and competition among plant functional types. Glob. Biogeochem. Cycles 10(4), 693–709 (1996).
ADS  CAS  Google Scholar 

38.
Haxeltine, A. & Prentice, I. C. A general model for the light-use efficiency of primary production. Funct. Ecol. 10, 551–561 (1996).
Google Scholar 

39.
Farquhar, G. D., Von Caemmerer, S. & Berry, J. A. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 plants. Planta 149, 78–90 (1980).
CAS  PubMed  PubMed Central  Google Scholar 

40.
Farquhar, G. D. & Von Caemmerer, S. Modelling of photosynthetic response to environmental conditions. In Physiological Plant Ecology II: Water Relations and Carbon Assimilation (eds Nobel, P. S. et al.) 549–587 (Springer, Berlin, 1982).
Google Scholar 

41.
Monteith, J. L. A reinterpretation of stomatal responses to humidity. Plant Cell Environ. 18, 357–364 (1995).
Google Scholar 

42.
Rothermel, R. C. A Mathematical Model for Predicting Fire Spread in Wildland Fuels (Vol. 115). Intermountain Forest and Range Experiment Station, Forest Service, US Department of Agriculture (1972).

43.
Sato, H., Kelley, D. I., Mayor, S. J., Cowling S. A., Calvo, M. M. & Prentice, I. C. Fire and low CO2 opened dry corridors in South America during the Last Glacial Maximum. Under Review for Nature Geosciences: NGS-2019–07–01558B (2020).

44.
Prentice, I. C., Harrison, S. P. & Bartlein, P. J. Global vegetation and terrestrial carbon cycle changes after the last ice age. New Phytol. 189(4), 988–998 (2011).
CAS  PubMed  Google Scholar  More

Cohort and data description
In order to examine seasonal changes of human molecular data, we leveraged the power of longitudinal multiomics data from profiling of 105 individuals (55 women and 50 men) with ages ranging from 25 to 75 years old (Fig. 1a; Supplementary Table 1). This cohort was generally healthy and well characterized for glucose dysregulation using annual oral glucose tolerance tests (OGTTs), insulin resistance measuring steady-state plasma glucose (SSPG), fasting glucose and hemoglobin A1c (HbA1c; an indicator of the average level of blood glucose over the past 100 days)19 as well as quarterly sample collections with measurements of transcriptomes (from peripheral blood mononuclear cells), proteome and metabolome from plasma, targeted cytokine and growth factor assays using serum. Nasal and gut microbiomes were analyzed using 16S rRNA sequencing providing information at the genus level and host exome sequencing was performed once from PBMCs (Fig. 1b). Moreover, 51 clinical laboratory tests were acquired on each visit and they were aligned to the meteorological data (e.g. air temperature), pollen counts (e.g. mold spores, grass pollens, tree pollens, weed pollens) and airborne fungi from the San Francisco bay area. In total, there were 902 visits (average across different types of omes‘) from which samples were drawn over up to 4 years (see “Methods”). The sample collections were generally evenly distributed throughout the year (Fig. 1b). Nearly all individuals lived in the San Francisco Bay Area with the exception of three individuals who lived in Southern California and frequented the Bay area (Supplementary Data 1). Participants in our study were well characterized for steady-state plasma glucose (SSPG) using the modified insulin suppression test20, in which 31 participants were insulin sensitive (SSPG  0.05, Supplementary Table 5, Supplementary Fig. 10). In our analysis we used subject ID as a random effect to account for different numbers of samples per subject. On the other hand, physical activity measured in total metabolic equivalent of task (MET) is significantly different between the IR and the IS groups in February, May, June, and August (P-value = 0.01787, Supplementary Fig. 11). However, a post-hoc analysis of all the omics features that were identified to be significantly different between the IR and the IS groups, are not associated with the physical activity. More

1.
United Nations, Department of Economic and Social Affairs, Population Division. World Population Prospects 2019—Data Booklet (ST/ESA/ SER.A/377), (2019). https://population.un.org/wpp/Publications/Files/WPP2019_DataBooklet.pdf
2.
Pimentel, D. & Burgess, M. Environmental and economic costs of the application of pesticides primarily in the United States. In Integrated Pest Management: Innovation-Development Process (eds Peshin, R. & Dhawan, A. K.) 47–71 (Springer, Dordrecht, 2014). https://doi.org/10.1007/978-1-4020-8992-3_4
Google Scholar 

3.
Devine, G. J. & Furlong, M. J. Insecticide use: Contexts and ecological consequences. Agric. Hum. Values 24(3), 281–306. https://doi.org/10.1007/s10460-007-9067-z (2007).
Article  Google Scholar 

4.
Sanchis, V. & Bourguet, D. Bacillus thuringiensis: Applications in agriculture and insect resistance management. A review. Agron. Sustain. Dev. 28(1), 11–20. https://doi.org/10.1051/agro:2007054 (2008).
Article  Google Scholar 

5.
WHO report. WHO specifications and evaluations for public health pesticides: Bacillus thuringiensis subspecies israelensis strain AM65-52. (World Health Organization, Geneva, 2007).

6.
Rizzati, V., Briand, O., Guillou, H. & Gamet-Payrastre, L. Effects of pesticide mixtures in human and animal models: An update of the recent literature. Chem. Biol. Interact. 254, 231–246. https://doi.org/10.1016/j.cbi.2016.06.003 (2016).
Article  PubMed  CAS  Google Scholar 

7.
Lacey, L. A. et al. Insect pathogens as biological control agents: Back to the future. J. Invertebr. Pathol. 132, 1–41. https://doi.org/10.1016/j.jip.2015.07.009 (2015).
Article  PubMed  CAS  Google Scholar 

8.
Adang, M. J., Crickmore, N. & Jurat-Fuentes, J. L. Diversity of Bacillus thuringiensis Crystal Toxins and Mechanism of Action. Adv. Insect Physiol. 47, 39–87. https://doi.org/10.1016/B978-0-12-800197-4.00002-6 (2014).
Article  Google Scholar 

9.
Crickmore, N. Bacillus thuringiensis toxin classification. In Bacillus thuringiensis and Lysinibacillus sphaericus. (eds Fiuza, L.M. et al.) ISBN 978-3-319-56677-1, 41-52, (Spinger, Cham, 2017).

10.
WHO report. Guideline specification for bacterial larvicides for public health use. WHO document WHO/CDS/CPC/WHOPES/99.2 (World Health Organization, Geneva, 1999).

11.
Bravo, A., Pacheco, S., Gomez, I., Garcia-Gomez B., Onofre, J., Soberon, M. Insecticidal Proteins from Bacillus thuringiensis and their Mechanism of Action. In Bacillus thuringiensis and Lysinibacillus sphaericus (eds Fiuza, L.M. et al.) ISBN 978-3-319-56677-1, 53–66, (Spinger, Cham, 2017).

12.
Palma, L., Muñoz, D., Berry, C., Murillo, J. & Caballero, P. Bacillus thuringiensis toxins: An overview of their biocidal activity. Toxins 6(12), 3296–3325. https://doi.org/10.3390/toxins6123296 (2014).
Article  PubMed  PubMed Central  CAS  Google Scholar 

13.
Ben-Dov, E. et al. Extended screening by PCR for seven cry-group genes from field-collected strains of Bacillus thuringiensis. Appl. Environ. Microb. 63(12), 4883–4890. https://doi.org/10.1128/aem.63.12.4883-4890.1997 (1997).
CAS  Google Scholar 

14.
Berry, C. et al. Complete sequence and organization of pBtoxis, the toxin-coding plasmid of Bacillus thuringiensis subsp. israelensis. Appl. Environ. Microbiol. 68(10), 5082–5095. https://doi.org/10.1128/aem.68.10.5082-5095.2002 (2002).
Article  PubMed  PubMed Central  CAS  Google Scholar 

15.
Bravo, A., Gill, S. S. & Soberon, M. Mode of action of Bacillus thuringiensis Cry and Cyt toxins and their potential for insect control. Toxicon 49, 423–435. https://doi.org/10.1016/j.toxicon.2006.11.022 (2007).
Article  PubMed  CAS  Google Scholar 

16.
Wei, J. et al. Activation of Bt protoxin Cry1Ac in resistant and susceptible cotton bollworm. PLoS ONE 11(6), e0156560. https://doi.org/10.1371/journal.pone.0156560 (2016).
Article  PubMed  PubMed Central  CAS  Google Scholar 

17.
Bravo, A., Likitvivatanavong, S., Gill, S. S. & Soberon, M. Bacillus thuringiensis: A story of a successful bioinsecticide. Insect Biochem. Mol. Biol. 41(7), 423–431. https://doi.org/10.1016/j.ibmb.2011.02.006 (2011).
Article  PubMed  PubMed Central  CAS  Google Scholar 

18.
Caccia, S. et al. Midgut microbiota and host immunocompetence underlie Bacillus thuringiensis killing mechanism. Proc. Natl. Acad. Sci. USA 113(34), 9486–9491. https://doi.org/10.1073/pnas.1521741113 (2016).
Article  PubMed  CAS  Google Scholar 

19.
Glare, T.R., O’Callaghan, M. Bacillus thuringiensis: Biology, Ecology and Safety. ISBN: 9780471496304, 350, (Wiley, New York, 2000).

20.
Rubio-Infante, N. & Moreno-Fierros, L. An overview of the safety and biological effects of Bacillus thuringiensis Cry toxins in mammals. J. Appl. Toxicol. 36, 630–648. https://doi.org/10.1002/jat.3252 (2016).
Article  PubMed  CAS  Google Scholar 

21.
EFSA Panel on Biological Hazards (BIOHAZ). Risks for public health related to the presence of Bacillus cereus and other Bacillus spp. including Bacillus thuringiensis in foodstuffs. EFSA J. https://doi.org/10.2903/j.efsa.2016.4524 (2016).
Article  Google Scholar 

22.
Amichot, M., Curty, C., Benguettat-Magliano, O., Gallet, A. & Wajnberg, E. Side effects of Bacillus thuringiensis var. kurstaki on the hymenopterous parasitic wasp Trichogramma chilonis. Environ. Sci. Pollut. Res. Int. 23, 3097–3103. https://doi.org/10.1007/s11356-015-5830-7 (2016).
Article  PubMed  CAS  Google Scholar 

23.
Renzi, M. T. et al. Chronic toxicity and physiological changes induced in the honey bee by the exposure to fipronil and Bacillus thuringiensis spores alone or combined. Ecotoxicol. Environ. Saf. 127, 205–213. https://doi.org/10.1016/j.ecoenv.2016.01.028 (2016).
Article  PubMed  CAS  Google Scholar 

24.
Caquet, T., Roucaute, M., Le Goff, P. & Lagadic, L. Effects of repeated field applications of two formulations of Bacillus thuringiensis var. israelensis on non-target saltmarsh invertebrates in Atlantic coastal wetlands. Ecotoxicol. Environ. Saf. 74, 1122–1130. https://doi.org/10.1016/j.ecoenv.2011.04.028 (2011).
Article  PubMed  CAS  Google Scholar 

25.
Duguma, D. et al. Microbial communities and nutrient dynamics in experimental microcosms are altered after the application of a high dose of Bti. J. Appl. Ecol. 52, 763–773. https://doi.org/10.1111/1365-2664.12422 (2015).
Article  CAS  Google Scholar 

26.
Venter, H. J. & Bøhn, T. Interactions between Bt crops and aquatic ecosystems: A review. Environ. Toxicol. Chem. 35(12), 2891–2902. https://doi.org/10.1002/etc.3583 (2016).
Article  PubMed  CAS  Google Scholar 

27.
van Frankenhuyzen, K. Specificity and cross-order activity of Bacillus thuringiensis pesticidal proteins. In Bacillus thuringiensis and Lysinibacillus sphaericus (eds Fiuza, L.M. et al.) ISBN 978-3-319-56677-1, 127–172, (Springer, Cham, 2017).

28.
Bizzarri, M. F. & Bishop, A. H. The ecology of Bacillus thuringiensis on the phylloplane: Colonization from soil, plasmid transfer, and interaction with larvae of Pieris brassicae. Microb. Ecol. 56(1), 133–139. https://doi.org/10.1007/s00248-007-9331-1 (2008).
Article  PubMed  CAS  Google Scholar 

29.
Raymond, B., Wyres, K. L., Sheppard, S. K., Ellis, R. J. & Bonsall, M. B. Environmental factors determining the epidemiology and population genetic structure of the Bacillus cereus group in the field. PLoS Pathog. 6(5), e1000905. https://doi.org/10.1371/journal.ppat.1000905 (2010).
Article  PubMed  PubMed Central  CAS  Google Scholar 

30.
Hendriksen, N. B. & Hansen, B. M. Long-term survival and germination of Bacillus thuringiensis var. kurstaki in a field trial. Can. J. Microbiol. 48(3), 256–261. https://doi.org/10.1139/w02-009 (2002).
Article  PubMed  CAS  Google Scholar 

31.
Hung, T. P. et al. Persistence of detectable insecticidal proteins from Bacillus thuringiensis (Cry) and toxicity after adsorption on contrasting soils. Environ. Pollut. 208, 318–325. https://doi.org/10.1016/j.envpol.2015.09.046 (2016).
Article  PubMed  CAS  Google Scholar 

32.
Hung, T. P. et al. Fate of insecticidal Bacillus thuringiensis Cry protein in soil: Differences between purified toxin and biopesticide formulation. Pest Manag. Sci. 72, 2247–2253. https://doi.org/10.1002/ps.4262 (2016).
Article  PubMed  CAS  Google Scholar 

33.
Enger, K. S. et al. Evaluating the long-term persistence of Bacillus spores on common surfaces. Microb. Biotechnol. 11(6), 1048–1059. https://doi.org/10.1111/1751-7915.13267 (2018).
Article  PubMed  PubMed Central  CAS  Google Scholar 

34.
Couch, T.L. Industrial fermentation and formulation of entomopathogenic bacteria. In Entomopathogenic Bacteria: From Laboratory to Field Application (eds Charles, J.-F. et al.) ISBN 978-90-481-5542-2, 297–316.43, (Springer, Dordrecht, 2000).

35.
Brar, S. K., Verma, M., Tyagi, R. D. & Valéro, J. R. Recent advances in downstream processing and formulations of Bacillus thuringiensis based biopesticides. Process Biochem. 41(2), 323–342. https://doi.org/10.1016/j.procbio.2005.07.015 (2006).
Article  CAS  Google Scholar 

36.
Setlow, P. Spore resistance properties. Microbiol. Spectr. 2(5), TBS-0003-2012. https://doi.org/10.1128/microbiolspec.TBS-0003-2012 (2014).
Article  CAS  Google Scholar 

37.
European Food Safety Authority. Conclusion on the peer review of the pesticide risk assessment of the active substance Bacillus thuringiensis subsp. Kurstaki (strains ABTS 351, PB 54, SA 11, SA 12, EG 2348). EFSA J. 10(2), 2540. https://doi.org/10.2903/j.efsa.2012.2540 (2012).
Article  CAS  Google Scholar 

38.
Bächli, G. TaxoDros: The database on Taxonomy of Drosophilidae: Database 2020/1.https://www.taxodros.uzh.ch. (1999–2020).

39.
Tennessen, J. M. & Thummel, C. S. Coordinating growth and maturation—Insights from Drosophila. Curr. Biol. 21(18), R750–R757. https://doi.org/10.1016/j.cub.2011.06.033 (2011).
Article  PubMed  PubMed Central  CAS  Google Scholar 

40.
Benz, G. & Perron, J. M. The toxic action of Bacillus thuringiensis “exotoxin” on Drosophila reared in yeast-containing and yeast-free media. Experientia 23(10), 871–872 (1967).
PubMed  CAS  Google Scholar 

41.
Saadoun, I., Al-Moman, F., Obeidat, M., Meqdam, M. & Elbetieha, A. Assessment of toxic potential of local Jordanian Bacillus thuringiensis strains on Drosophila melanogaster and Culex sp. (Diptera). J. Appl. Microbiol. 90, 866–872. https://doi.org/10.1046/j.1365-2672.2001.01315.x (2001).
Article  PubMed  CAS  Google Scholar 

42.
Khyami-Horani, H. Toxicity of Bacillus thuringiensis and B. sphaericus to laboratory populations of Drosophila melanogaster (Diptera: Drosophilidae). J. Basic Microbiol. 42(2), 105–110. https://doi.org/10.1002/1521-4028(200205)42:23.0.CO;2-S (2002). 
Article  PubMed  Google Scholar 

43.
Obeidat, M. Toxicity of local Bacillus thuringiensis isolates against Drosophila melanogaster. WJAS 4(2), 161–167 (2008).
Google Scholar 

44.
Obeidat, M., Khymani-Horani, H. & Al-Momani, F. Toxicity of Bacillus thuringiensis β-exotoxins and δ-endotoxins to Drosophila melanogaster, Ephestia kuhniella and human erythrocytes. Afr. J. Biotechnol. 11(46), 10504–10512 (2012).
Google Scholar 

45.
Cossentine, J., Robertson, M. & Xu, D. Biological activity of Bacillus thuringiensis in Drosophila suzukii (Diptera: Drosophilidae). J. Econ. Entomol. 109(3), 1–8. https://doi.org/10.1093/jee/tow062 (2016).
Article  CAS  Google Scholar 

46.
Biganski, S., Jehle, J. A. & Kleepies, R. G. Bacillus thuringiensis serovar israelensis has no effect on Drosophila suzukii Matsumura. J. Appl. Entomol. 142, 33–36. https://doi.org/10.1111/jen.12415 (2017).
Google Scholar 

47.
Haller, S., Romeis, J. X. R. & Meissle, M. Effects of purified or plant-produced Cry proteins on Drosophila melanogaster (Diptera: Drosophilidae) larvae. Sci. Rep. 7(1), 11172. https://doi.org/10.1038/s41598-017-10801-4 (2017).
ADS  Article  PubMed  PubMed Central  CAS  Google Scholar 

48.
Benado, M. & Brncic, D. An eight-year phenological study of a local drosophilid community in Central Chile. J. Zool. Syst. Evol. Res. 32, 51–63. https://doi.org/10.1111/j.1439-0469.1994.tb00470.x (1994).
Article  Google Scholar 

49.
Nunney, L. The colonization of oranges by the cosmopolitan Drosophila. Oecologia 108, 552–561. https://www.jstor.org/stable/4221451 (1996).
ADS  PubMed  Google Scholar 

50.
Mitsui, H. & Kimura, M. T. Coexistence of drosophilid flies: Aggregation, patch size diversity and parasitism. Ecol. Res. 15, 93–100.  https://doi.org/10.1046/j.1440-1703.2000.00328.x (2000).
Google Scholar 

51.
Withers, P. & Allemand, R. Les drosophiles de la région Rhône-Alpes (Diptera, Drosophilidae). Bull. Soc. Entomol. Fr. 117(4), 473–482. https://www.persee.fr/doc/bsef_0037-928x_2012_num_117_4_3076 (2012).
Google Scholar 

52.
Stevens, T., Song, S., Bruning, J. B., Choo, A. & Baxter, S. W. Expressing a moth abcc2 gene in transgenic Drosophila causes susceptibility to Bt Cry1Ac without requiring a cadherin-like protein receptor. Insect Biochem. Mol. Biol. 80, 61–70. https://doi.org/10.1016/j.ibmb.2016.11.008 (2017).
Article  PubMed  CAS  Google Scholar 

53.
George, Z., Crickmore, N. Bacillus thuringiensis applications in agriculture. In Bacillus thuringiensis Biotechnology (ed Sansinenea, E.) 392, (Springer, Dordrecht, 2012).

54.
Nepoux, V., Haag, C. R. & Kawecki, T. J. Effects of inbreeding on aversive learning in Drosophila. J. Evol. Biol. 23, 2333–2345. https://doi.org/10.1111/j.1420-9101.2010.02094.x (2010).
Article  PubMed  CAS  Google Scholar 

55.
Vantaux, A., Ouattarra, I., Lefèvre, T. & Dabiré, K. R. Effects of larvicidal and larval nutritional stresses on Anopheles gambiae development, survival and competence for Plasmodium falciparum. Parasite. Vector. 9, 226. https://doi.org/10.1186/s13071-016-1514-5 (2016).
Article  CAS  Google Scholar 

56.
Moret, Y. & Schmid-Hempel, P. Survival for immunity: The price of immune system activation for bumblebee workers. Science 290(5494), 1166–1168. https://doi.org/10.1126/science.290.5494.1166 (2000).
ADS  Article  PubMed  CAS  Google Scholar 

57.
Kutzer, M. A. & Armitage, S. A. O. The effect of diet and time after bacterial infection on fecundity, resistance, and tolerance in Drosophila melanogaster. Ecol. Evol. 6(13), 4229–4242. https://doi.org/10.1002/ece3.2185 (2016).
Article  PubMed  PubMed Central  Google Scholar 

58.
Andersen, L. H., Kristensen, T. N., Loeschcke, V., Toft, S. & Mayntz, D. Protein and carbohydrate composition of larval food affects tolerance to thermal stress and desiccation in adult Drosophila melanogaster. J. Insect Physiol. 56, 336–340. https://doi.org/10.1016/j.jinsphys.2009.11.006 (2010).
Article  PubMed  CAS  Google Scholar 

59.
Rion, S. & Kawecki, T. J. Evolutionary biology of starvation resistance: What we have learned from Drosophila. J. Evol. Biol. 20(5), 1655–1664. https://doi.org/10.1111/j.1420-9101.2007.01405.x (2007).
Article  PubMed  CAS  Google Scholar 

60.
Burger, J. M. S., Buechel, S. D. & Kawecki, T. J. Dietary restriction affects lifespan but not cognitive aging in Drosophila melanogaster. Aging Cell 9, 327–335. https://doi.org/10.1111/j.1474-9726.2010.00560.x (2010).
Article  PubMed  CAS  Google Scholar 

61.
Khazaeli, A. A. & Curtsinger, J. W. Genetic analysis of extended lifespan in Drosophila melanogaster III. On the relationship between artificially selected and wild stocks. Genetica 109, 245–253. https://doi.org/10.1023/a:1017569318401 (2000).
Article  PubMed  CAS  Google Scholar 

62.
Atkinson, W. & Shorrocks, B. Breeding site specificity in the domestic species of Drosophila. Oecologia 29(3), 223–232. https://www.jstor.org/stable/4215461 (1977).
ADS  PubMed  CAS  Google Scholar 

63.
Walsh, D. B. et al. Drosophila suzukii (Diptera: Drosophilidae): Invasive pest of ripening soft fruit expanding its geographic range and damage potential. J. Integr. Pest Manag. https://doi.org/10.1603/IPM10010 (2011).
Article  Google Scholar 

64.
Delbac, L. et al. Drosophila suzukii est-elle une menace pour la vigne?. Phytoma 679, 16–21 (2014).
Google Scholar 

65.
Poyet, M. et al. Invasive host for invasive pest: When the Asiatic cherry fly (Drosophila suzukii) meets the American black cherry (Prunus serotine) in Europe. Agric. For. Entomol. 16(3), 251–259. https://doi.org/10.1111/afe.12052 (2014).
Article  Google Scholar 

66.
Poulin, B., Lefebvre, G. & Paz, L. Red flag for green spray: Adverse trophic effects of Bti on breeding birds. J. Appl. Ecol. 47, 884–889. https://doi.org/10.1111/j.1365-2664.2010.01821.x (2010).
Article  Google Scholar 

67.
Zeigler, D.R. Bacillus genetic stock center catalog of strains, 7th edition. Part 2: Bacillus thuringiensis and Bacillus cereus. http://www.bgsc.org/_catalogs/Catpart2.pdf (1999).

68.
Gonzales, J. M. Jr., Brown, B. J. & Carlton, B. C. Transfer of Bacillus thuringiensis plasmids coding for δ-endotoxin among strains of B. thuringiensis and B. cereus. Proc. Natl Acad. Sci. USA 79, 6951–6955. https://doi.org/10.1073/pnas.79.22.6951 (1982).
ADS  Article  Google Scholar 

69.
Santos, M., Borash, D. J., Joshi, A., Bounlutay, N. & Mueller, L. D. Density-dependent natural selection in Drosophila: Evolution of growth rate and body size. Evolution 51(2), 420–432. https://doi.org/10.2307/2411114 (1997).
Article  PubMed  Google Scholar 

70.
Bradberry, S. M., Proudfoot, A. T. & Vale, J. A. Glyphosate poisoning. Toxicol. Rev. 23(3), 159–167. https://doi.org/10.2165/00139709-200423030-00003 (2004).
Article  PubMed  CAS  Google Scholar 

71.
R Development Core Team. R: A language and environment for statistical computing. ISBN 3-900051-07-0 https://www.R-project.org (R Foundation for Statistical Computing, Vienna, 2008).

72.
Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01 (2015).
Google Scholar 

73.
Kosmidis I. brglm: Bias Reduction in Binary-Response Generalized Linear Models. R package version 0.6.1, https://www.ucl.ac.uk/~ucakiko/software.html, (2017).

74.
Horton, T., Bretz, F. & Westfall, P. Simultaneous inference in general parametric models. Biometrical J. 50(3), 346–363. https://doi.org/10.1002/bimj.200810425 (2008).
MathSciNet  Article  Google Scholar 

75.
Therneau, T.M., Grambsch, P.M. Modeling Survival Data: Extending The Cox Model. ISBN 0-387-98784-3 (Springer, New York, 2000).

76.
Therneau, T.M. coxme: Mixed Effects Cox Models. R package version 2.2-5. https://CRAN.R-project.org/package=coxme (2015). More