Dynamical analysis of the global elevational diversity gradient in passerine birds reveals a prominent role for highlands as species pumps

Low elevation regions harbor the majority of the world’s species diversity compared to high elevation areas. This global elevational diversity gradient, suggests that lowland species have had more time to diversify, or that net diversication rates have been higher in the lowlands (either due to higher ecological limits or intrinsically higher diversication rates). However, highlands seem to be cradles of diversity as they contain many young endemics, suggesting that their rates of speciation are exceptionally fast. Here, we use a phylogenetic diversication model that accounts for the dispersal of species between different elevations to examine the evolutionary dynamics of the elevational diversity gradient in passerine birds, a group that has radiated globally to occupy almost all elevations and latitudes. We nd strong support for a model where passerines diversify at the same rate in the highlands and the lowlands but where the rate of dispersal from high to low elevations is more than twice as fast as in the reverse direction. This suggests that while there is no consistent trend in diversication across elevations, highland regions act as species pumps because the diversity they generate migrates into the lowlands, thus setting up the observed gradient in passerine diversity. This species pump is particularly strong in the tropics, where the inferred rate of speciation is 1.4 times faster than in the temperate zone. We conclude that despite their lower diversity, highland regions are disproportionally important for maintaining diversity in the adjacent lowlands. The extinction of species in the tropical highlands due to rapid climate change this century could thus have major and long-lasting impacts on global passerine diversity.


