Phylogeny

Our sampling of species for acoustic communication was guided by which species were included in time-calibrated trees and had diel activity data. Time calibration is essential for combining trees, estimating diversification rates, and reconstructing trait evolution. We started with a tetrapod phylogeny that included 1824 species22. This tree had four main advantages. First, all species already had diel activity data. Second, all major tetrapod clades were included (i.e., amphibians, mammals, birds, lepidosaurs, turtles, and crocodilians), with species sampled within each clade in proportion to the clade’s overall species richness (see details below). Proportional sampling is important for analyses relating traits and diversification. The sampling of species within these major clades was also designed to represent major groups (e.g., orders, families), with sampling roughly proportional to their species richness. The tree was generated by combining trees from studies focused on each major clade (amphibians39, birds31, crocodilians40, lepidosaurs41, mammals42, turtles43), linked by an overall phylogeny among major clades44. This latter tree has an estimated topology and divergence times that are similar to those from an extensive phylogenomic study45 but has better sampling of major clades. There were two versions of the overall tree22, differing in their backbone trees within birds (i.e., Ericson29 and Hackett30 trees31). Both trees were used for all analyses here. Third, although the tree has extensive taxon sampling, it is not so large as to make computationally intensive analyses impractical (e.g., HiSSE). Fourth, because the phylogeny was assembled for analyzing diel activity, the species sampling should be unbiased regarding whether species have or lack acoustic communication. Note that species sampling in the diel activity dataset was not designed to favor any particular diel activity state but rather to represent higher taxa (e.g., classes, orders, families) in proportion to their species richness (see above). The trees are available as Supplementary Data 1 and 2.

We recognize that there are more recent phylogenies available within some groups. However, many of these have more limited taxon sampling. Furthermore, we found the presence of acoustic communication to be largely invariant within most major clades (e.g., birds; Table 1). Therefore, alternative phylogenies within these groups should have limited impact on the overall results.

Acoustic communication data

