Acoustic communication is crucial for survival across the animal kingdom, with acoustic signals being shaped by the interaction of producer and receiver selective pressures. While spectral features’ variation reflects species-specific selection, the evolutionary history of acoustic communication rhythms, i.e., the rhythmic modulations of acoustic signals, remains unknown. Using data from 98 species spanning primarily mammals and birds, with additional representation from amphibians, reptiles, fishes, and insects, we investigate the origins of acoustic communication rhythms, notably whether they are shaped by the producer’s anatomical characteristics, environmental constraints, or social complexity. Regression models which controlled for phylogenetic relatedness did not support an influence of these species-specific selective forces; instead, explicit phylogenetic models of trait evolution showed that most species’ rhythms are conserved around an evolutionary optimum of 2.7 Hz that falls within the neural delta range (1–4 Hz) and predates mammalian divergence. Given the known conserved brain oscillations across species and delta involvement in active sensing, we propose that, unlike spectral features, acoustic rhythm could be governed by a universal neural mechanism facilitating effective intra and interspecific communication via a shared channel that has persisted through evolutionary times.

Rhythm computation

To explore the evolution of communication rhythm in animals, we analyzed acoustic sequences from 98 species (58 birds, 28 mammals, 4 amphibians, 4 insects, 1 reptile, and 1 fish). We calculated rhythm by analyzing the low-pass variations in the signal amplitude, allowing broad applicability across species (Fig 1a–1d). Validation against conventional methods showed consistent results (Fig 1e: F_2.171 = 0.33, p = 0.72), confirming the robustness of our approach. Control of the impact of signal-to-noise ratio (SNR) and sequence length revealed no significant relationship with rhythm (S1a and S1b Fig). Additional investigation of the allometric relationship between animal weight and the dominant frequency of acoustic signals (S2 Fig), and of the effect of context on rhythm (S4a Fig) further confirmed the validity of our acoustic database.

Species-specific selective pressure

To investigate which factors influence the evolution of rhythm, we fitted two phylogenetic regressions in a Bayesian multilevel framework. A full model investigating the impact of weight, mastication status, and living environment while controlling for phylogenetic relatedness, and a null model controlling for phylogenetic relatedness only. Because insects, fish and anurans are underrepresented in our dataset, we restricted the comparative evolutionary analysis to birds and mammals to ensure robust phylogenetic inference. We extended the analyses to the full dataset in the supplementary materials as an exploratory approach, providing an indication of potential generalization. We compared models via their leave-one-out expected log pointwise density (ELPD) and stacking weight. Due to heteroskedasticity, distributional (scale-location) models better fit the data, leveraging over 90% of the stacking weights (S3 Fig). Including the predictors did not improve predictive performance (Fig 2a), and the null model leveraged the highest stacking weight (Fig 2b). The posterior distribution of the regression coefficients of the full model revealed that none of our predictors had any decisive effect on rhythm, given that their 95% credible intervals (CI) all contained zero. Taken together, these results indicate that including these predictors does not add decisive explanatory value to phylogenetic history. The supplementary analyses including all species yielded comparable results, suggesting that this pattern may generalize to more distantly related taxa (S4 Fig). As vocal complexity, our proxy for social complexity in vocalizing species, and beak morphology in birds were not available for all species, these parameters could not be included in the main phylogenetic model. We therefore conducted separate linear mixed model analyses, regressing rhythm on each factor, with phylogenetic class and order as random effects to account for evolutionary influence. These analyses revealed no significant relationship between rhythm and either vocal complexity (S5d Fig: t = −0.75, p = 0.46) or beak morphology (S5e–S5g Fig: t = −1.4, 0.98, 0.35; p = 0.17, 0.33, 0.73).

Phylogenetic history of rhythm

