Latitudinal variation in mode of phenotypic evolution

We tested whether modes of phenotypic evolution varied with latitude using 2 types of models. First, we tested whether support for various “single-regime” models that estimate a single set of parameters on the entire avian phylogeny [26] varied according to a clade-level index of tropicality. Second, we developed and used “2-regime” models with distinct sets of parameters for tropical and temperate species and tested whether these latitudinal models were better supported than single-regime models.

Across single-regime fits, we found no evidence for a latitudinal trend in the overall support for any model of phenotypic evolution (Fig 1A–1F and S4 Table), with one exception: There was an increase in model support for the MC model in tropical lineages for the locomotion phylogenetic principal component (pPC) 3 (Fig 1F and S4 Table). Similarly, there was no evidence that the overall support for models incorporating competition (i.e., MC or DD models) is higher in tropical clades (Fig 1G and S4 Table). Models with latitude (i.e., 2-regime models) were not consistently better supported than models without latitude, for any model or trait (S5 Table). Indeed, single-regime models were the best-fit models across 86% of individual clade-by-trait fits (S7 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Model support for single-regime models reveal little impact of latitude on the mode of phenotypic evolution in birds (66 clades with ≥50 species, with data from 7,163 species).

There is no relationship between the proportion of taxa in a clade that breed in the tropics and statistical support (measured as the Akaike weight) for (a) BM, (b) OU, (c) EB models, (d) DD_exp models, or (e) DD_lin models. In MC models (f), there is an increase in model support for locomotion pPC3 (solid line). The relative support for a model incorporating competition (i.e., MC or DD models) does not vary latitudinally for any trait (S4 Table). Each point represents the mean Akaike weight across clade-by-trait fits to stochastic maps of biogeography (i.e., each clade contributes a point for each of 7 traits; see S2 and S3 Datas). BM, Brownian motion; DD, diversity-dependent; DD_exp, exponential diversity-dependent; DD_lin, linear diversity-dependent; EB, early burst; MC, matching competition; OU, Ornstein–Uhlenbeck; pPC, phylogenetic principal component.

https://doi.org/10.1371/journal.pbio.3001270.g001

Latitudinal variation in the effect of interactions on phenotypic evolution

We found no evidence for a latitudinal trend in the slope estimated from single-regime DD models (Fig 2C and 2D and S6 Table). However, the strength of repulsion estimated from single-regime MC models increased in more tropical families for locomotion pPC3 (Fig 2B and S6 Table). Parameter estimates from 2-regime models with competition (i.e., MC or DD models) do not support a stronger effect of biotic interactions on phenotypic evolution in the tropics (Fig 3B–3D)—in most traits, there is no consistent difference between estimates of the impact of competition on tropical and temperate lineages, and in one case (bill pPC2), there is evidence that competition impacts temperate lineages to a larger degree than tropical lineages (Fig 3B–3D and S7 Table). In all cases, there was substantial variation in the fits, and the overall magnitude of differences between tropical and temperate regions was rather small (Fig 3B–3D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Parameter estimates from single-regime models reveal varying impacts of latitude.

There is no impact of latitude on the effect of competition on trait evolution as measured by the slope of (a) DD_exp models or (b) DD_lin models. (c) The effect of competition on trait evolution as measured by the repulsion parameter (“S”) from the MC models increases with the index of tropicality (the proportion of species in the clade with exclusively tropical breeding distributions) for locomotion pPC3 but not for other traits. (d) There is no relationship between the proportion of taxa in a clade that breed in the tropics and the estimated rate of trait evolution from BM models. Solid lines represent statistically significant relationships (S6 and S13 Tables). For (a–c), each point represents the mean across clade-by-trait fits to stochastic maps of biogeography (for all families with at least 50 species), and for (d), each point represents the MLE for each clade-by-trait fit (see S2 and S3 Datas). BM, Brownian motion; DD, diversity-dependent; DD_exp, exponential diversity-dependent; DD_lin, linear diversity-dependent; MC, matching competition; MLE, maximum likelihood estimate; pPC, phylogenetic principal component.

https://doi.org/10.1371/journal.pbio.3001270.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Parameter estimates from 2-regime models reveal varying impacts of latitude.

Estimates of slopes from (a) DD_exp models and (b) DD_lin models are not consistently different in tropical regions in any trait. (c) MC models estimated a decreased effect of competition in the tropics on bill pPC2. (d) Estimates of evolutionary rates from BM models show accelerated rates of locomotion pPC3, but not other functional traits, in temperate regions. Asterisks indicate statistical significance (S7 and S14 Tables). For (a–c), each point represents the mean across clade-by-trait fits to stochastic maps of biogeography and of tropical/temperate membership (for all families with at least 50 species), and for (d), each point represents the mean across stochastic maps of tropical/temperate membership maximum (see S4 and S5 Datas). BM, Brownian motion; DD, diversity-dependent; DD_exp, exponential diversity-dependent; DD_lin, linear diversity-dependent; MC, matching competition; pPC, phylogenetic principal component.

https://doi.org/10.1371/journal.pbio.3001270.g003

Impact of assuming continental scale sympatry

Phylogenetic models of competitively driven trait evolution rely on reconstructions of ancestral ranges to delimit the pool of potential species interaction at each point in the evolutionary history of a clade. Given the scale of our analyses and the computational limits of existing models of ancestral range estimation, we assumed that species occurring on the same continent were able to interact with one another. On average, species in our analyses are sympatric with 50% of clade members at the continental level, although there are differences across continents (mean range: 34% to 74%; S5 Fig and S9 and S10 Tables). Notably, we also found that temperate species are more likely to coexist in sympatry with family members than tropical species (S11 Table). To determine the impact of assuming continental scale sympatry, we investigated whether we would detect a latitudinal difference in the effect of competition on phenotypic evolution if it existed, even if competition occurs among only truly sympatric species rather than among all species occurring on the same continent. Simulations examining the impact of the continental scale sympatry assumption on the statistical power of 2-regime MC models demonstrate that, even for relatively small clades, large but biologically plausible latitudinal differences in the effect of competition should be detectable, even when sympatry is overestimated (S8 Fig). Nevertheless, there is evidence that this assumption can impact the power to detect subtle differences between regions and for smaller trees, the estimated direction of the difference (S8 Fig). However, restricting our empirical analyses to large clades (N ≧ 100), we still find no support for a consistently stronger impact of competition on phenotypic evolution in tropical lineages (S8 Table).

Latitudinal variation in tempo of phenotypic evolution

Evolutionary rates estimated from single-rate models did not vary according to clade-level index of tropicality (Figs 2 and S9 and S13 Table). Similarly, estimates of rates from latitudinal models were neither consistently lower nor higher in tropical regions (Figs 3D and S10 and S14 Table). We did find lower rates of locomotion pPC3 (Figs 3D and S10 and S14 Table) and bill pPC2 evolution in tropical lineages (S10 Fig and S14 Table), but the difference between regions was small, and the overall strength of this relationship was weak. Observational error contributed to these patterns: We found a significant negative correlation between observational error and the clade-level index of tropicality for body mass (S11 Fig and S15 Table), and we also found that there was a correlation between rates of body mass and locomotion pPC3 evolution in standard single-regime BM models excluding error (S12 Fig and S16 Table) and that the magnitude of the difference between tropical and temperate rates of trait evolution was higher in analyses of 2-regime fits excluding error (S12 Fig and S17 Table).

Predictors of support for models with an effect of competition on phenotypic evolution

We found no evidence that territoriality or diet specialization are useful predictors of support for models that incorporate the impact of co-occurring species on phenotypic evolution (S18 Table). We did, however, find that the maximum proportion of species co-occurring on a continent (i.e., the maximum number of extant lineages on a single continent divided by the total clade size) had a pronounced impact on model selection—clades with a high proportion of lineages occurring on the same continent were more likely to be best fit by the MC model, whereas clades with a low proportion of co-occurring lineages were more likely to be best fit by the exponential DD model (S13 and S14 Figs and S18 Table). In addition, we found that the MC model was less likely to be favored in clades with many members living in single-strata habitats (S18 Table).

Two-regime phylogenetic models of phenotypic evolution

One approach to analyze gradients in phenotypic evolution is to fit phylogenetic models of phenotypic evolution that allow model parameters (e.g., evolutionary rates) to vary across the phylogeny; such models are already available for the simplest models of trait evolution such as BM and Ornstein–Uhlenbeck (OU) models [68,69]. To explore the effects of species interactions, we developed further extensions to early burst (EB), DD, and MC models, allowing parameters to be estimated separately in 2 mutually exclusive groups of lineages in a clade. Generalizing these new models to estimate parameters on more than 2 groups, or on non-mutually exclusive groups, is straightforward.

We began by developing a 2-regime version of the EB model in which rates of trait evolution decline according to an exponential function of time passed since the root of the tree [70]. We used this model here to ensure that the DD models, which incorporate changes in the number of reconstructed lineages through time, are not erroneously favored because they accommodate an overall pattern of declining rates through time. To estimate rates of decline separately for mutually exclusive groups, we formulated a 2-regime EB model with 4 parameters (Table 1): z₀ (the state at the root), σ₀² (the evolutionary rate parameter at the root of the tree), r_A (controlling the time dependence on the rate of trait evolution in regime “A”), and r_B (time dependence in regime “B”). This model can be written as follows: (Eq 1) where is the trait value of lineage j at time t, and dW_t represents the BM process (S1 Fig). This model is the 2-regime equivalent of the EB model, where σ²(t) = σ₀²e^rt; the (1/2) factor in Eq 1 comes from taking the square root of the rate.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Parameters of models used in analyses.

https://doi.org/10.1371/journal.pbio.3001270.t001

DD models represent a process where rates of trait evolution respond to changes in ecological opportunity that result from the emergence of related lineages [19,20]. When the slope of these models is negative, this is interpreted as a niche filling process where rates of trait evolution slow down with the accumulation of lineages. We considered 2 versions of DD models, with either exponential (DD_exp) or linear (DD_lin) dependencies of rates to the number of extant lineages. The 2-regime model has 4 free parameters (Table 1): z₀ (the state at the root), σ² (the evolutionary rate parameter), r_A (the dependence of the rate of trait evolution on lineage diversity in regime “A”), and r_B (diversity dependence in regime “B”). For the exponential case, this model can be written as follows: (Eq 2) for the exponential case, where and are the number of lineages in regime A and B at time t. This model is the 2-regime equivalent of the DD_exp model, where σ²(t) = σ₀²e^rn(t); the (1/2) factor in Eq 2 comes from taking the square root of the rate. For the linear case, this can be written as follows: (Eq 3) This model is the 2-regime equivalent of the DD_lin model, where σ²(t) = σ₀² + bn_t and b denotes the slope in the linear model. Standard DD models ignore whether lineages coexist, yet only those lineages likely to encounter one another in sympatry are able to compete with one another. Thus, we extended our model to incorporate ancestral biogeographic reconstructions to identify which species interactions are possible at any given point in time (i.e., which species co-occur [23]). With biogeography, these become (Eq 4) for the exponential case and (Eq 5) for the linear case, where A is a matrix denoting biogeographical overlap, such that A_j,l = 1 if lineages j and l coexist in sympatry at time t, and 0 otherwise (S1 Fig).

The MC model is a model of competitive divergence [22,43], wherein sympatric lineages are repelled away from one another in trait space. We formulated the 2-regime MC model, which has 4 parameters (Table 1): z₀ (the state at the root), σ² (the evolutionary rate parameter), S_A (the strength of repulsion in regime “A”), and S_B (the strength of repulsion in regime “B”). This model can be written as follows: (Eq 6) Incorporating biogeography, this becomes (Eq 7) We developed inference tools for fitting the 2-regime MC and DD models to comparative trait data, following the numerical integration approach used previously [43,56]. For the EB model, we developed a branch transformation approach similar to the one used in mvMORPH [71]. In all model fits, we incorporated the possibility to account for deviations between measured and modeled mean trait values for each species [72–74] (see S1 Text for details). These deviations are of 2 types: the “known” deviation associated with estimating species means from a finite sample and the “unknown” deviation linked to intraspecific variability unrelated to the trait model (e.g., instrument errors and phenotypic plasticity). We follow the common practice of lumping these 2 sources of deviations (often called “measurement error”) and referring to them as “observational error.” A simulation study demonstrated the reliability of estimates using these tools (S1 Text and S7 Data). Functions to simulate and fit these phenotypic models are available in the R package RPANDA [86].

Phylogeny and trait data

We obtained phylogenies of all available species from birdtree.org [26] and created a maximum clade credibility tree in TreeAnnotator [75] based on 1,000 samples from the posterior distribution (S13 and S14 Datas). Since the MC and DD models require highly sampled clades [43], we used the complete phylogeny including species placed based on taxonomic data [26] and the backbone provided by Hackett and colleagues [76]. We then extracted trees for each terrestrial (i.e., non-pelagic) family with at least 10 members (n = 108). As island species are generally not sympatric with many other members of their families (median latitudinal range of insular taxa = 1.28° and non-insular taxa = 15.27°), we further restricted our analyses to continental taxa, excluding island endemics and species with ranges that are remote from continental land masses. We gathered data on the contemporary ranges of each species from shapefiles [77].

Mass data were compiled from EltonTraits [78] (n = 9,442). In addition, we used a global dataset based on measurements of live birds and museum specimens [48] to compile 6 linear morphological measurements: bill length (culmen length), width, and depth (n = 9,388, mean = 4.5 individuals per species), as well as wing, tarsus, and tail length (n = 9,393, mean = 5.0 individuals per species). These linear measurements were transformed into pPC axes describing functionally relevant variation in bill shape and locomotory strategies (S1 Text and S2 and S3 Tables and S1 Data).

Biogeographic data and reconstruction

Phylogenetic models that account for species interactions require identifying lineages that are likely to encounter one another [43]. To discretize the contemporary ranges of each species, we classified them as being present or absent in 11 different global regions [79]: Western Palearctic, Eastern Palearctic, Western Nearctic, Eastern Nearctic, Africa, Madagascar, South America, Central America, India, Southeast Asia, and Papua New Guinea/Australia/New Zealand. To assign each species to the global region(s) they occupied, we used several approaches. As a first pass, we used the maximum and minimum longitude and latitude for species’ (nonbreeding) ranges. When the rectangle formed by these values fell entirely within the limits of a given global region, we assigned that region as the range for the focal species. Next, for species that did not fall entirely into one region, we compiled observation data from eBird.org [80] to identify all of the regions that a species occupies using country-level observations. Finally, for species whose ranges could not be resolved automatically using these techniques, we manually inspected the ranges.

We incorporated estimates of the presence/absence of each lineage in each range through time using ancestral range estimation under the dispersal–extinction–cladogenesis (DEC) model of range evolution [81]. We fit DEC models to range data and phylogenies for each family with the R package BioGeoBEARS [81,82]. Since the continents have changed position over the course of the time period of family appearance (clade age range = 12.84 to 71.88 Mya), we ran a stratified analysis with adjacency and dispersal matrices defined for every 10 my time slice [79]. Using the maximum likelihood (ML) parameter estimates for the DEC model, we then created stochastic maps for each family in BioGeoBEARS, each representing a single hypothesis for which ranges each lineage occupied from the root to the tip of the tree.

Tropical and temperate breeding habitats and reconstruction

To investigate the impact of latitude on trait evolution in 2-regime models, we assigned each species to either the “tropical” or “temperate” regime, based on its breeding range (i.e., a species that breeds exclusively in the temperate zones but migrates to the tropics when not breeding is assigned to the temperate zone). We focused on the breeding ranges of all species as they are likely to be the arena of strongest competition over territorial space and food. To do this, we first assigned each species to either “tropical,” “temperate,” or “both” based on breeding range limits extracted from range data in shapefiles and defining the tropics as the region between −23.437° and 23.437° latitude. We then fit a continuous-time reversible Markov model where transitions between all categories were allowed to occur at different rates, using make.simmap in phytools [83] on the MCC tree. We then used the ML transition matrix to create a bank of stochastic maps under this model, each indicating a possible historical reconstruction of tropical versus temperate habitats through time from the root to tips (S1 Fig). In each stochastic map, we collapsed the “both” category and the “temperate” category to compare lineages with exclusively tropical ranges to lineages with breeding ranges that include temperate regions. Therefore, our “tropical” category indicates that a species breeds exclusively in the tropics, and our “temperate” category contains all species with breeding ranges that include the temperate zone (S4 Fig).

We note that this is a relatively simplistic way of categorizing tropical and temperate membership, and we hope that future methods will enable more sophisticated inferences of historical biogeography alongside paleolatitude and/or paleoclimate. However, given the scope of our analyses, and the emerging evidence that many tropical species ranges have shifted over the timescale of this study [84,85], we opted to keep the results of the historical biogeographical inference and the latitudinal regime reconstruction independent. Future extensions may accommodate the development of more sophisticated paleolatitude models, as well as interactions between various abiotic (e.g., global climate fluctuation [57]) and biotic factors.

Accounting for uncertainty in historical biogeography and latitude

We accounted for uncertainty in ancestral reconstructions by fitting phenotypic models on at least 20 stochastic maps of ancestral tropical/temperate range membership (for all 2-regime models) and/or biogeography (for all models incorporating competition, in both single- and 2-regime versions). For the single-regime model fits that included competition (i.e., DD and MC models), we computed model support and parameter estimates as means across fits conducted on stochastic maps of ancestral biogeography. For the 2-regime model fits, we computed model support and parameter estimates as means across fits conducted on stochastic maps of ancestral tropical/temperate range membership. For the 2-regime model fits with competition, these means also account for variation in estimates of ancestral biogeography (S1 Fig).

Given the scope of these analyses, we chose to account for uncertainty in the biogeographic reconstructions and in the ancestral reconstruction of tropical/temperate living while keeping the topology fixed under the MCC tree. A previous study with a similar model fitting approach found that results on MCC trees were highly concordant with results fit to trees sampled from the posterior distribution [56]. Moreover, there is no reason, to our knowledge, why basing inferences on the MCC tree would bias conclusions about latitude in any systematic way.

Latitudinal variation in mode of phenotypic evolution

We tested whether modes of phenotypic evolution varied with latitude in several ways. First, we used “single-regime” models (Table 1), i.e., models that estimate a single set of parameters on the entire phylogeny regardless of whether lineages are tropical or temperate. We tested whether support for each of these single-regime models varied according to a clade-level index of tropicality (i.e., the proportion of species in each clade with exclusively tropical breeding ranges). Second, we used our newly developed “2-regime” models (Table 1) with distinct sets of parameters for tropical and temperate species and tested whether these latitudinal models were better supported than models without latitude.

We used ML optimization to fit several “single-regime” models of trait evolution to the 7 morphological trait values described above. For all families, we fitted a set of 6 previously described models [43] that include 3 models (BM, OU, and EB) of independent evolution across lineages, implemented in the R-package mvMORPH [71], and 3 further models (DD_exp, DD_lin, and MC) that incorporate competition and biogeography, implemented in the R-package RPANDA [86]. For details of reconstruction of ancestral biogeography, see S1 Text. In the DD models, the slope parameters can be either positive or negative, meaning that species diversity could itself accelerate trait evolution (positive diversity dependence), with increasing species richness driving an ever-changing adaptive landscape [4,67], or, alternatively, increasing species diversity could drive a concomitant decrease in evolutionary rates (negative diversity dependence), as might be expected if increases in species richness correspond to a decrease in ecological opportunity [87].

In cases where families were too large to fit because of computational limits for the MC model (>200 spp., n = 13), we identified subclades to which we could fit the full set of models using a slicing algorithm to isolate smaller subtrees within large families. To generate subtrees, we slid from the root of the tree toward the tips, cutting at each small interval (0.1 Myr) until all resulting clades had fewer than 200 tips. We then collected all resulting subclades and fitted the models separately for each subclade with 10 or more species separately, resulting in an additional 28 clades (n = 136 total).

In addition to this set of models, we fitted a second version of each of these models where the parameters were estimated separately for lineages with exclusively tropical distributions and lineages with ranges that include the temperate region (i.e., “2-regime” models; S1 Text and S2 Fig), limiting our analyses to clades with trait data for more than 10 lineages in each of temperate and tropical regions (S1 Fig; for details of ancestral reconstruction of tropical and temperate habitats, see S1 Text and S4 Fig). The BM and OU versions of these latitudinal models were fit using the functions mvBM and mvOU in the R package mvMORPH [71], and the latitudinal EB, MC, and DD models were fitted with the newly developed functions available in RPANDA [86].

We examined model support in 2 ways. First, we calculated the Akaike weights of individual models [88], as well as the overall support for any model incorporating species interactions and overall support for any 2-regime model. Second, we identified the best-fit model as the model with the lowest small-sample corrected Akaike information criterion (AICc) value, unless a model with fewer parameters had a ΔAICc value <2 [88], in which case we considered the simpler model with the next lowest AICc value to be the best-fitting model.

Latitudinal variation in strength of interactions and tempo of phenotypic evolution

We tested for latitudinal variation in the effect of species interactions on trait evolution using both our single- and 2-regime model fits. With the first class of model, we tested whether parameters that estimate the impact of competition on trait evolution (i.e., the slope parameters of the DD models and the S parameter from the MC model) estimated from our single-regime models varied according to the proportion of lineages in each clade that breed exclusively in the tropics. With the second class of models, we tested whether 2-regime models estimated a larger impact of competition on trait evolution in tropical than in temperate lineages.

Similarly, we tested whether lineages breeding at low latitudes experience lower or higher rates of morphological evolution compared to temperate lineages using our 2 types of models. First, we tested whether rates of morphological evolution varied according to the proportion of lineages in each clade that breed exclusively in the tropics. We estimated this rate directly as the σ² parameter from the single-regime BM model. For the single-regime EB and DD models, we calculated estimates of evolutionary rates at the present from estimates of the rate at the root and the slope parameters. Second, we compared rates estimated separately for tropical and temperate lineages from the 2-regime implementations of the BM, EB, and DD models. We also examined the impact of observational error on rate estimates by fitting single-regime and 2-regime BM models without accounting for observational error.

Examining the potential impact of assuming continental scale sympatry

Our biogeographical reconstructions add important realism into models of species interactions. Nevertheless, species that occur on the same continent do not necessarily interact with one another. We conducted a simulation analysis to determine how our ability to detect the impact of competition on trait evolution may be impacted by the fact that only a subset of the species occurring in a given continent are actually sympatric.

First, we determined the proportion of species that are sympatric within each continent. We calculated range-wide overlap for all family members that ever coexist on the same continent from BirdLife range maps [77] (S6 Data). We defined sympatry as 20% range overlap according to the Szymkiewicz–Simpson coefficient (i.e., overlap area/min(sp1 area, sp2 area)). We also determined if overall levels of sympatry vary latitudinally; to do so, we subset pairs of taxa whose latitudinal means are separated by less than 25° latitude [36] and calculated the midpoint latitude for each pair.

Next, we conducted a simulation study to determine whether competition unfolding between ‘truly’ sympatric species only (i.e., at a level finer than the course continental scale we employed) would systematically impact the fit (i.e., model selection) or performance (i.e., parameter estimation) of the 2-regime competition (MC) models for which we used continental-level sympatry (as in the empirical analyses). To do this, we selected 3 clades spanning the range of tree sizes, each with some traits best fit by single-regime MC model, but none best fit by 2-regime MC model (Cracidae.0 [N = 50, N_tropical = 38, N_temperate = 12], Nectariniidae.0 [N = 122, N_tropical = 89, N_temperate = 33], Picidae.1 [N = 190, N_tropical = 86, N_temperate = 104]). For each of these clades, we simulated 2 biogeographic scenarios reflecting empirical levels of sympatry (see above). In the first, we downsampled the continental biogeography such that 50% of tropical and 50% of temperate taxa that were estimated to occur in the same continent were sympatric (see S1 Text for more details). In the second scenario, to reflect the observed latitudinal variation in sympatry, we downsampled the continental biogeography such that 33% of tropical and 50% of temperate taxa that were estimated to occur in the same continent were sympatric (see S1 Text for more details).

With these downsampled biogeographic histories, representing hypothetical range overlap that is more realistic than our continental-level assumption of sympatry, we simulated trait evolution under the 2-regime MC model. For each clade, we used the mean σ² value estimated under the single-regime MC model in empirical fits of a trait that was best fit by the single-regime MC model. We then varied the ratio of the S_tropical:S_temperate within the range of values in other trait-by-clade combinations where the 2-regime MC model was the best-fit model (S12 Table). For each clade, parameter combination, and downsampled biogeographic scenario, we simulated 100 datasets, for a total of 3,000 simulated datasets. Finally, we fit the same 12 models that were used in empirical analyses. We conducted model selection to identify the best-fit model for each simulated dataset and assessed whether the estimated ln(|S_tropical|/|S_temperate|) had the sign expected given the simulated ratio of S_tropical:S_temperate (S9 Data).

Predictors of support for models with competition

To identify factors other than latitude which influence whether models with competition were favored by model selection, we examined the impact of habitat (the proportion of species in single-strata habitats), territoriality (the proportion of species with strong territoriality), diet specialization (calculated as the Shannon diversity of diets among species in a clade), clade age, clade richness, and the maximum proportion of species co-occurring on a continent.

Statistical approach

We tested for an impact of the proportion of species in a clade that breed exclusively in the tropics on model support and parameter estimates in single-regime models by conducting phylogenetic generalized least squares (PGLS) using the pgls function in the R package caper [89], estimating phylogenetic signal (λ) using ML optimization, constraining values to 0 ≤ λ ≤ 1. We tested support for the 2-regime versions of each model type (BM, OU, EB, DD, and MC) across families for a given trait by fitting intercept-only PGLS models with support for latitudinal models as the response variable. We conducted similar analyses to test overall support for latitudinal models across families for each trait and for differences in parameter estimates for tropical and temperate taxa. We found that statistical support for models incorporating competition was relatively rare in small clades (S6 Fig). As this pattern could be related to lower statistical power in smaller datasets [43], we focused all analyses of evolutionary mode (i.e., model support and parameter estimates from models incorporating competition) on clades with at least 50 species (n = 66 for single-regime fits and n = 59 for 2-regime fits).

For analyses of predictors of support for models with competition, we used the R package MCMCglmm [90] to fit phylogenetic generalized linear mixed models with categorical response variables indicating whether MC or DD_exp models were chosen as the best-fit model (S12 Data).