Introduction
The striking differences in the form and diversity of life as one travels up mountains is one of the most prominent and long recognized patterns in biogeography (Lomolino 2001). Across animals and plants, and across mountains globally, diversity tends to peak at low or intermediate elevations and then declines towards mountain summits (McCain 2010;Quintero and Jetz 2018). Thus, while mountain regions have been regarded as hotspots of diversity (Orme et al. 2005;Rahbek et al. 2019b) at a global scale, most species occur in the lowlands with relatively few species in the highlands. The possible evolutionary origins of this global elevational diversity gradient (EDG) can be divided into two nonmutually exclusive explanations. One possibility is that for most organism groups, there has simply been more time for species to accumulate at low to intermediate elevations, with highland environments only colonized more recently (Wiens et al. 2007). Another possibility is that net rates of diversi cation are faster in the lowlands perhaps because the greater area or energy availability increases the ecological limit to diversity or promotes faster intrinsic rates of speciation and lower extinction (Pigot et al. 2016).
Yet, despite the higher diversity of species in the lowlands, there is also evidence that highland regions may be cradles of species diversity, characterized by exceptionally fast rates of speciation (Körner and Spehn 2002;Merckx et al. 2015). In other words, speciation may be fastest at high elevations, where species diversity is currently lowest (Quintero and Jetz 2018).
Macroecological studies of current species distribution patterns have long supported the idea that highlands may be hotbeds of evolution, especially in the tropics (Fjeldsa 1994). Large numbers of endemic species are concentrated on tropical mountains, far more than would be expected by chance or current climate (Jetz et al. 2004). More recent phylogenetic studies have shown that many of these endemic species are comparatively young, pointing towards rapid speciation (Hughes and Eastwood 2006;Weir 2006). There are various reasons why rates of speciation may be faster at high elevations, including the fragmentation of habitats on different mountain summits (Hughes and Eastwood 2006), ecological opportunity as new habitats were made available during recent mountain uplift (Cozzarolo et al. 2019), exposure to higher levels of ultraviolet radiation boosting rates of mutation (Davies et al. 2004), and the susceptibility of species inhabiting narrow thermal bands to become isolated by the expansion of glaciers during recurring ice ages (Graves 1988;Weir 2006). All these factors are thought to be particularly important in the tropics because limited seasonal variation enables the greater thermal strati cation of species across mountain slopes, enhancing opportunities for geographic isolation and ecological selection (Janzen 1967). However, while some studies have supported the idea that diversi cation is faster in the highlands (Cai et al. 2018), others studies have found no evidence for differences in diversi cation rates across elevations (Rana et al. 2019), and thus whether rates of speciation vary consistently across elevation remains debated.
One explanation for this lack of resolution is that the dynamics of the EDG may be highly variable, differing between mountain regions with contrasting geological histories and geographic properties (Rahbek et al. 2019a). For instance, some highland lineages in the Andes have undergone explosive radiations, but such cases appear to be less common in mountain systems in the Afrotropics (Rahbek et al. 2019a). In addition, gradients in diversity -as well the phylogenetic branching times used to model these gradients -arise not only from differences in speciation and extinction but also the dispersal of lineages across elevations. Yet, models accounting for both differences in diversi cation as well as the movement of species between elevational bands have rarely been applied. Some verbal models predict that highland communities are evolutionary sinks, derived primarily from the lowland taxa either invading or being passively transported to higher elevations during mountain uplift (Ribas et al. 2007). Other verbal models predict that highland regions act as species pumps, with lineages arising at higher elevations moving downslope (Bates and Zink 1994;Roy 1997;Garcia-Moreno et al. 1998;Ribas et al. 2007;Schwery et al. 2015), providing an important contributor to the greater diversity of the adjacent lowlands, and leading to high species richness at the ecotone between mountains and lowlands, particularly in the tropics (Herzog et al. 2005;McCain 2010;Quintero and Jetz 2018). Given this complexity, establishing the evolutionary dynamics of the EDG requires developing models that integrate the processes of diversi cation and dispersal while also accounting for potential differences in diversi cation dynamics across species, latitudes and regions.
Here, we study the dynamics of the EDG of passerine birds, a global 'super radiation' of exceptional diversity, including approximately 5700 species, found at almost all elevations, latitudes and biogeographic realms (other than Antarctica). Passerines are a relatively recent radiation, having arisen in approximately the last 47 million years (Oliveros et al. 2019), and represent an ideal study group.
Passerines show a clear EDG, with over 8 times as many species present in the lowlands (n = 5095, lower elevation bound < 1500 m) as there are mid (> 1500 m) or high (> 3000 m) elevation specialists combined (n = 611) (Figure 1; Extended Data Table 1). The sheer diversity of passerines, combined with the availability of phylogenetic trees (Jetz et al. 2012), provides substantial power for comparing competing scenarios of diversi cation and dispersal between different elevation bands ( Figure 1). Importantly, although characterizing the distribution of species in topographically diverse regions is challenging, the elevational ranges of passerines are relatively well known, with estimates of upper and lower elevation limits available for almost all (96%) species.
Using this dataset, we rst examine the relationship across passerines between the elevational state of a species and its speciation rate using the tip-DR metric, which has previously been applied to test for latitudinal and elevation gradients in avian speciation rates (Quintero and Jetz 2018). Because this metric does not account for the movement of species between elevational bands, we then apply a recently developed dynamic phylogenetic model (secsse, Herrera-Alsina et al. 2019), that in addition to testing for elevation-dependent rates of speciation and extinction, also enables us to account for the transition of species between different elevational states over evolutionary time. Using this dynamic model, we are thus also able test whether species movement downhill is faster or slower than the rate of movement uphill, and thus whether highlands act a source or sink of species diversity respectively. Finally, because there is evidence that rates of speciation in passerines may vary between the tropics and temperate zone (Weir and Schluter 2007) and between the New and Old World (Jetz et al. 2012) we allow rates of speciation (or extinction) to vary between different longitudinal (Old and New World) or latitudinal (tropics and temperate) regions, when testing for an effect of elevation ( Figure 2).