Visual inspection of the median rhythm for our tested species shows that in every class, acoustic rhythm spans mainly the lower rates <5 Hz (Fig 2d). To understand if the produced rhythms have randomly evolved within this range, favoring species-specific rhythm, or have been maintained around an optimum value, we fitted models of the evolution of vocal rhythm under Brownian Motion (BM) and OU processes [10], representing each evolutionary scenario, respectively. Comparisons of ELPD values and model stackings show that the OU model best fits the data (Fig 3a). The median posterior estimate of rhythm is 2.7 Hz, with 95% CI of [0.45, 4.99] Hz (Fig 3b). This represents both the optimum to which the OU process reverts over time and the likely state at the root of the phylogeny. The posterior distribution of sigma represents the stochastic volatility (Fig 3c) and has a median of 0.68 and a 95% CI of [0.35,2.37]. The posterior distribution of alpha represents the strength of attraction (Fig 3d) and has a median of 1.95 and a 95% CI of [0.42,6.01]. The proportion of posterior half-life estimates lower than 1 (height of the tree) is 95% (Fig 3e), supporting strong selection with fast reversals to the optimum. Consistent with this, rhythm values close to the optimum are also reconstructed for most interior nodes of the phylogeny (Fig 3f). In summary, the model suggests that there is a phylogeny-wide evolutionary pressure towards an optimal rhythm, to which species that would deviate quickly revert to. This result is further supported by models fitted separately for mammals and birds, where the OU model consistently received the highest stacking weight (birds: wOU = 1, mammals: wOU = 0.925) and showed a similar pattern of slow rhythm conservation (S6 Fig). In addition, a supplementary model that included anuran, insect, and fish species yielded comparable outcomes (S7 Fig), suggesting that the conservation of slow acoustic rhythm may have an even older evolutionary origin.

Discussion

In this exploratory study, we report an optimum rhythm for acoustic communication in birds and mammals with credible generalization to additional clades (insects, anurans, fish), challenging the idea that it would be shaped by biomechanical constraints. A biomechanical influence would predict distinct rhythms between masticating and non-masticating species, along with a strong negative allometric relationship between rhythm and weight in masticating species. Instead, our analyses show that rhythm evolution cannot be explained by weight, environmental pressures, or social complexity (Figs 2c and S5d). This surprising absence of effect indicates that biomechanical constraints exert only a limited influence on acoustic communication rhythms at the macroevolutionary scale. By contrast, control analyses on the signals’ dominant frequency (S2 Fig) reconfirmed that spectral characteristics primarily depend on the specifics of the production organ, which chiefly vary with the animal’s weight [11]. Given the minimal contribution of physical characteristics, the lack of explanatory power from environmental and social variables, and the diversity of auditory structures and production mechanisms across taxa, it is unlikely that this shared rhythm results from similarities in anatomy alone.

A key insight comes from comparing potential evolutionary scenarios for rhythm. The phylogenetic analyses clearly favored an OU model over a BM model, indicating strong stabilizing selection rather than unconstrained divergence (Fig 3). This result strengthens the view that species-specific factors, although capable of introducing local variation, do not account for the large-scale conservation of rhythm across distant taxa. With regard to its evolutionary origins, the presence of a similar rhythm across birds and mammals suggests that it was already present in their last common ancestor roughly 340 million years ago (Fig 3). Although insects, amphibians, and fishes are represented by fewer species in our dataset, the presence of slow rhythms in these taxa may hint at an even earlier emergence (S7 Fig). While other factors such as temperature or reverberation are also known to influence communication rhythms in ectotherms, such short-term and context-specific effects are unlikely to generate the macroevolutionary pattern observed here. These findings prompt the consideration of alternative explanations, potentially rooted in neural rather than anatomical, environmental or social factors. As neural mechanisms are essential for survival, they tend to be highly conserved across species and may offer a more plausible basis for a widespread optimal rhythm in communication. Although the size of neurons varies across species, their conductance is kept relatively constant through compensatory variations in ion channel density [22]. This results in conserved time constants for neural responses, including oscillatory phenomena [23].

The conserved rhythm found here around 2.8 Hz (95% interval 0.45–4.99 Hz) best matches delta brain oscillations (1–4 Hz). Interestingly, delta-range oscillations have been described across a wide range of taxa, including mammals, birds, reptiles and insects. In mammals, conserved oscillatory patterns have been reported from humans, primates, carnivores and rodents, with delta among the most stable across species [23]. In reptiles, such as turtles, delta activity dominates brain rhythms and supports ancient functions such as selective attention [24]. Even in birds and insects, delta oscillations have been observed during sleep or quiescent states, indicating that slow oscillatory activity is widespread throughout the animal realm [25]. Delta oscillations support slow perceptual integration and are engaged during active sensing behaviors such as sniffing or blinking [26], suggesting that they provide a fundamental temporal framework through which organisms sample and organize incoming sensory information. Because acoustic signals contain information distributed across both slow and fast temporal scales, as shown in our control analyses on context (S4 Fig), a conserved slow rhythm implies the need for a complementary faster mechanism capable of resolving finer temporal details. Within this neural framework, slow and fast timescales likely serve complementary functions. Delta oscillations would provide a long integration window that supports the identification of acoustic structure, while faster processes (most likely low-gamma) would enable rapid detection and fine temporal discrimination. This division of labor echoes the Fourier uncertainty principle, which states that long temporal windows allow precise frequency resolution at the expense of temporal precision, whereas short windows offer the opposite trade-off [27]. In this view, delta rhythms provide the temporal window best suited for integrating the broad acoustic patterns that carry communicative meaning, while faster oscillations handle local, transient features. Because the present study relies on large-scale comparative acoustic data, it naturally captures the slow integrative timescale rather than the fast detection range. Understanding how these slow and fast processes interact to shape acoustic communication will require dedicated neurophysiological work in future studies.

