Phylogenetic trees often depart from the expectations of stochastic models, exhibiting imbalance in diversification among lineages and slowdowns in the rate of lineage accumulation through time. Such departures have led to a widespread perception that ecological differences among species or adaptation and subsequent niche filling are required to explain patterns of diversification. However, a key element missing from models of diversification is the geographical context of speciation and extinction. In this study, we develop a spatially explicit model of geographic range evolution and cladogenesis, where speciation arises via vicariance or peripatry, and explore the effects of these processes on patterns of diversification. We compare the results with those observed in 41 reconstructed avian trees. Our model shows that nonconstant rates of speciation and extinction are emergent properties of the apportioning of geographic ranges that accompanies speciation. The dynamics of diversification exhibit wide variation, depending on the mode of speciation, tendency for range expansion, and rate of range evolution. By varying these parameters, the model is able to capture many, but not all, of the features exhibited by birth–death trees and extant bird clades. Under scenarios with relatively stable geographic ranges, strong slowdowns in diversification rates are produced, with faster rates of range dynamics leading to constant or accelerating rates of apparent diversification. A peripatric model of speciation with stable ranges also generates highly unbalanced trees typical of bird phylogenies but fails to produce realistic range size distributions among the extant species. Results most similar to those of a birth–death process are reached under a peripatric speciation scenario with highly volatile range dynamics. Taken together, our results demonstrate that considering the geographical context of speciation and extinction provides a more conservative null model of diversification and offers a very different perspective on the phylogenetic patterns expected in the absence of ecology.
Ecological processes have long been thought to be of fundamental importance in explaining the patterns of speciation and extinction among species and through time (Darwin 1859; Simpson 1953; Stanley 1979; Schluter 2000). For instance, competition and adaptation to available niches are thought to underpin the adaptive radiation of Darwin's finches (Grant P.R. and Grant B.R. 2007). However, it is also recognized that there may be a stochastic element to speciation and extinction rates (Raup et al. 1973; Gould et al. 1977), and disentangling the roles of ecological and stochastic processes in species diversification has become a central theme in macroevolution (Nee 2006). Studies invoking an influence of ecology on diversification are therefore often based on the statistical rejection of a purely stochastic model of diversification (Nee 2006).
The pure birth (Yule 1925) and birth–death (Kendall 1948) models have been the dominant null models in both paleontological and phylogenetic studies over the last 3 decades (Raup 1985; Nee et al. 1994; Ricklefs 2007). These models describe the branching of phylogenetic trees in which the probabilities of speciation and extinction are equal across lineages and constant through time (Raup et al. 1973; Gould et al. 1977; Nee 2006). Empirical studies have consistently revealed that observed clades of similar age often show a greater variance in species richness than expected under an equal rates model (i.e., phylogenetic trees are unbalanced) (Slowinski and Guyer 1993; Mooers and Heard 1997; Blum and Francois 2006). This finding has, therefore, often been interpreted as evidence that ecological differences among lineages must underpin variation in speciation and extinction rates (Dial and Marzluff 1989; Slowinski and Guyer 1993; Owens et al. 1999; Ricklefs 2004). Another set of predictions, which arise from the constant rate model, concerns the accumulation of lineages on a reconstructed phylogeny through time. On a logarithmic scale the number of species increases linearly with time under a pure birth model or shows an upturn toward the present under a birth–death model where death is non-zero (Nee et al. 1992; Nee 2006; Ricklefs 2007). In contrast, numerous studies of real reconstructed phylogenies have found that, in many clades, lineage accumulation was initially rapid but has declined toward the present (Nee et al. 1992; Zink and Slowinski 1995; Seehausen 2006; Weir 2006; McPeek 2008; Phillimore and Price 2008; Rabosky and Lovette 2008). Further evidence for temporal heterogeneity in speciation and extinction rates comes from the lack of correlation between clade age and clade richness in many taxa (Rabosky 2009). Although a multitude of processes may explain heterogeneity in diversification rates, this departure from the predictions of constant rate models is often taken as evidence for the saturation of niche space and limits to diversity (Weir 2006; Phillimore and Price 2008).
The use of the birth–death process, however, extends beyond null models, permeating many aspects of macroevolutionary research (Nee 2006). Estimating the history of diversification from neontological phylogenies (i.e., those based on extant lineages alone) requires an underlying model of diversification (Ricklefs 2007) and the birth–death model has often been used to infer speciation, extinction and diversification rates using branching times (Nee 2001; FitzJohn et al. 2009), the age and diversity of the clade (Yule 1925; Nee 2001; Ricklefs 2007), or terminal branch lengths (Weir and Schluter 2007). Finally, the birth–death model is also assumed in some Bayesian phylogenetic reconstruction approaches in which a constant rate process provides a prior distribution for branching rates (Drummond and Rambaut 2007).
Although the simplicity, mathematical tractability, and apparent transparency of the birth–death model have probably contributed to its continued popularity (Nee 2006), a key element missing from current null models of diversification is the geographical context of speciation. For most animal species at least, a key requirement for speciation is the occurrence of geographical isolation and reduction in gene flow between populations (allopatric or peripatric speciation) (Coyne and Orr 2004). The splitting of a lineage in a phylogenetic tree occurs alongside the splitting of its distribution in space. Speciation therefore not only adds to the number of species within a clade but also simultaneously alters the distribution of geographic range sizes across lineages and through time (Gaston 1998; Waldron 2007). This is important because the geographic range size of a lineage is likely to constrain the probabilities of further speciation and extinction (Rosenzweig 1975). All else being equal, larger areas will offer greater opportunities for geographical isolation due to a higher incidence of dispersal barriers, greater habitat heterogeneity, and the limits to gene flow (Rosenzweig 1978; Gavrilets et al. 2000; Kisel and Barraclough 2010). Large ranges also provide a buffer against stochastic or environmentally driven fluctuations in size that may lead to extinction (Jablonski 1995). Here, we develop a model of diversification that explicitly incorporates the processes of geographical isolation and range dynamics underlying the birth and death of lineages and examine the extent to which the inclusion of these processes can account for the “nonrandom” patterns observed in phylogenetic trees. We first explore the phylogenetic and geographic patterns of clades expected under this model and then compare these patterns with those 1) generated under birth–death models and 2) observed in real avian phylogenies.
MATERIAL AND METHODS
To investigate the effects of the geography of speciation and extinction on cladogenesis, we simulated a process of species diversification using spatially explicit models of species' ranges through time (Fig. 1). Our approach builds on earlier simulation studies that have addressed the geography of speciation (Barraclough and Vogler 2000; Phillimore et al. 2008) and the heritability of geographic range size (Waldron 2007). Our model is both neutral and noninteractive in that the dynamics of species ranges are specified from a single distribution (Hubbell 2001) and we assumed no interactions among taxa or limits to the number of potentially coexisting species (MacArthur and Wilson 1967). We regard these conditions as the most sensible and transparent starting point for a geographic model of diversification in which the role of ecology is excluded. Each simulation started with a continuous square domain with side lengths of 500 units containing a single species, represented by a square range with sides of 150 units (9% of the domain). The position of the range was randomly assigned by drawing the coordinates of the range centroid from a uniform distribution and discarding those ranges which crossed the domain boundary (Fig. 1a). We also explored the effects of assuming a larger starting range size of 350×350 units (49% of the domain) (Supplementary Figs. 4 and 5, available from http://www.sysbio.oxfordjournals.org/).
Simulating Range Dynamics
Geographic ranges are dynamic through time, shifting in their size, shape, and position (Ricklefs and Cox 1972; Liow and Stenseth 2007). We simulated these dynamics by adding a random normal deviate to each edge of each rectangular species range at each time step with each simulation lasting for 2000 time steps, allowing the boundaries of the species range to drift independently of one another across the landscape (Barraclough and Vogler 2000; Phillimore et al. 2008) (Fig. 1b). We initially assume that range boundaries fluctuate idiosyncratically through time, and thus in our baseline simulations set the mean of the normal deviate (μ) to 0. We therefore model the location of range boundaries as evolving under a Brownian process, a simplification that may be appropriate if ranges are limited by the independent variation in a large number of environmental and evolutionary factors (Colwell and Lees 2000; Ricklefs 2008). Fluctuations in species ranges can occur rapidly (Davis and Shaw 2001) relative to the lengths of time required for speciation (Weir and Schluter 2007), but there is little empirical evidence on typical rates of range size evolution across clades. We therefore explored arbitrary levels of range dynamism by modifying the variance of the random normal deviate (σ2) from which changes in range boundary position are drawn (values: 0.05, 2.5, 5, 7.5, and 10): when σ2 = 10, 50% of range edges advance or retreat by more than 2.13 spatial units per time unit. We note that although faster rates of range evolution are possible, the effects of range dynamism on the patterns of diversification were evident within the range of parameter values explored. Given the suggested importance or range expansions in determining slowdowns in diversification rate (Phillimore and Price 2009), we also explored scenarios where ranges had a tendency to expand by altering the mean of the random normal deviate from 0.025 to 0.6 (values: 0.025, 0.05, 0.075, 0.1, 0.2, 0.4, 0.6).
Random fluctuations in species range boundaries lead to changes in geographic range size, which evolves according to the product of the summed normal deviates of opposite range boundaries. This may eventually result in the range size of a species drifting to zero, whereupon the species becomes extinct (Fig. 1d). Extinction is therefore an emergent property of the range dynamics and is most likely when μ is low and σ2 is high. Because the dynamics of species' ranges occurred in discrete time, this will tend to underestimate the extinction rate compared with when ranges fluctuate continuously because in the latter case species could have drifted to extinction during the intervening period. To minimize the extent to which extinctions were underestimated, we therefore simulated range movement at a temporal resolution that exploratory simulations showed was sufficient for the number of missed extinctions to stabilize at a low value under the range of parameter values we explored.
The edge of the domain represented a hard boundary, and ranges that extended beyond the domain boundary were immediately truncated (Barraclough and Vogler 2000). We regard this scenario as more closely resembling the geometric constraints on species' ranges than assuming either a torus or an infinite domain (Hubbell 2001). In our simulations, time is not absolute but is instead set by the relative rates of range movement and speciation, and doubling the length of the simulation, for instance, is equivalent to doubling the rate of range movement and speciation. However, we envisage a single time unit as corresponding to a period on the order of 105 years, and other parameters were chosen to give biologically plausible rates of extinction and speciation within this time frame. We did not investigate other, more complex, models of range size evolution (Liow and Stenseth 2007) because theoretical explanations for these dynamics are generally based on nonneutral, ecological processes (Ricklefs and Bermingham 2002).
For most animal species, some form of geographical isolation between populations appears to be required for speciation (allopatric or parapatric speciation) (Coyne and Orr 2004). New species arose in our simulations either by a barrier bisecting the existing distribution of a species (vicariance) or from dispersal to a new location (dispersal or peripatric speciation). Under both scenarios, the probability of speciation arose as an emergent property of the model. Multiple species could speciate in a given time step either through the bisection of multiple species' ranges or through the simultaneous production of dispersal events, with speciation treated as an instantaneous event. Although the assumption of instantaneous speciation is a simplification (Losos and Adler 1995), this approach can be thought of as modeling only those instances where geographical isolation persisted for long enough to cause speciation. Following speciation, the ranges of the daughter species moved independently of one another and could each undergo subsequent speciation events or go extinct (Fig. 1e). Although here we focus on an allopatric model of speciation, we expect that a parapatric scenario, in which reproductive isolation arises in the absence of complete geographic isolation (Coyne and Orr 2004), would produce results qualitatively similar to our vicariance model. This is because both models result in the splitting of an ancestral species range regardless of whether this occurs due to a dispersal barrier (vicariance) or selection across a steep environmental gradient (parapatry). Under both speciation models we discarded simulations in which fewer than 8 species survived to the last time step, equal to the minimum number of species observed within our selection of bird genera.
In the vicariance scenario, we modeled geographic barriers by randomly dropping line segments onto the domain (Rosenzweig 1978) (Fig. 1b). Barriers could miss a species range, pass only part of the way through it, or bisect it, with only the latter resulting in speciation (Rosenzweig 1978) (Fig. 1c). Barriers were dropped on the domain following a waiting time drawn from an exponential distribution (mean = 20 time units), rounded to the nearest time step, consistent with barriers arising via a Poisson process of constant rate per unit time. Barriers were removed immediately so that range dynamics could proceed unimpeded. The coordinates for the centroid of each barrier within the domain were randomly drawn from uniform distributions. The orientation of the barrier was assigned as horizontal or vertical randomly and with equal probability; all ranges consequently retain a rectangular form. Barriers were allowed to extend beyond the domain to minimize the tendency for a mid-domain effect in the frequency of barriers (Colwell and Lees 2000). Although barriers will still tend to overlap more at the center of the domain than at the edges, this is mitigated by the random placement of the species range at the start of the simulation.
The length of the barriers was drawn from a uniform distribution, in the range 0–550 so that species covering the entire domain had a nonzero probability of speciation (P =0.0003). In accordance with previous modeling results (Rosenzweig 1978), the probability of speciation in our vicariance model peaks at intermediate range sizes (Fig. 2). For square ranges, the peak probability (P = 0.005) of speciation occurs when ranges occupy approximately 20% of the domain (Fig. 2a). This is because, as the size of the range increases, species are more likely to encounter barriers but fewer of the barriers will be sufficiently long to completely bisect the range (Fig. 2a) (Rosenzweig 1978). The shape of the relationship between range size and the probability of bisection is sensitive to the distribution of barrier lengths and so we also examined the consequences of using an exponential distribution of barrier-lengths whereby most barriers were short (mean length =110 spatial units) (Supplementary Fig. 1a). In this case, the probability of speciation is reduced, with the peak probability (P = 0.0008) corresponding to a smaller range size (occupying approximately 10% of the domain) (Supplementary Fig. 1a). The probability of range bisection is also modulated by range shape; increasing the length to width ratio of the range increases the chance of being bisected (Fig. 2a). Range shape is particularly important when most barriers are small relative to the size of species' ranges (Supplementary Fig. 1a).
We simulated dispersal speciation such that, at each time step, each extant species could give rise to a new peripatric species. The probability of producing a peripatric species is expected to be linearly related to either the perimeter or the area of the geographic range, depending on the relative range width and the average dispersal distance (MacArthur and Wilson 1967). Because the size and shape of the species' ranges varied in our simulations, for consistency we modeled the probability of dispersal as the range perimeter expressed as a proportion of the domain perimeter. Qualitatively similar probabilities were obtained if the proportional area was used (Supplementary Fig. 1b). When a species was found to have dispersed, a dispersal direction was selected perpendicular to a randomly selected point along the range perimeter, thus favoring dispersal perpendicular to the long axis of a range. The dispersal distance away from the parent range was drawn from an exponential distribution with a mean of 100 and hence a median dispersal distance of approximately 69.3 spatial units. The point identified was used as the centroid of a new species range with edge lengths drawn from a normal distribution (mean = 25, standard deviation = 5). If the colonist occurred near the edge of the domain, we truncated its distribution so that it occurred entirely within the domain; if the new range fell entirely outside the domain, the speciation event was assumed to have failed. Ranges occurring near the edge of the domain have a lower chance of producing peripheral isolates because they are more likely to lose colonists beyond the domain edges. The random placement of the species range at the start of the simulation means that this will not lead to bias across models.
In the peripatric model, speciation probability is also a peaked function of range size, with the maximum probability of speciation (P = 0.01) occurring when ranges covered approximately 28% of the domain (Fig. 2b). This occurs because although dispersal probability approaches 1 for the largest ranges, an increasing proportion of these events fall beyond the boundary of the domain. We also investigated the case where the probability of dispersal was proportional to the area of the range rather than the perimeter (Supplementary Fig. 1b). Where the probability of speciation is proportional to perimeter length, range shape alters the speciation probability because elongated ranges have a higher perimeter to area ratio (Fig. 2b).
Phylogenetic Tree Shape and Clade Patterns
For each of the 240 combinations of parameter values, we performed 200 simulations, recording the timing of all speciation and extinction events; the coordinates of the species' ranges at the end of the simulation; and the phylogenetic relationships among species (Fig. 1f). From the complete phylogeny including extinctions, we extracted the reconstructed phylogeny by retaining only those branches leading to extant tips. Clade size was recorded as the number of extant species and varied greatly between simulations. We therefore measured the imbalance of the reconstructed trees using the β tree-splitting parameter, which is independent of clade size and, in contrast to other clade imbalance statistics, does not assume a particular model of cladogenesis (Aldous 1996). Code to implement the β tree-splitting parameter was provided by M. Blum (Blum and Francois 2006). Finally, we recorded the range size of the extant species and calculated the proportion of the domain occupied by each species, the skew (g) in the range size frequency distribution for each clade, and the phylogenetic signal in range size using Pagel's lambda (λ) (Pagel 1999): λ varies between 0 and 1, with λ = 1 consistent with a Brownian model and λ = 0 indicating that the trait varies at random with respect to phylogeny (Freckleton et al. 2002).
We measured temporal shifts in diversification rate (ρ) as the proportional difference in rate between the first (r1) and second (r2) half of the reconstructed phylogeny (ρ =(r2 −r1)/(r1 +r2)). The rate within each half, split at the midpoint between the tips and the crown group divergence, was calculated as r =(log(n2) −log(n1))/t (Magallon and Sanderson 2001), where n1 and n2 are the number of extant species at the start and finish of the time period and t is its length. The parameter ρ varies from −1 to 1, with ρ = 0 indicating constant diversification rates, ρ > 0 indicating a speedup in diversification rate, and ρ < 0 indicating a slowdown in diversification rate towards the tips of the tree. We implemented this method in preference to the γ statistic (Pybus and Harvey 2000) because γ exhibits scaling with clade size when changes in diversification rate are not constant (McPeek 2008).
Bird Trees and Birth–Death Models
To examine how our geographical models departed from constant rate models, we compared β and ρ from our simulations with that predicted by the birth–death model, under a variety of extinction rates (b =0.2,d =0,0.05,0.2), run for 20 time steps. Finally, we compared the results of our geographical model with the imbalance and temporal diversification in 41 species-level bird phylogenies taken from Phillimore and Price (2008). These phylogenies, at or around the genus level, typically contained more than 80% of the recognized species from the clade. To examine the effects of incomplete sampling, we randomly removed 20% of the species from our simulated trees and recalculated tree imbalance and the shift in diversification rate. This did not qualitatively alter our conclusions and so we focus on the results from our complete clades (Supplementary Figs. 2 and 3). Following Phillimore and Price (2008), only lineages up to the last bifurcation event prior to 2 Ma were counted when estimating the shift in diversification rate because more recent splitting events may not be recorded as different species (Purvis et al. 2009). We calculated the range size of the species in our bird trees using a global database on avian breeding distributions (Orme et al. 2005) according to the standard avian taxonomy of Sibley and Monroe (1990). We used the combined area of the biogeographic realms (Olson et al. 2001) within which a species was distributed as a measure of domain size, and for each clade calculated the average proportion of the domain occupied by each species. Finally, we calculated the skew in the range size frequency distribution and the phylogenetic signal in range size (Pagel 1999) of each clade.
Temporal Dynamics of Species Richness and Range Size
In our model, the phylogenetic patterns of diversification are determined by an interaction between the dynamics of species ranges and the geographic mode of speciation. Whether through the bisection of geographic ranges (vicariance model) or the budding of small peripheral isolates (peripatric model), speciation leads to a reduction in geographic range size between parent and the average of the daughter lineages. When ranges have either no or little tendency to expand (μ =0 −0.05), this successive apportioning of species distributions results in a rapid decline in range size through time. Clades that by chance experience rapid diversification become characterized by more restricted ranged species, resulting in a strong negative relationship between clade richness and geographic range size (Fig. 3). By the end of the simulations, the median range of the extant lineages typically covers less than 5% of the available area (Fig. 4g–j). Consequently, although speciation probability is a peaked function of range size, the restricted nature of most geographic ranges means that for the majority of species, reductions in range size generally lead to lower probabilities of speciation.
Declining range sizes lead to negative values of ρ, indicating strong temporal slowdowns in the net rate of diversification. This is particularly apparent under the vicariance model (Fig. 5i), where levels of ρ are similar to the strongest declines in diversification rate found in our selection of bird genera (Fig. 5g). In contrast, under the birth–death model only constant or accelerating diversification rates are predicted, depending on the relative extinction rate (Fig. 5h). When ranges have a stronger tendency to expand, this counteracts the trend toward smaller range sizes, increasing the diversification rate and leading to larger clade sizes (Fig. 3). When the rate of expansion balances the decline in range size caused by speciation, temporally homogenous diversification rates occur (Fig. 5i,k). At the highest rates of range expansion that we explored (μ =0.4,0.6), however, geographic ranges tend to fill the domain, reducing the probability of speciation and leading to strong temporal slowdowns in the net rate of diversification (Fig. 5i,k).
When rates of range movement are high (σ2 = 5), species with restricted ranges rapidly become extinct, resulting in smaller clade sizes (Fig. 4b–e) and an increase in average range size (Fig. 4g–j). These high rates of extinction lead to an apparent acceleration in the rate of diversification through time (Fig. 5j,l). A strong tendency for range expansion prevents small ranged species from going extinct, so that diversification returns to a more constant rate through time (Fig. 5j,l). Under both speciation models, values of ρ most consistent with that observed amongst bird genera, occur under low rates of range expansion (μ = 0.05) and little stochasticity in range edge movement (σ2 = 0.05).
Phylogenetic Tree Imbalance and the Distribution of Range Sizes
Under most combinations of speciation mode and range dynamics our model predicts positively skewed range size frequency distributions (RSFDs) (Fig. 4l–o), mirroring the strong skew in range sizes observed across bird genera (median g =1.43) (Fig. 4k). In our vicariance model, this occurs due to the tendency for barriers to result in the asymmetric division of species ranges (Fig. 4l). The peripatry model produces particularly strongly skewed RSFDs due to the persistence of a few species with large geographic ranges and production of many small peripheral isolates (Fig. 4n). The closest approximation to the skew observed in the RSFDs of bird genera occurs when speciation is by vicariance and rates of range expansion are low (μ = 0.05), a result that is insensitive to the degree of stochasticity in range movement. Realistic levels of skewness also occur in the peripatric model, but the conditions are more restrictive, requiring both low rates of range expansion (μ = 0.05) and highly stochastic range movement (σ2 =5).
The differences in range sizes among species lead to variation across lineages in the probability of diversification. Under a peripatric model, with no tendency for ranges to expand (μ =0), highly unbalanced trees are produced (median β = −0.71) (Fig. 5e) compared with that expected under a birth–death model (Fig. 5b). Indeed, levels of imbalance are comparable to that observed in real bird trees (median β = −0.51, Fig. 5a). This is because the production of small peripheral isolates does not reduce the range size of the ancestral lineage, which is able to continue speciating at the same rate. In contrast, despite variation in range size among lineages, trees produced under a vicariance model are highly balanced (median β =1.8) (Fig. 5c,d) compared with observed clades of birds (Fig. 5a) and the birth–death model (Fig. 5b). This occurs because, although range sizes are skewed, lineages that experience rapid diversification also exhibit a reduction in range size, reducing opportunities for further cladogenesis. An additional factor leading to more balanced tree topologies in the vicariance models is the tendency for barriers to simultaneously bisect multiple species ranges.
The dynamics of ranges following speciation has complex effects on tree imbalance and range size distributions. When rates of range movement are low (σ2 = 0.05), but ranges have a tendency to expand, this leads to more balanced trees because even species with initially restricted distributions are likely to speciate (Fig. 5c,e). However, when rates of range movement are high (σ2 =5), the effect of range expansion on tree balance is reversed. This is because when there is no tendency for range expansion, species with small ranges rapidly go extinct leading to less skewed RSFDs and more balanced trees (Fig. 5d,f). In contrast, when ranges have a tendency to expand, species with small ranges are saved from extinction leading to more positively skewed RSFDs (Fig. 4m,o) and more imbalanced trees (Fig. 5d,f), converging on the expectation of the equal rates birth–death model (median = 0.19) (Fig. 5b).
Phylogenetic Signal in Range Size
Within bird genera, range size tends to exhibit a low phylogenetic signal (λ <0.2 in 29 of 41 cases) (Fig. 4p). Similarly, low levels of λ are obtained under both the vicariance and the peripatric model (Fig. 4q–t). A weak phylogenetic signal is found regardless of the rate at which ranges fluctuate through time, occurring under both highly volatile and highly stable conditions (Fig. 4q–t). The main cause of low λ estimates is therefore not the rate at which ranges are stochastically changing in size through time (simulations show that the product of two Brownian processes on its own leads to high estimates of λ, similar to a Brownian process). Rather, the weak phylogenetic signal arises because of the highly non-Brownian way ranges are “inherited” during speciation under the vicariance and peripatry models. High rates of range expansion lead to higher values of λ, by increasing the variance of range sizes among the ancestors of present-day species (Fig. 4q–r). Eventually, however, under very high rates of expansion, species tend to fill the domain and the phylogenetic signal subsequently declines (Fig. 4q–r).
The Geography of Diversification
A fundamental assumption of the birth–death model is that the probabilities of speciation and extinction are independent of the history of cladogenesis (Losos and Adler 1995). However, the bifurcation of a lineage in a phylogenetic tree will usually arise because of the splitting of its geographic range in space (Coyne and Orr 2004). Newly formed sister lineages will, therefore, occupy only a fraction of the parental range and may differ markedly in size depending on the geometry of range division and mode of geographic isolation (Anderson and Evensen 1978; Barraclough and Vogler 2000; Waldron 2007). Consequently, if the probabilities of speciation and extinction rates are dependent on the geographical extent of a species (Rosenzweig 1975), variation in the rates of diversification through time and among lineages are expected (Losos and Adler 1995). Here we show that considering the geography of speciation and extinction has a profound influence on the phylogenetic patterns of diversification expected in the absence of ecology.
Highly unbalanced trees and strong temporal slowdowns in diversification are often observed in real phylogenies but do not arise under the standard birth–death model, leading to the widespread perception that species biological traits and ecological interactions underlie the dynamics of species radiations (Slowinski and Guyer 1993; Mooers and Heard 1997; Weir 2006; Phillimore and Price 2008). However, our models show that these patterns can also arise simply because of the feedback on speciation and extinction rates imposed by the division of geographic ranges during speciation. These results suggest that even for groups representing classic cases of adaptive radiations (Losos et al. 1998; Schluter 2009), we should not discount the possibility that the rate of diversification has been limited by the geographical opportunity for isolate formation and range expansion rather than the availability of unutilized niche space (Rundell and Price 2009) (e.g., many of the Greater Antillean Anolis sister species are in fact allospecies found on separate mountains on a single island; Losos and Parent 2009). This neutral perspective does not preclude the role of niche expansion and adaptation in initiating species radiations (Nee et al. 1992). Indeed, if range sizes decline due to speciation, then rare events that promote range expansion, such as key innovations, may be essential in renewing the process of diversification.
The Geographical Mode of Speciation
The way in which geographical isolation occurs in our model has a strong effect on the dynamics of cladogenesis. Under a vicariance scenario, speciation leads to successive reductions in range size and a declining rate of diversification, particularly in those branches where speciation has been most rampant. The result is a strong temporal decline in the diversification rate and highly balanced topologies. When speciation occurs via dispersal, the production of small peripheral isolates also leads to a decline in the average range size and thus slowdown in speciation rate. However, in contrast to the vicariance model, the range size of one of the daughter lineages is not eroded by the process of speciation. Lineages that have been prolific speciators will therefore continue to be so, producing extremely unbalanced tree topologies. A similar result is obtained for Hubbell's point mutation speciation model because here most newly formed species are extremely rare, that is, singletons (Mooers et al. 2007).
Variation in range shape is also an important determinant of speciation in our model. For a given range size, linear distributions are more prone to bisection by geographic barriers (Rosenzweig 1978; Graves 1988) and more likely to give rise to dispersal events. Mathematical models of the genetic basis of speciation also predict that one-dimensional geographic distributions should lead to higher speciation rates because of lower connectivity and gene flow between populations (Gavrilets 2004). The dynamic relationship between range shape and diversification rate in our models is similar to that of range size: in the vicariance model, the random division of a linear range tends to give rise to daughters with more compact distributions, further confounding the potential for future speciation, whereas no such trend occurs in the dispersal model.
Range Evolution and the Role of Geographical Speciation
For continental radiations, it has been argued that the most likely way in which niche saturation impinges on diversification rates is through a heightening of competition and the progressive limits to range expansion and thus opportunities for speciation (Rosenzweig 1975; Price 2008; Phillimore and Price 2009). Such a model implies that in the absence of competition, ranges would rapidly expand following speciation. Indeed, when ranges have a strong tendency to expand in our models, rates of diversification are more constant through time. However, species distributions, particularly at large spatial scales (Russell et al. 2006), are often thought to be set for reasons other than competition, including physiological limits (Root 1988), habitat structure (Holt et al. 2005), and dispersal barriers (Goldberg and Lande 2007). These limits may be stable over evolutionary time due to a lack of genetic variability (Kellermann et al. 2009), developmental constraints, or genetic swamping preventing adaptation (Kirkpatrick and Barton 1997). Although species may be able to adapt and subsequently expand their distributions, the timescales for niche evolution may be comparable to that typically assumed for speciation (Wiens and Donoghue 2004). A signal of slowdowns in the diversification rates of large clades with small geographic ranges is therefore, by itself, insufficient evidence for competition in limiting range expansion and speciation. One way in which the role of competition in generating slowdowns could be tested would be to modify our models to prevent sympatry among lineages that have recently split/arisen.
Previous studies have shown that when range dynamics are highly volatile, the signal of speciation in the distributions of extant species is rapidly eroded, making the geographical mode of speciation difficult to detect (Barraclough and Vogler 2000; Losos and Glor 2003; Phillimore et al. 2008). Here we show that when range volatility is high, the reconstructed patterns of diversification are also drastically altered. Only under conditions of relatively stable range sizes should we expect strong slowdowns in diversification rate and high tree imbalance. Changes in species distributions can be highly volatile, occurring over short timescales compared with that required for speciation (Davis and Shaw 2001; Lessa et al. 2003). However, the extent to which such high rates of range evolution are typical of clades in general is unclear. Phylogenetic signal in range size is typically weak, with most variation occurring among low taxonomic levels (Waldron 2007), and this is often interpreted as evidence that range size is highly labile (Gaston 1998; Freckleton and Jetz 2009). However, our models show that even when range sizes vary little through time according to a random-walk process that on its own leads to range size having a strong phylogenetic signal, the apportioning of geographic ranges during speciation, and the corresponding substantial departure from a Brownian model of evolution, ensures that the phylogenetic signal is low. A faster rate of range evolution on its own leads to similarly high estimates of phylogenetic signal, as is the case with a Brownian model (Revell et al. 2008). However, as faster range evolution leads to a greater frequency of extinctions, this censoring may dampen phylogenetic signal in addition to the effect of the speciation process. Taken together, a perhaps surprising finding of our study is that low phylogenetic signal does not therefore provide strong evidence for the lability of geographic ranges.
Incorporating Geography into Models of Diversification
In our model, altering the range dynamics and mode of speciation gives rise to wide variation in the temporal dynamics and shape of phylogenetic trees. As a result, the appropriate null expectation of imbalance and slowdowns for a clade will differ depending on the conditions characterizing its evolution. For instance, if vicariance has been the predominant mode of speciation within a clade, then a highly unbalanced topology would be even less explicable by a random process than comparison with a birth–death model would suggest. Conversely, if most speciation has occurred via peripatry, then high imbalance would be expected by chance alone and would not require explanations based on ecological differences among species. Moreover, our models show that a pattern of equal and constant rates of diversification is only expected under particular scenarios of range evolution and speciation mode; namely, when rates of range expansion and dynamism are high and speciation occurs by peripatry. When these conditions do not accurately characterize the history of a clade, the use of an equal and constant rates birth–death model would therefore be inappropriate.
A number of studies have examined the effects of altering the assumptions underlying the birth–death model on patterns of phylogenetic tree imbalance. For example, Chan and Moore (1999) showed that when speciation results in a decline in the probability of diversification for both daughter lineages, as in our vicariance model, this leads to more balanced trees. When only one of the daughter lineages exhibits a reduction in speciation rate, as in our dispersal model, more unbalanced trees are produced (Rogers 1996; Chan and Moore 1999). The causes of such refractory periods have generally been cast in a geographic context, for example, the time required for ranges to expand following speciation (Losos and Adler 1995). Heard (1996) and Rogers (1996) found that when per-lineage speciation rates evolved according to either a gradual or a punctuated model, the trees produced were more unbalanced; a scenario approximating the peripatric model of speciation that we investigated. By simulating the underlying dynamics of species ranges and the process of geographic isolation our model, generalizes and extends these findings. For instance, although the impact on speciation of a short refractory period is equivalent to a scenario in which ranges have a strong tendency to expand (Losos and Adler 1995), range expansion also reduces the extinction rate with complex effects on the distribution of geographic range sizes and subsequent dynamics of diversification.
Nevertheless, our model makes a number of simplifications which future studies should address. Vicariance was simulated using a simple model of linear barrier formation, but other ways of modeling range division may lead to differences in how the probability of speciation scales with range size and in the resulting asymmetry in daughter species' range sizes (Rosenzweig 1978; Waldron 2007). In the absence of empirical information on the geographic barriers found in nature, we chose arbitrary barrier length distributions, with the barriers dropped at random onto the domain, whereas in reality barriers may be more frequent at continental margins (e.g., mountain ranges and rivers). A further simplification in our model was the assumption that species range margins shift independently through time. In particular, biogeographic settings, directions of range expansion, and thus patterns of geographic range size and shape may be more constrained (Pigot et al. 2010). Moreover, if changes in the position of range boundaries are correlated, then large shifts in species distributions may occur without substantial changes in geographic range size (e.g., parallel shifts in the northern and southern boundary of a range in response to climate change). Such a scenario may extend the breadth of parameter space under which a geographical model can account for observed patterns of diversification.
The imbalance and temporal dynamics of the bird phylogenies we examined seem to be best explained by our dispersal model with little post-speciation changes in geographic range size. This model, however, predicts unrealistically high levels of skew in the distribution of range sizes, suggesting that the patterns in clade shape across bird genera are unlikely to be explained by a single region of the parameter space we explored. This is perhaps not unexpected given the disparate histories of real clades, including their ancestral range sizes, geographical locations, and opportunities for range expansion (Ricklefs 2006; Weir et al. 2009). Moreover, we did not explore the more realistic scenario where vicariance and dispersal each contributes to the origin of new species. In the future, it may be possible to estimate the combination of geographic parameters that are most likely to have given rise to observed phylogenetic and distribution patterns in a likelihood framework or using approximate Bayesian techniques. For instance, clades diversifying predominantly by peripatry may be characterized by more asymmetric ranges and shorter lived species that those produced by vicariance.
Although a model of equal rates amongst species and constant rates through time is regarded as the most parsimonious model of cladogenesis (Nee 2006), from a geographic perspective, such a scenario requires that the range sizes of daughter species are similar to both that of the parent and to each other, implying either that range size is strongly heritable (Jablonski 1987) or that the geographical signal of lineage splitting is completely scrambled by post-speciation range evolution. Accordingly, in our simulations, the closest approximation to a birth–death model is reached under a peripatric model of speciation with high rates of range growth and highly volatile fluctuations in range boundaries. When viewed in a geographical context, then the birth–death model represents only a small region of parameter space, with the occurrence of equal and constant rates requiring a rather restrictive set of biological conditions to be fulfilled.
This work was funded by a Grantham Studentship to A.L.P. from the Grantham Institute for Climate Change, Imperial College London. A.B.P. is funded by an Imperial College Junior Research Fellowship. C.D.L.O. is funded by an RCUK Fellowship.
We are grateful to Tim Barraclough, Rob Ewers, Trevor Price, Robert Ricklefs, James Rosindell, Gavin Thomas, and two anonymous referees for their helpful comments on earlier drafts of the manuscript and to Yael Kisel, Lynsey Mcinnes, Andy Purvis, and Sarah Whitmee for discussion on these topics. Many of the ideas developed in this paper benefited from discussions at an NERC Centre for Population Biology workshop on Outstanding questions in speciation research.