Diversi cation and dispersal across elevations
When comparing the tip-DR of species currently occurring in the lowlands and the highlands, rates of speciation are inferred to be slightly but signi cantly faster in the highlands. This is consistent with previous assemblage-level analyses, indicating that the average tip-DR of highland bird assemblages is greater than that of lowland bird assemblages (Quintero and Jetz 2018). However, closely related species show a strong tendency to share similar tip-DR values (Pagel's λ = 0.987, 95% CI: 0.984-0.99), and thus any association between elevation and speciation rate could be driven by phylogenetic nonindependence. In accordance with this, when we accounted for shared ancestry using a phylogenetic generalized linear model (PGLS), the relationship between elevation and tip-DR disappeared (Extended Data Table 2).
While this result suggests that there is no association between elevation and speciation rate, this analysis must be interpreted with caution because it assumes that species distributions are static and thus does not account for the movement of species between different elevational bands over evolutionary time. To address this shortcoming, we next tted a dynamic model that estimates the speciation (or extinction rate) associated with each elevation state as well as the transition rates between these states. We compared the t of a model in which rates of speciation (or extinction) vary across elevation to two alternative null models. First, a constant rate model, in which all lineages share an identical rate of speciation (and extinction). Second, a concealed-trait dependent model, in which rates of speciation (or extinction) are allowed to vary across lineages due to another 'hidden' trait, but they do so independently of elevation. This second null model is more realistic because it accounts for the possibility that rates of diversi cation have been heterogeneous across passerines, and means that support for elevationdependent rates of speciation (or extinction) does not arise simply because of the unrealistic assumption of a constant rate null model applied to such a large clade.
When compared to a standard constant rate (CR) null model, in which rates of speciation are equal across elevational states, an Elevation-dependent (ED) model, in which rates of speciation increase with elevation, was more strongly supported. This result seemingly supports previous analyses suggesting faster rates of bird speciation in the highlands (Quintero and Jetz 2018). However, when we applied a concealed trait-dependent (CTD) null model, in which rates of speciation are allowed to vary across lineages but independently of elevation, the CTD null model was overwhelmingly supported (AICw 0 .99). Thus, while our results provide substantial support for a scenario in which speciation rates vary across lineages, this variation is likely due to other factors than elevation (Extended Data Table 3). The best supported model was a latitudinal-CTD model, in which rates of speciation were faster in the tropics than the temperate zone. Across all models, extinction was estimated to be low (latitudinal-CTD model, μ = 0.0001/myr) and a latitudinal-CTD model with heterogeneity in rates of speciation, was more strongly supported than a CTD model with heterogeneity in rates of extinction across lineages and/or regions. Finally, we found no evidence for differences in rates of speciation between the Old and New World (Extended Data Table 3).
While we found no effect of elevation on rates of speciation or extinction, we found substantial differences in rates of downhill and uphill dispersal. According to the best model, the estimated per lineage rate of dispersal uphill is 0.072/myr and the rate of downhill dispersal is 0.189/myr ( Figure 2). This model of faster downhill dispersal was much more strongly supported than an alternative expansion-contraction scenario, in which rates of elevational range expansion and contraction may differ, but there is no difference between uphill and downhill movement ( Figure 2). The best-supported model allowed transition rates between elevational states to differ from those of the transition rates between the concealed trait states, rather than assuming these are identical (Extended Data Table 3), even though the former model has more parameters. We found, as may be expected, that lineages disperse between continents at rates that are smaller than dispersal across elevational bands. The best model shows that the rate of lineage exchange between the tropics and temperate zone is 0.0187/myr with lineage dispersal between the Old Word and New World occurring at a rate of 0.00016/myr ( Figure 2). Note that all these rates are per-lineage rates.
Geographic origin of passerines and accumulation of lineages over time Because the elevational origin of songbirds can have an important in uence on the current gradient of species richness, we extracted the probabilities (using the best supported model) of each state (i.e., the combination of elevation and latitude) at the most-basal node of the phylogeny to estimate where the clade rst appeared. These probabilities indicate the most likely state of the ancestral species just before it split at the crown of the phylogeny, and so the true origin of the clade (at the stem age) may have been different. Regardless of whether rates of speciation are allowed to vary with latitude or longitudinal region, a highland origin of passerines was equally well supported as a lowland origin. This lack of strong support for either elevational state as the origin for passerines re ects the relatively rapid transitions across elevation inferred by our analysis and which have likely erased this historical signal.

Discussion
In spite of clear differences in passerine diversity at different elevations, we nd no differences in speciation or extinction rates across elevational zones. Instead, the best supported scenario is a Latitudinal-Concealed-Trait Dependent model with higher rates of downhill dispersal than uphill dispersal. This CTD model indicates that there is substantial variation in speciation rates across lineages and latitude, but that this variation is not related to elevation. Thus, while previous more taxonomically or geographically focused studies have found evidence for differences in diversi cation across elevation, this is not supported by our global analysis across all passerines.
The absence of a consistent effect of elevation on diversi cation is unlikely to be explained by a lack of statistical power. Our analysis contains thousands of species and our model did detect differences in diversi cation rate between the temperate and tropical zone. This suggests that if there was a globally coherent effect of elevation on speciation or extinction, this would also have been detected in our analysis. Furthermore, the absence of any effect of elevation on diversi cation cannot be explained by a failure of our model to detect ner-scale heterogeneity in speciation rates within the tropics or temperate zone. The strong support for the Concealed Trait-Dependent model con rms that there is indeed substantial heterogeneity in rates of speciation across passerines within these regions, but that this variation is not aligned with elevation. Critically, our results show that failure to account for this background heterogeneity in rates of speciation would have led to the spurious conclusion that rates of speciation do increase with elevation. We therefore suggest that application of such a Concealed Trait-Dependent model is an important advance that should also be applied when examining other putative drivers of diversi cation rates (e.g. body size).
One possible explanation for why our global analysis does not detect an elevational gradient in speciation or extinction, is that different dynamics may prevail in different regions obscuring the overall importance of elevation. Differences across latitude may be particularly important in this regard, because it is mainly highlands in the tropics where rates of speciation are thought to be promoted (Fjeldså et al. 2012). However, this is unlikely to explain our results because most of the species in our analysis (85%) occur in the tropics, and thus it is the dynamics of the EDG in this region which dominate our model inferences. Instead, our results suggest that net rates of diversi cation are boosted in the tropics, but this occurs regardless of elevation.
Another possibility is that although diversi cation and elevation may be related in some mountain systems, these effects of elevation are highly context dependent, determined by the particular geological history and environment in different mountains (Quintero and Jetz 2018). Accounting for such variability is challenging because it would require a substantial increase in model complexity and the number of parameters that need to be inferred. We note that support for such a complex model would not alter our main conclusion that, at a global scale, there is no consistent trend in diversi cation across elevation.
Our nding that net rates of speciation in passerines do not vary consistently with elevation contrasts with that of a recent study showing that speciation rates increase with elevation (Quintero and Jetz 2018). We note a number of important methodological differences that may explain these contrasting ndings. Quintero and Jetz (2018) calculated speciation rate for each extant passerine species (tip-DR) and then averaged this across all the species inhabiting each elevational band across major mountain systems. This approach strongly differs from ours, because it mainly captures recent speciation and because it links the speciation rate of a species to its current elevation, thus ignoring dispersal across elevational bands and the possibility that the speciation events leading to a particular lineage occurred at a different elevation to where that species currently resides. Furthermore, while the shared phylogenetic history of each species is automatically accommodated in our dynamic model (i.e. each branch of the tree is only counted once), our results show that treating each species as statistically independent when analyzing tip-DR may lead to spurious results. On the basis of these results, we suggest that accounting for the biogeographic dynamics and phylogenetic non-independence of lineages is likely to be critical when making inferences of how rates of speciation or extinction vary across elevation or any other gradient.
Many of the world's current highland regions did not develop until relatively recently in geological history suggesting that lowland areas may have had a head start in accumulating species. The Andes did not undergo signi cant uplift until ~15 mya (Yin and Harrison 2000), the Himalayan-Tibetan orogenic system's higher elevations are of similar age (Yin and Harrison 2000) and the New Guinea Highlands built up to current-elevation levels at ~5 mya (van Ufford and Cloos 2005). Many low-elevation regions are much older and more climatically stable, including those associated with cratonic systems like the Guiana Shield (Lujan and Armbruster 2011) and the large pan-African craton (Rino et al. 2008).
Much of the research on the evolutionary dynamics of the EDG has focused on explaining patterns at relatively ne taxonomic scales or particular mountain systems. While these have provided evidence for both uphill (Chazot et al. 2016) and downhill (Elias et al. 2009) movement of species, they have often supported the idea that most clades originated in the lowlands and then invaded higher elevations (McGuire et al. 2007). For example among avian genera, Leptopogon ycatchers (Bates and Zink 1994), Chat-tyrants (Garcia-Moreno et al. 1998), Andropadus greenbuls (Roy 1997) and Thamnophilus antshrikes (Brum eld and Edwards 2007) show a repeated movement of species into the highlands followed by in situ diversi cation. Larger and older clades, however, paint a more complex picture, with McGuire et al (2007) providing evidence that hummingbirds (Trochilidae) have undergone multiple colonizations of Andes from the lowlands, but also the reverse. Our ndings across passerines, place these previous genus-and family-level studies in a broader context, showing that while there are many cases of uphill dispersal followed by radiation, overall, on a per-lineage basis, the dispersal from high to low elevations, occurs at a much faster rate.
Why per-lineage rates of dispersal from highlands to lowlands are faster than the reverse remains unclear. Indeed, most studies have suggested a net movement of species uphill, either because of ecological opportunity (Cozzarolo et al. 2019) or because species are passively transported to higher elevations with their preferred habitats during mountain uplift (Ribas et al. 2007). We suggest that on a global scale the faster transition to the lowlands may be explained by the contrasting geographic structure of highland and lowland regions. Highlands cover a much smaller area than the lowlands, are fragmented or form narrow mountain chains. As a result, almost all highland species occur in close geographic proximity to the adjacent lowlands and thus have the potential to colonize these habitats. In contrast, many lowland species have geographic ranges that are remote from any highland region and thus have no opportunity to colonize these habitats. Future studies could test this hypothesis by extending our modelling framework to include an additional geographic state indicating whether lowland species are close or distant from the nearest highland. We note, however, that we are already close to the current computational limit: run-times of the analyses (even when performed in parallel on a computer cluster) were in the order of weeks.
It is important to note that we are reporting the rates of successful dispersal events between elevational bands and that dispersal events not leading to establishment are not recorded. Hence, another explanation for higher rates of dispersal into the lowlands than into the highlands, is that high elevation habitats may be less conducive to colonization over macro-evolutionary timescales (Schumm et al. 2020), either because of a lower availability or resources (Pigot et al. 2016) or because species exhibit strong ecological and physiological constraints that prevent adaptation to colder climates ). In other words, it is possible that lineages are equally likely to disperse regardless the elevation they arise from, but those elevational bands where the niche space is wider, are more likely to accommodate colonizers. Thus, our results are potentially consistent with the idea that higher ecological limits to coexistence in the lowlands ultimately underlie the EDG in passerine birds ).
Our results provide support that highland regions have acted as a species pump for passerine birds, not because of faster total rates of diversi cation (although they may have higher intrinsic rates), but because species arising in the highlands disperse downhill. Thus while highlands constitute only a fraction of passerine diversity they play a disproportionate role in boosting the diversity of the lowlands.
Highland species are at particular risk of extinction from anthropogenic warming this century due to a lack of available cooler habitats to which species can disperse (Freeman et al. 2018). Our results suggest that the loss of highland lineages will also have long term impacts on the diversity of the lowlands and passerines globally.