The independence and parallel functioning of slow and fast scales acoustic processing is attested by human psychophysics findings. When time-compressed unintelligible speech is chunked and repackaged by interspersing silent gaps to restore regular syllable and word rates, speech becomes intelligible again. These results show that speech comprehension is limited by the packaging of information at slow rhythm, rather than by the absolute capacity to decode accelerated segmental (phonemic) cues [6,28]. This emphasizes the relevance of slow rhythms in maintaining acoustic signal decodability. Human speech operates at a mean rate that lies at the higher end of the median universal rhythm identified in this study, clustering around 5 Hz across all languages [8]. This reflects humans’ tendency to convey more information per unit of time compared to other species. Dogs, for instance, despite millennia of domestication and a demonstrated ability to understand many aspects of human language [29], process speech at a slower rhythm, typically between 1 and 3 Hz, than the human average of 5 Hz [30]. While they do not track the syllabic rhythm of speech, dogs tend to process it more globally, at the word rate, which aligns with their own vocal production rhythm (e.g., bark rate). These comparative findings highlight species-specific variations in preferred communication rhythms, while also underscoring a broader pattern: a conserved and effective temporal resolution for communication that enables interspecies sensitivity to acoustic signals. In sum, the conserved production rhythm identified here aligns closely with an equally conserved neural rhythm that facilitates the assignment of meaning to acoustic input, ultimately supporting an effective and evolutionarily robust communicative system.

In conclusion, the maintenance throughout evolution of a common slow rhythm for acoustic signal production and perception de facto results in a common communication channel across coexisting species, offering possibilities for interspecies signaling (e.g., a common danger) and/or eavesdropping, and thus conferring evolutionary advantages. Spectral (frequency) and temporal (rhythm) features, both key to effective communication, seem to have followed different trajectories. While spectral traits diversify with the production and hearing organ anatomy, acoustic rhythms might be shaped by basic neural factors essential for survival that have led to a slow conserved optimum. Whether these neural factors were maintained because they amounted to an optimal communication design, or whether this design arose from recycling vital and highly conserved neural traits (neural rhythmicity via maintained neuronal conductance) remains to be elucidated. A further interesting avenue will be to probe for an optimum that combines fast and slow timescales in acoustic communication, alongside testing the delta wave hypothesis by assessing whether interspecific variation in rhythm is accompanied by corresponding variation in neural timescales.

Vocal sequences

To perform an extensive phylogenetic comparison on a balanced phylogeny, we collected acoustic and biological data for at least one species per infra-order of tetrapods, when data were available, as well as a few species of insects and fishes, to obtain a good representation of rhythm throughout the phylogeny. Acoustic sequences were gathered from public and private databases (Fonozoo, Cisro, Berlin Museum fur Naturkunde and Xeno-canto), online videos platform (Youtube, Dailymotion), and from different research groups that kindly shared audio files.

Following a cross-species literature survey [31–34], we defined a sequence as a recording of an acoustic display produced by a single individual, containing more than two calls separated by less than two seconds of silence. For computational purposes, we only retained sequences lasting more than one second.

Data collection then followed a strict and standardized protocol. For each order, we first screened public databases and ranked species by the number of available sequences. We then selected the first species fulfilling the following criteria: at least five clean sequences (non-noisy, with a single individual vocalizing at a time) originating from different individuals. Once such a species was identified, we proceeded to the next infraorder. When no species meeting these criteria was available, we extended the search to private databases and directly contacted researchers. If no data could be obtained, that infraorder was excluded.

