in

Evolution of communication signals and information during species radiation

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 <20 strikes/s.

In order to make an initial assessment of drum types’ discriminability (which, if strong, can suggest a potential to encode species-specific information), we visualized the acoustic space occupied by the drumming of different species by plotting their spatial distribution in the 3d acoustic subspace spanned by the first 3 PCs; Supplementary Fig. 1a). We proceeded similarly using the linear discriminants (LDs) resulting from the DFA carried out to calculate species-specific drumming information content (Supplementary Table 12; Supplementary Fig. 1b; see next section for details on information calculation). This approach was conducted in addition to using the PCs to further validate our assessment of species information encoding in woodpeckers’ drumming.

Life-history data

We used Gorman’s specific description44 to distinguish between species that produce drumming behaviour and those that do not (Fig. 1b). We attributed a ‘drummer’, a ‘non-drummer’ or an ‘occasional drummer’ status when this was clearly stated, and an ‘unknown’ status when the case seemed ambiguous (e.g. when conditional tense was used or when no clear report could be documented). We could thus define a drumming status for each of the 209 species used in reconstructing the ancestral state of this trait.

Similarly, distribution areas in square kilometres were obtained using the same source44 in combination with the website ‘https://www.daftlogic.com/projects-google-maps-area-calculator-tool.htm’ for each of the 92 species included in our acoustic and phylogenetic analyses. Based on these distribution areas, a species pair was defined as sympatric as soon as an overlap was found between the pair’s distributions areas, even if this was only at the range edge. We are aware that higher encounter rates (and thus potentially larger overlapping areas) are more likely to trigger a significant selection pressure for signal divergence between two species. Yet, this approach allows us to be conservative in the criteria used to defined sympatry (e.g. the difficulty of estimating an overlap percentage is much more likely to induce biases), and was supported by matching sympatry levels between our definition and the composition of the communities used in this study. Therefore, the ‘sympatry level’ used in our PGLS regressions corresponds to the number of species sympatric to a given species (e.g. Veniliornins callonotus’ (V.cal) sympatry level is ‘4’, i.e. 4 other species in our remaining sample (n = 91 other species) have distribution areas overlapping with that of V.cal (see Supplementary Data 2).

To assess whether morphological features determine drumming acoustic structure, we collected anatomical measurements on specimens from the Muséum national d’Histoire naturelle (MNHN, Paris, France) and the Natural History Museum (NHM, Tring, UK). We measured beak length, width and height (standardized, as measured at the most posterior point of the beak opening), the wing chord (from the most prominent point of the wrist joint to the most prominent point of the longest primary feather), and tarsus length on its inferior side. After initial inspection of inter-variable correlation, we retained beak length as the single beak measurement, and wing length as the single body size proxy measure78, both to be used in further phylogenetic analyses (see section on PGLS). Wing chord was measured with an Ecotone ornithological ruler, while beak and tarsus measurements were collected using digital calipers (±0.02 mm accuracy for <10 mm measurements and ±0.03 mm accuracy for >10 mm measurements).

We calculated the ‘beak length/wing length’ ratio as a proxy for mechanical constraints on drumming. Drumming can indeed be physically considered as an ‘oscillating spring’, whose motion can be influenced both by beak length and body size53. To standardize our approach and match it to what was done for extraction of drumming acoustic features, we collected measurements from 3 specimens per species and computed median specific value for later analyses. Note that only one specimen was available for Celeus spectabilis (the type specimen) and Picumnus nebulosus.

Body size and body mass data were collected using literature data44,79 and the following websites/archives for species with missing weight data: National Geographic, May 2015 (Dendrocopos noguchii); http://portal.vertnet.org/search?q=Campephilus+pollens (Campephilus pollens). Since we already had a measured proxy for body size with wing chord, we retained only body mass for later PGLS analyses.

Calculation of information

We quantified the species-specific information encoded within each specific drum. To this end, we used two classification algorithms to evaluate the actual discriminative power of drumming signals (and not only that of the classification algorithms), namely a Random Forest classification (RFC) and a DFA. Although the Random Forest algorithm can capture arbitrary groupings of acoustic features, it behaved with lower efficiency in cross-validation (Supplementary Fig. 14). We therefore chose the DFA to generate the confusion matrix of the posterior probability of each drum in our sample (n = 736) as belonging to one of the 92 species being studied. These posterior probabilities were generated in a leave-one-out cross-validation procedure (i.e. 736 different linear discriminant classifiers were trained based on 735 calls to classify the one call that was left out). To prevent overfitting, the DFA was based on the 6 PCA-scores and to ensure equal weighing of each species the DFA was trained with a uniform prior. This resulted in a confusion matrix (Fig. 2b) with 16.5% of correct classification (average of the diagonal), which was significantly higher than expected by chance with 1.09% (pDFA: p < 0.001; Supplementary Fig. 3). From this matrix, we calculated the local mutual information value (measured in bits) based on Shannon’s Entropy42 and following the Eq. (1):

$$MI_L(X_A) = mathop {sum}nolimits_M {left[{pleft( {X_Mleft| {X_A} right.} right)} ast log _2frac{{pleft( {X_M|X_A} right)}}{{pleft( {X_M} right)}}right]}$$

(1)