Phylogenetic framework
We used the Bayesian pseudo-posterior distribution of time-calibrated phylogenies provided by (Jetz et al. 2012), which includes 9993 species of the world's birds out of a total of 10473 species. Using the program TreeAnnotator from the BEAST2 package v2.4.2 (Bouckaert et al. 2014) we produced a maximum-clade credibility tree from all available stage 2 trees with the Hackett backbone. We pruned this tree to the level of Passeriformes (5966 species). We are aware of the shortcomings of using a megaphylogeny that does not include sequence data for each taxon (Braun et al. 2019), but the statistical power achieved by the large number of tips and branching events compensates, at least in part, for the possible lack of taxonomic precision.

Elevation data and large-scale realms
We compiled elevation data for passerines, recording lower and upper elevation bounds of their distribution, based on descriptions in the Handbook of the Birds of the World (del Hoyo et al. 2016). We did not include occasional records at extreme elevation, and elevational distributions are based on breeding range, thus excluding wintering and transient elevation records. Species without known elevational distribution were assigned NA in the data set (n = 260) rather than removed, because this would bias diversi cation rates, and our analysis can account for this missing data (see below).
Our model of elevation-dependent diversi cation and dispersal requires treating elevation as a categorical state. To do this we de ned 3 elevational bands: lowlands (from sea level up to 1500 meter above sea level), mid-elevations (1500 -3000 m), and highlands (> 3000 m). These categories broadly agree with those established by (Chapman 1926) for Neotropical montane birds based on dominant vegetation zones associated with the tropics, subtropics, and alpine zones, respectively. Species may inhabit multiple bands (indeed some species span the entire elevation gradient) and accounting for this variation in elevational range size is important when determining the dynamics of the EDG. We therefore de ned three additional categories for species inhabiting more than one elevational band: low-mid, mid-high and full range (i.e., species living from lowlands to high montane areas). The number of species in each category is shown in Table 1. From hereon we refer to these 6 elevational categories as the elevational state of the species.
We rst performed a global analysis in which rates of diversi cation rates depend on elevation, resulting in 6 states. We then accounted for the possibility that rates of diversi cation may differ across latitude and longitude in the following way. We distinguished species occurring in the tropics from species in the temperate (latitudinal analysis), so that latitude and elevation can differentially affect diversi cation resulting in a model with 12 states (e.g., a species could be a tropical lowland species, or in a temperate mid-high state). We classi ed species as tropical when most of the latitudinal span of their breeding distribution (del Hoyo et al. 2016) lies between tropics of Cancer and Capricorn. Finally, we distinguished Old World from New World species (longitudinal analysis, again resulting in 12 possible states (e.g. Old World low-mid state), in which longitude and elevation can differentially affect diversi cation. We did not perform an analysis where diversi cation rates depend on elevation, latitude and longitude, because the large state space required for such an analysis was computationally unfeasible and numerically unstable given the size of our phylogeny (but see below).

State-dependent diversi cation analysis and parameterization
We used the SSE framework (State-dependent Speciation and Extinction) which allows determining whether diversi cation rates of a clade are associated with an evolving trait (Maddison et al. 2007). In this model, the speciation rate (λ i ) or extinction rate (μ i ) of a lineage depends on its trait state i (here elevation, or a combination of elevation and latitude (tropical or temperate) or longitude (New World or Old World)). In order to keep the number of estimated parameters as low as possible during the likelihood optimization, we only optimized the speciation (or extinction) rate of low, mid and high elevation and used these values to average the states that are a combination of those elevations (e.g. the low and mid elevation rates of speciation are averaged to yield the rate of low-mid elevation).
The state of a species is not static, with species switching to a different state at rate q ij , where i and j represent the state of origin and the state of destination, respectively. This allows us to use trait and branching patterns simultaneously to study macroevolutionary dynamics. In other words, when lineages living in different elevation states experience different speciation/extinction regimes, the shift from one elevational state to another will have an effect on diversi cation rates. Statistical support for elevation affecting diversi cation rates is found when the likelihood of a model where speciation (or extinction) differs across elevation states (Elevation-Dependent model, ED) is higher -after correcting for differences in numbers of parameters -than a model where rates depend on an unknown (hidden or concealed) trait (Concealed Trait-Dependent, CTD) and a model with constant rates (CR) (Beaulieu and O'Meara 2016;Herrera-Alsina et al. 2019). The comparison in terms of likelihood between ED and CTD models is important to prevent spurious conclusions on associations between heterogeneity in diversi cation rates across lineages and the evolution of the trait. We used the R package secsse (Several Examined and Concealed States-Dependent Speciation and Extinction; (Herrera-Alsina et al. 2019) which computes and optimizes the likelihood of the model with 2 or more states.
We assume that species transition between elevation states via the expansion into an adjacent elevation band (i.e., from low to low-mid or from high to mid-high) or from contraction at the edge of the range (i.e. full range to low-mid or mid-high). In other words, we do not allow disjunct elevational ranges as these are rarely if ever observed (Pigot et al. 2016). We considered two alternative scenarios for how rates of transitions between states may vary. In both cases we estimate two transition rates. First, under an updown model, uphill transitions have a different rate than "downhill" transitions. Note that under the updown scenario, changing from low-mid elevation to low elevation means going downhill. Second, under an expansion-contraction model, all expansion rates (i.e., change from one single band to two, or from two to three) are equal but different from contraction rates (i.e., change from two bands to a single band, or from three to two bands) which are also equal. We implemented the expansion-contraction model because it allows us to test a scenario where the rate at which species colonize or become extinct at a particular elevation can differ, but where these rates are independent of elevational direction. In other words, in contrast to the up-down model, colonizing (or becoming extinct at) a lower or higher elevation band is equally likely. As recommended by Herrera-Alsina et al. (2019), for the concealed trait we assumed an identical model structure, in terms of the number of states and possible transitions between these. We implemented two versions of the CTD model: one where the transition rates between concealed states are the same as those between elevation states, and one where we relax this assumption and allow the transition rates of the concealed trait to be different from the transition rates between elevational bands. Because differences in diversi cation across states could be due to either differential speciation or extinction rates, we tested both cases for all the models.
For the longitudinal and latitudinal analyses, in addition to allowing different rates of speciation (or extinction) across elevational bands, we also allow the overall rate of speciation (or extinction) to differ between regions (i.e. tropical vs. temperate regions, Old World vs. New World). We did so by multiplying the rates in one region by a constant factor to give rates in the other region and this factor was optimized. To avoid models with many parameters we did not consider the more complex scenario where the elevation-dependence in rates of speciation (or extinction) varies between regions. We assumed that the transition of species between regions (region exchange) happens to and from the same elevational band (low, mid and high only). While we assume that the transition rates between regions is the same for all elevational bands, these rates are different from the rates of transition between elevational bands within regions.
To prevent nding only local optima during the likelihood optimization, we used ve different initial parameter sets. The rst set of parameters were the estimates of speciation and extinction from a birthdeath model t to the branching times and transition rate was assumed to be a fth of speciation rate. For the second set we doubled the speciation rates of the rst set, and halved the transition rates. In the third, we halved the speciation rates of the rst set and doubled the transition rates. Similarly, the fourth had doubled extinction rates and halved transition rates, and the fth had halved extinction rates and doubled transition rates compared to the rst set. The highest likelihood of the ve starting points was taken as the global optimum and used to compare across models. We used AIC weights -thus penalizing the number of free parameters -to select the best models per analysis.
Our global, latitudinal and longitudinal analyses differ in their assumptions on what factors diversi cation rates depend on (elevation only, elevation + latitude, elevation + longitude). Using only the data necessary to study these dependencies would prevent model comparison, because the data sets would differ. We had six states in the global analysis (as there are six elevational bands) whereas the longitudinal analysis has 12 states (six bands in combination with tropical and temperate regions) which are different from the 12 states in latitudinal analysis (six bands in combination with New and Old worlds). Therefore, we made the AIC values comparable by adding an extra likelihood term to the likelihood computed by secsse that covers the transitions not covered in secsse, using the function tDiscrete from the R package geiger (Harmon et al. 2008). That is, for the longitudinal analysis we added the (maximum) loglikelihood of a simple model of transitions in latitude which uses the phylogenetic tree. In this way, the total likelihood of longitudinal analysis incorporates the likelihood of a model with transitions across elevations + longitude (computed by secsse) + model of transitions between latitude (computed by geiger). Note that only the component calculated using secsse handles diversi cation rates simultaneously with state transition. Similarly, for the latitudinal analysis we added the (maximum) loglikelihood of a model of transitions in longitude given the phylogenetic tree, and for the global analysis we added the (maximum) loglikelihood of a model of transitions in both longitude and latitude -which in fact is the sum of the two previous loglikelihoods.
Finally, to provide a more direct comparison with previous studies, we also examined the association between elevation and tip diversi cation rate (tip-DR) using ANOVA and phylogenetic generalized linear model (PGLS) tted in the R package 'caper'. For this analysis we labelled species as either 'highland' (High, Mid-High) or 'lowland' (Low, Low-Mid) and excluded species limited to mid-elevations (Mid) or spanning the entire gradient (Full-Range), resulting in n = 5186 species. We calculated tip-DR using the 'evol_distinct' function of the 'phyloregion' R package (Daru et al. 2020).
Declarations which adds three more categories: low-mid, mid-high and full range (i.e., species living from lowlands to highlands).