Only recordings with minimal background noise and a clearly identifiable focal individual were retained; any containing rhythmic background noise, overlapping individuals, or ambiguous sound sources were excluded. This careful filtering minimized the risk of confounding effects and ensured that the extracted acoustic envelopes reliably reflected the focal signal across different sources.

The final dataset comprised 98 species: 58 birds, 28 mammals, 4 amphibians, 4 insects, 1 reptile, and 1 fish. Species with fewer than five sequences were included only when they represented the sole available species for their infraorder.

Weight

Weight presents an allometric relationship with morphological features and physiological processes involved in acoustic productions across various species. Thus, we collected biological data on the mean weight for each species, to serve as a proxy for heart rate, breathing rate, and metabolism. This involved calculating the average weight by considering both the minimum and maximum weights recorded for each species, irrespective of sex. These data were primarily obtained from the handbook of mammals of the world [35], and the handbook of birds of the world [36]. When data were unavailable from these sources, we looked for reference articles.

Beak size

Just as masticatory abilities may have influenced rhythm in mammals, a similar proposition could be made regarding beak morphology in birds. Unlike other animals, the morphological traits of a bird’s beak do not consistently adhere to an allometric relationship with its weight [37]. Thus we collected information on beak length, width and depth of our species. As these measures were not available for all species, and only applicable to birds, we could not include them in the phylogenetic model. We therefore built an additional linear mixed model investigating the variation of rhythm, including phylogenetic class and order as random effects to account for evolutionary influence, and beak length, width, depth, and their interactions as fixed effects. Data are available in the study github.

Living environment

As environmental conditions can impact vocal communication [38], we also collected data on the typical habitat of each species. We used a five level categorical classification, with habitats being either classified as closed (defined as habitats with heavy tree coverage), semi-closed (defined as habitats with light three coverage or human cities), open (defined as fully open habitat with no three coverage or obstacle), shore (for species living near a significant amount of water such as lakes, rivers or seas) or water (for species living below the water surface). These data were primarily obtained from the handbook of mammals of the world, and the handbook of birds of the world. When data were unavailable from these sources, we looked for reference articles. Data and linked references are available in the study github.

Mastication status

As some have proposed that rhythmic communication in vocalizing animals may be linked to mastication regime [14–16], we also classified each species according to their mastication status (yes or no).

Social complexity

As social complexity increases, individuals may need to communicate more information in a given time, and therefore speed up their communication. As communication signals have been linked to social complexity [20], we also gathered species vocal repertoire complexity (number of distinct calls in the species vocal repertoire) when this information was available. As this measure was not available for all species, we could not include it in the phylogenetic model. We therefore built an additional linear mixed model investigating the variation of rhythm, including phylogenetic class and order as random effects to account for evolutionary influence, and vocal complexity as fixed effect. These data were primarily obtained from the handbook of birds of the world, and reference articles. Data and linked references are available in the study github.

Rhythm analyses

To quantify rhythm in these acoustic sequences, we decided to adapt the method developed by Tilsen and colleagues to compute rhythm in human speech production [39]. This method uses the signal amplitude to automatically compute the rhythmic component of a sequence, without making any assumption on the components’ size, and is thus widely applicable across all species regardless of variations in unit size or spectral characteristics. First, we denoised the sequences using a first-order Butterworth filter, with a bandpass filter between a minimal frequency (minF), defined as 200 Hz below the minimum frequency of the animal call, and a maximum frequency (maxF), defined as 200 Hz above the maximum frequency of the animal acoustic signal obtained from reference articles. When this information was not available, we applied a large range filter with a 100 Hz minF and 10,000 Hz maxF. We then computed the normalized envelope of the denoised sequences using the Hilbert Transform. Next, we low-pass-filtered this envelope with a fourth-order Butterworth filter with a 20 Hz cut-off frequency to obtain the slow changes in acoustic energy. Before further analysis, we downsampled the resulting signal at 150 Hz for computational purposes, and then applied a continuous wavelet transform using the Morlet wavelet to obtain a time-frequency representation of the amplitude envelope. We replaced the Fourier transform with a wavelet transform, to allow for more flexibility with regards to the variation in sequence length present in our dataset. We finally analyzed that representation’s power spectrum to extract the five frequency peaks of highest amplitude in the power spectrum and used the time-frequency representation to select the main rhythmic component conserved across the entire sequence.