We obtained acoustic communication data from books, original papers, and online databases (e.g., xeno-canto: https://www.xeno-canto.org/). Data on individual species and supporting references are given in Supplementary Data 3. We searched the literature for data on each species in the tree using Web of Science and Google Scholar between December 1, 2017 and January 15, 2019. The following keywords were used: “acoustic communication,” “call,” “vocal communication,” “vocalization,” “song,” and “sound.” Each keyword was used in combination with the species name for every search.

Acoustic communication can be defined as the transmission of intraspecific messages via airborne sound waves46. We assigned species to one of the two states following this definition (0 for absence of acoustic communication, 1 for presence). Species that produce sound but are not known to exchange information from sound with conspecifics were considered to lack acoustic communication (e.g., rattling in rattlesnakes). Vibrations transmitted through a fluid (air or water) are generally defined as “sound,” whereas those in solids are generally referred to as “vibrations”47. Species utilizing only substrate-borne vibration communication were also characterized as lacking acoustic communication.

We recognize that there may be some cases in which sounds are produced but communication with conspecific individuals is present but not yet documented. We also recognize that conclusively documenting the absence of a behavior is challenging. We address why our conclusions should not be an artifact of these factors in the “Discussion” section.

There were 25 species initially included in the trees for which we could not find adequate data on the presence or absence of acoustic communication. These were pruned from the trees using Mesquite48 version 3.51. We avoided removing two mammal species by replacing species lacking data on acoustic communication data with congeners having these data. Specifically, we replaced Heteromys anomalus with H. desmarestianus and Liomys adspersus with L. pictus.

We obtained data for a total of 1799 tetrapod species (Table 1). We assembled estimates of the total described species richness of amphibians49, non-avian reptiles50, birds51, and mammals52. Based on these estimates, our sampling within clades was strongly proportional to each clade’s total richness (r2 = 0.936; P = 0.0016). Nevertheless, amphibians were somewhat oversampled relative to other clades (Table 1). However, described amphibian richness is still increasing rapidly, with ~700 species added in the past 3 years49 (2016–2019). Thus the actual number of amphibian species might be similar to the numbers of described lepidosaurs and birds.

Testing relationships between acoustic communication and diel activity

We used diel activity data assembled in a previous study22, following their definitions. In brief, species were classified as being diurnal (primarily active by day), nocturnal (primarily active by night), crepuscular (active primarily at twilight and/or early morning), or arrhythmic (cathemeral; similarly active during the day and night or changing day–night activity seasonally). Many analyses required treating diel activity as a two-state character, and most tetrapod species22 seem to be primarily nocturnal (41%) or diurnal (51%). Therefore, they alternatively coded all crepuscular and arrhythmic species as either nocturnal or diurnal. Under “maximum diurnal” coding, crepuscular and arrhythmic species were coded as diurnal. Under “maximum nocturnal” coding, these species were coded as nocturnal. We used these strategies here, given that our analyses also require two-state coding. Diel activity data are given in Supplementary Data 4.

We used two approaches to test whether the presence of acoustic communication is associated with nocturnal activity. First, we compared maximum-likelihood models of independent and dependent evolution between diel activity and acoustic communication23. These models test how one trait affects the transition rates of the other trait. Analyses were performed using the R package phytools53 version 0.6-60. We primarily tested whether acoustic communication depends on diel activity, but we also tested whether diel activity depends on acoustic communication, and whether both depend on each other. Prior to these analyses, we compared the fit of models with equal transition rates (i.e., gains and losses, 0 to 1 and 1 to 0, respectively) between states of each character (ER) and models with all transition rates different between states (ARD). Models were compared with the function fitDiscrete in GEIGER54 version 2.0 using AIC values (lowest AIC considered best-fitting28). For maximum diurnal coding, the best fit model was ARD, but AIC differences between models for maximum nocturnal coding were <4 (Supplementary Table 2). For acoustic communication, AIC differences were also <4 (Supplementary Table 1). Therefore, we used both ARD and ER models for both variables for these tests. Likelihood-ratio tests were used to evaluate whether a model of correlated evolution explained the data better than the null model of independent evolution of the two traits. In addition, models were compared using the AIC and AIC w statistics28. The AIC w measures the weight of evidence in favor of a model, on a scale from 0 to 1. The model with the highest AIC w had the best fit.

For the second approach, we used phylogenetic logistic regression24. This approach can be used to test for relationships between two categorical variables. We utilized the function phyloglm in the R package phylolm55 version 2.6. We used the IG10 phylogenetic generalized linear model24. This method modifies a GLM (generalized linear model) framework for binomial distributions to incorporate a correlation matrix representing relationships among taxa. We treated diel activity as the independent variable and acoustic communication as the dependent variable.

Diversification

To test whether acoustic communication increased diversification rates, we used the BiSSE (Binary State Speciation and Extinction) and HiSSE approaches in the R packages diversitree56 version 0.9–7 and HiSSE57 version 1.8. Importantly, HiSSE incorporates unmeasured factors (i.e., hidden states) that could impact diversification rates besides the included character. We compared the fit of five different models. We first ran two BiSSE models: (i) a full BiSSE model with two observed states (acoustic communication absent or present, 0 and 1, respectively), with different rates of transitions (q 01 , q 10 ), speciation (λ 0 , λ 1 ), and extinction (μ 0 , μ 1 ) associated with each state. (ii) Second, a restricted BiSSE-equivalent model57 with two observed states and no hidden states, constrained speciation and extinction rates (λ 0 = λ 1 , μ 0 = μ 1 ), and two transition rates (HiSSE two-state model). We then ran a full HiSSE model (iii) with two hidden states (A, B) contained within each observed state (i.e., states 0A, 1A, 0B, 1B), speciation and extinction rates varying independently across all four states, and transition rates between all observed and hidden states free to vary except for dual transitions (e.g., q1A to q0B, q0A to q1B). These dual transitions were considered unlikely and were disallowed57. In addition, we tested two null HiSSE models: (iv) HiSSE Null2 (two hidden states), with different speciation and extinction rates only between the hidden states and different transition rates between observed and hidden states, and (v) HiSSE Null4 (four hidden states) with two additional hidden states (C, D), different speciation and extinction rates for each hidden state but not for the observed states, and equal transition rates across all the hidden and observed states. The first model was implemented in diversitree, and models ii–v were implemented in HiSSE. All analyses were conducted with default root setting “madfitz”38. Comparisons were conducted for both trees. Models were evaluated based on AICc scores28.

Note that BiSSE/HiSSE analyses can account for incomplete sampling of species by correcting for biased sampling of states among species38. However, we considered our sampling of states among species to be unbiased, and so no correction was considered necessary.

Testing for phylogenetic signal

Prior to reconstructing the evolution of acoustic communication across the tree, we evaluated whether this trait was phylogenetically conserved. Phylogenetic signal was first evaluated using Pagel’s26 lambda test, utilizing the fitDiscrete function in GEIGER54. Low lambda values (close to 0) indicate weak phylogenetic signal, whereas higher values (close to 1) indicate strong signal26. In addition, we also compared the fit of a model based on the estimated lambda value (EL) to a white-noise (WN) model, with no phylogenetic signal. Likelihoods were compared using the AIC. AIC differences >4 between models were considered strong support for the best-fitting model28.

Phylogenetic signal was also assessed using the D-statistic27. The D-statistic was designed to test for phylogenetic signal in binary traits. We used the phylo.d function in the R package caper58 version 0.5.2. The D-statistic compares the observed D-value to alternative D-values generated with simulated data based on Brownian motion (BM; strong phylogenetic signal) and WN models (no signal). Negative D-values indicate traits are strongly clumped (i.e., more conserved than BM), D = 0 supports the BM model, and D = 1 supports the WN model27. We also tested whether the D-statistic differed significantly from the BM and WN models. Non-significant P values for BM indicate conserved trait evolution. Note that phylogenetic signal will reflect low evolutionary rates (and phylogenetic conservatism) when analyzing discrete variables59, as done here.

Ancestral-state reconstruction

We reconstructed the evolution of acoustic communication across the tree using maximum-likelihood ancestral-state reconstruction. Different diversification rates associated with different states can potentially impact reconstructions60. However, the HiSSE Null4 model had the strongest support for both trees, suggesting that acoustic communication did not significantly impact diversification rates (Table S5). Therefore, we used simpler models that did not incorporate diversification (ARD and ER, see above). Ancestral reconstructions were performed using both the ARD and ER models, given the similar fit between these models (Table S1).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.