where MIL is the local mutual information for a given species (XA); (pleft( {X_Mleft| {X_A} right.} right)) is the conditional probability of classifying a drum as belonging to the species XM (M for the model, here the DFA) given that the actual species is XA; and p(XM) is the unconditional probability distribution of predicted species. MIL quantifies the discriminability of one particular species, XA, by taking into account not only the probability of correct classification but also the distribution of classifications both correct and incorrect for that species, (pleft( {X_M|X_A} right),) in comparison to that obtained for the entire dataset, p(XM). The presence of systematic errors can provide additional information that is taken into account in information theory.

The overall mutual information, MI, is then given by the average MIL over species, following Eq. (2):

$$MI = mathop {sum}nolimits_A {[pleft( {X_A} right) ast MI_Lleft( {X_A} right)]} = frac{1}{{n_S}}mathop {sum}nolimits_A {MI_L}$$

(2)

where nS is the number of species. MI can also be related to an average probability of correct detection, pc assuming equal probability of misclassification across all species (X_M ,ne, X_A) (i.e. for scenarios where there is not systematic errors) by inverting the relationship shown in Eq. (3):

$$MI = log _2n_s – left( {1 – p_c} right) ast log _2left( {n_s – 1} right) + p_c ast log _2p_c + left( {1 – p_c} right) ast log _2(1 – p_c)$$

(3)

We also used Eq. (3) to calculate information through evolutionary time (evolution-through-time plots) for hypothetical scenarios (see below ‘Analytical simulations of selection for information’ and Fig. 3c, d) assuming different time courses for pc.

MIL and MI can be normalized by their ‘ceiling value’, namely the maximum amount of species-specific mutual information potentially encoded while discriminating nS species: (MI_{{mathrm{ceil}}} = log _2n_s). Here, (MI_{{mathrm{ceil}}} = log _2(92) = 6.52,{mathrm{bits}}). Ceiling information is reached when the percent of correct classification is 100% for all species. Normalized MI (both overall and local) values range between 0 and 100%. Comparing the normalized MIL values across drumming types showed significant differences (Supplementary Fig. 2b), but these could be the result of the unequal number of species within each drumming type. To control for unequal sample size for that analysis, we also calculated the MIL values selecting randomly 5 species per drumming type (corresponding to the maximum number of species available for the RS and IS drumming types) and iterated this computation 100 times. Comparison of the mean (over n = 100 iterations) normalized MIL across drumming types (each comprising n = 5 species) showed similar results to those found using the full number of species available, with IS significantly encoding more species-specific information than the other drumming types, followed by RS-AC, and then by DK, SS and SF (Fig. 2c).

We also estimated a CI, also obtained from the DFA output and defined by Eq. (4). CI ranges from −100 (minimum classification: for A and B a given pair of species, A is never correctly classified into A and always misclassified into B, and B is never correctly classified into B and always misclassified into A) to 100 (maximum classification: A is always correctly classified into A and never misclassified into B, and B is always correctly classified into B and never misclassified into A).

$${mathrm{CI}}left( {{mathrm{A}},{mathrm{B}}} right) = frac{{left[ {pleft( {X_{mathrm{A}}left| {X_{mathrm{A}}} right.} right) – pleft( {X_{mathrm{B}}left| {X_{mathrm{A}}} right.} right)} right] + left[ {pleft( {X_{mathrm{B}}left| {X_{mathrm{B}}} right.} right) – pleft( {X_{mathrm{A}}left| {X_{mathrm{B}}} right.} right)} right]}}{2} times 100$$

(4)

CI is preferred over normalized MIL for investigating the effect of sympatry on signal information because it allows considering species pairwise discrimination that can be directly compared to pairwise acoustic or phylogenetic distances.

Evolutionary analyses

Phylogenetic generalized least squares

As morphological and ecological factors can play a direct or indirect role in the evolutionary changes in the acoustic structure, we used phylogenetic generalized least squares (PGLS) to examine the current relationships between life-history variables and drumming’s acoustic structure and amount of information. PGLS allow the quantification of these relationships after accounting for effects that could simply be the result of phylogenetic closeness80. Since Miles et al.54 found body size to influence drumming speed and sexual dimorphism to influence drumming length, we gave particular attention to the effect of physical traits (wing length, the ratio of beak length to wing length, and body mass) on drumming structure as well as the effect of geographical distribution traits (sympatry level, size of distribution area) on drumming information (see ‘Methods, Life-history data’ section above for details on life-history variables). We used PC1-PC6 (the components of the PCA carried out on drumming acoustic parameters) as proxies for drumming acoustic structure and the normalized mutual information to quantify the information content about species identity.

PGLS regressions were fitted using restricted maximum likelihood (REML). Model comparison was based on inspection of the Akaike Information Criterion corrected for sample size (AICc), using the null model’s AICc as reference and stepwise forward selection. The improvement of a model was deemed significant only for a decrease in AICc >2 (from the AICc of the null model to the AICc of the fitted model). Model summaries can be found in Supplementary Tables 3, 5 and 6). For variable standardization, prior to running any PGLS model, all life-history variables were z-scored.

For models showing an improvement compared to the null model (i.e. ∆AICc >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.


Source: Ecology - nature.com

Antarctic sea ice may not cap carbon emissions as much as previously thought

A champion of renewable energy