To assess the validity of the proposed methodology, we conducted a comparative analysis between the calculated vocal rates of a subset comprising 10% of our database and those derived from two widely accepted conventional approaches: (1) by counting the number of elements per second (Count) and (2) by computing the inter-element interval (IEI). The three methods gave sensibly similar results (Fig 1d: F_2.171 = 0.33, p = 0.72), hence validating our rhythm quantification method.

Signal-to-noise ratio (SNR) and length

As further control analyses, we quantified recording durations in seconds and SNR in decibels. To control for the effect of both factors and their interaction on rhythm, we build a linear mixed model investigating the variation of rhythm including group and order as random effects, and SNR, length and their interaction as fixed effects.

Dominant frequency analyses

To determine dominant frequencies, which unlike the fundamental frequency are measurable in all types of communicative signals [40], we isolated the first acoustic unit in each denoised sequence. We then applied a single discrete Fourier transform to compute the power spectrum of these units, and extract the peak of highest amplitude. The obtained results were also visually controlled in Praat, to make sure that the extracted dominant frequency matched the acoustic energy present in the unit. If the first unit had poor SNR leading to inaccurate computation of the dominant frequency, we selected the next unit in the sequence.

Context effect control

While the existing literature highlights the importance of context and its correlated arousal levels on acoustic signal rate [1], for most species we were not able to obtain these data. Nevertheless, whenever possible we selected recordings of different call types for each species. Further, we performed separated analysis of variance (ANOVA) to control for the effect of call type on rhythm in three species: one avian and two mammalian, from which we could obtain different call types, including contact calls, alarm calls, songs, agonistic and antagonistic vocal displays.

Phylogenetic tree sample reconstructed from genetic sequences

To test our hypotheses of interest we used phylogenetic comparative methods (PCM), a broad family of methodological tools for characterizing and controlling for the evolutionary dynamics thought to give rise to the data under study. PCMs require a representation of the relatedness of the taxa under study in the form of a phylogenetic tree sample. To represent the tree topology of the species we performed a phylogenetic analysis based on comparable genetic sequences, using a Bayesian framework to infer a posterior tree sample. We first matched each species in the sample with their closest genetic proxies in GenBank [41], and extracted mitochondrial DNA for the corresponding species. For 54 of the 98 species, matching mitochondrial genomes were available from literature and deposited in GenBank. For the 44 remaining species we chose proxies from another closely related species. We took a species within the same genus when possible; if this was not available, we chose species within the family of the target species, after confirming that no more than one species per family was included in the original list. Only for four target species—Correlophus ciliates (Squamata), Galbula ruficauda (Aves), Leptosomus discolor (Aves), Phaethon rubricauda (Aves) —we did not find proxies within the family and we had to find a proxy within the order. To choose the best mtDNA proxy with those deposited in GenBank, we considered completeness of the available mitochondrial sequences and comparable average size, weight, and environment of the target species. Maximum missing data is 100 base pairs, for an average size of 16706 base pairs. MtDNA genomes were aligned with MAFFT software [42] and standard settings. The alignment was manually screened in BioEdit (version 7.2, https://thalljiscience.github.io/) for spotting irregularities and potential outlier sequences. Sequences were then cut to keep only the coding region, which is more conserved across species, using the Homo sapiens sequence as a reference. The final alignment consisted of 21860 base pairs, which include large INDELs sections to accommodate alignment between the most divergent species (e.g., Apis mellifera).

We used BEAST2 to generate the trees, running 10′000′000 iterations of Markov chain Monte Carlo (MCMC) with a thinning interval of 1,000. We used the following settings to approximate the broad evolutionary range of the species considered: assuming an HKY substitution model, a strict clock (Uniform rates across branches), and a Birth-Death tree prior with a Yule birth rate. This resulted in 10′000 trees, of which we use 50 for phylogenetic comparative analyses.

Bayesian multilevel models for phylogenetic regression

To assess the impact of several predictors of interest on different properties of variation in vocal rhythm and dominant frequency, we used phylogenetic regression modeling, a comparative method that assesses the effect of predictors on a response while controlling for the phylogenetic relatedness of the taxa. Due to heterogeneity in the number of datapoints and individuals in each species, we employed both non-distributional regression models, which model the mean of the response variable as a function of predictors, and distributional (scale-location) regression models, which model both the mean and standard deviation of the response variable as a function of predictors. We control for species-level idiosyncrasies in both the median and (in some cases) standard deviation of rhythm via phylogenetic random intercepts and slopes (phylogenetic random effects are similar to the standard random effects used in hierarchical regression modeling, but are generated by a Gaussian Process with a covariance kernel that is a function of the phylogenetic patristic distances between species under study rather than independently and identically distributed with diagonal variance).

We fitted four phylogenetic regression models using brms [21] for each response variable (vocal rhythm, dominant frequency), resulting in eight models. Due to the constraints of bayesian multilevel modeling, particularly its sensitivity to missing data across hierarchical levels, we restricted the main model to predictors available for the majority of species. Specifically, we included average species weight (log-transformed, centered, and standardized) [43], mastication status, and species’ living environment. We also modeled the interaction between mastication status and weight. Additional predictors that were only available for a limited subset of species were analyzed separately in independent linear mixed models to avoid excessive data imputation or loss of statistical power.

Our distributional models have the following basic generative process (below, is shorthand for all model predictors, including fixed and random effects):

Non-distributional models have the following structure:

We employ the default priors of brms.

The first of the four models was a full distributional one that modeled both the expected median and variance of the response variable as a function of these predictors, while controlling for species-level idiosyncrasies in both the median and variance of rhythm via phylogenetic random intercepts and slopes. The second was a full non-distributional model that treated only the expected median rhythm as a function of the predictor variables as well as phylogenetic random intercepts and slopes. The third of these was a null distributional model that included only phylogenetic random intercepts and slopes for the expected median and variance. The final model was a null non-distributional that included only phylogenetic random intercepts and slopes for mean rhythm. We ran each of these models for 4,000 iterations of the no U-turn sampler over 4 chains with a log-normal link function and discarded the first half of samples, aggregating posterior samples across the retained sampled trees.

Models in brms have the following formulae:

Full, distributional

Full, non-distributional

Null, distributional

Null non-distributional

We compared fitted models via their leave-one-out ELPD values [44] and stacking [45], which average predictive distributions of different models to generate weights representing their relative predictive power. We used the function loo_compare to measure differences in ELPD across models. Finally, we inspected posterior distributions of regression coefficients of the full distributional model to assess the effects of predictors of interest.

Evolutionary dynamics of vocal rhythm

We further investigate the properties of rhythm across species using two Gaussian Process models of continuous trait evolution, asking specifically whether the evolution of vocal rhythm is characterized by a random process of drift (characterized by BM) or whether selective forces draw rhythm values toward an optimal value over time (a mean-reverting scenario characterized by an OU process). Under BM, the displacement of a continuous trait at time s has a variance proportional to the amount of time elapsed over the course of displacement (below denoted as t), where represents the scale of the drift process:

Under an OU process, the displacement of a character has the following formula:

In the first component of the sum, represents the strength of selection to the optimal value . The second component represents a process of BM, with representing the scale of drift. Thus, the OU process allows for both selective and random forces in character evolution.

An standard way to interpret is to transform it to the phylogenetic half-life, ln 2/α [45]. This is interpreted as the average time for a trait to evolve halfway from an ancestral state toward a new optimum, indicating how long it will take before adaptation to a new regime is more influential than constraints from the ancestral state. If half-life values are greater than the height of the phylogeny (1 in our case, as the tree length is scaled to unit height), the process increasingly resembles BM and involves a slower adaptation speed.

As above, we employ a distributional approach, allowing species-level mean rhythm values and species-level standard deviations of rhythm values to evolve over the phylogeny according to BM or OU processes.

The distributional BM process has the following generative process:

represents the trait value at the root of the tree, while is a matrix of the shared history (the time between the root age of the tree and the most recent common ancestor) of each pair of nodes in the tree and is the positive scale of drift. Conventions are as above.

The distributional OU process has the following generative process:

represents the trait value at the root of the tree, is the positive scale of drift, α is the positive strength of selection, and Δ is a matrix of pairwise cophenetic distances between species in the phylogeny, scaled to a maximum distance of 1. Conventions are as above. We place Normal(0,1) priors over unconstrained parameters and Gamma(1,1) priors over positive parameters. In addition to running these models on all species in our sample, we validate results by running models on bird and mammal species alone.

Phylogenetic reconstruction of rhythm values

Rhythm values were reconstructed to internal nodes of the maximum clade credibility (MCC) tree of the phylogeny, using ggtree package [46] by drawing 10 draws from each of the posterior distributions inferred from the 50 different trees in the tree sample and sampling values at internal nodes of the tree from the normal distribution parameterized by the OU process, conditioned at the expected tip values (Fig 3f).

Software

All analyses and visualization were done using Stan and R version 4.1.2 (2021-11-01) with the following packages Seewave [43], Soundgen [47], DoBy [48], Lme4 [49], MuMYn [50], brms [21], ggplot2 [51].

S1 Fig. Signal-to-noise ratio and sequence length effect on rhythm.

a) Rhythm (Hz) as a function of signal-to-noise ratio (SNR) showing the absence of relationship between SNR and Rhythm (t = −0.21, p = 0.84, R² = 0.02). b) Rhythm (Hz) as a function of recording length showing the absence of relationship between sequence length and rhythm (t = −1.45, p = 0.15, R² = 0.02). The data underlying this Figure can be found in https://zenodo.org/records/19816728.

https://doi.org/10.1371/journal.pbio.3003798.s002

(TIF)

S2 Fig. Dominant frequency phylogenetic regression.

a) Leave-one-out expected log pointwise density difference (ELPD) between the full and the null distributional models. b) Stacking weight of the models. c) Posterior credible interval (95% and 85%) of the full model. d) Dominant frequency plotted on a logarithmic scale as a function of log-transformed weight with predicted slopes from the full distributional model and their standard error. The data underlying this Figure can be found in https://zenodo.org/records/19816728.

https://doi.org/10.1371/journal.pbio.3003798.s003

(TIF)

S4 Fig. Bayesian multilevel model including all species.

a) Leave-one-out expected log pointwise density difference (ELPD) between the null and the full distributional (scale-location) models. b) Stacking weight of the models. c) Posterior credible intervals (95% and 85%) of the full distributional model (W*M = Weight*Mastication). d) Rhythm plotted as a function of log-transformed weight with predicted slopes from the full distributional model and their standard error for masticating and non-masticating species. e) Distribution of raw median rhythm (in Hz) across birds and mammals. Each vertical black line represents the median rhythm of a single species. The data underlying this Figure can be found in https://zenodo.org/records/19816728.

https://doi.org/10.1371/journal.pbio.3003798.s005

(TIF)

S5 Fig. Additional analysis on rhythm.

a) Rhythm in sequences of different call types (1 = affiliative grunt, 2 = scream, 3 = threat grunt) in olive baboons (Papio anubis). b) Rhythm in sequences of different call types (1 = bark, 2 = growl, 3 = howl, 4 = snarl, 5 = whine) in dogs (Canis lupus familiaris). c) Rhythm in sequences of different call types (1 = alarm call, 2 = flight call, 3 = song) in Eurasian stone-curlews (Burhinus oedicnemus). d) Rhythm plotted as a function of vocal repertoire size. e) Rhythm plotted as a function of beak length in birds. f) Rhythm plotted as a function of beak depth in birds. g) Rhythm plotted as a function of beak width in birds. The data underlying this Figure can be found in https://zenodo.org/records/19816728.

https://doi.org/10.1371/journal.pbio.3003798.s006

(TIF)

S6 Fig. Independent phylogenetic modeling for birds (a–e) and mammals (f–j).

a,f) Stacking weight of the OU and BM distributional models. b,g) Posterior distribution of rhythm in the OU distributional model. c,h) Posterior distribution of sigma, the scale of the drift process. d,i) Posterior distribution of alpha, the strength of selection. e) Posterior distribution of the phylogenetic half-life. e,j) Visualization of reconstructed median posterior values using a maximum clade credibility (MCC) tree. The data underlying this Figure can be found in https://zenodo.org/records/19816728.

https://doi.org/10.1371/journal.pbio.3003798.s007

(TIF)

S7 Fig. Phylogenetic modeling including all species.

a) Stacking weight of the OU and BM distributional models. b) Posterior distribution of rhythm in the OU distributional model. c) Posterior distribution of sigma, the scale of the drift process. d) Posterior distribution of alpha, the strength of selection. e) Posterior distribution of the phylogenetic half-life. f) Visualization of reconstructed median posterior values using a maximum clade credibility (MCC) tree. The data underlying this Figure can be found in https://zenodo.org/records/19816728.

https://doi.org/10.1371/journal.pbio.3003798.s008

(TIF)