A Conceptual and Statistical Framework for Adaptive Radiations with a Key Role for Diversity Dependence

Introduction

Adaptive radiations are expected to leave two distinct signatures on diversification (defined as origination/speciation minus extinction) patterns (Yoder et al. 2010): a large initial increase in diversification rate and a subsequent decrease in the diversification rate. The initial, sudden positive rate shift is generally thought to be triggered by ecological opportunity (Simpson 1949, 1953; Schluter 2000; Yoder et al. 2010). There are three types of ecological opportunity (Simpson 1949, 1953) through which species can enter the “adaptive zone”: (1) extinction of antagonists; (2) appearance of a new environment, either due to dispersal to a new area or due to external forces changing the environment; or (3) a key innovation, that is, the appearance of aspects of organismal phenotype that promote diversification (Hunter 1998) by enabling escape from competition for niche space, by reducing the probability of extinction by increasing population density via increased individual fitness or by favoring reproductive or ecological specialization (Heard and Hauser 1995). All types of ecological opportunity result in ecological release. Types 1 and 2 require (ecological) changes in the environment of the focal species, whereas type 3 requires an evolutionary change in the focal species itself. However, the three types of ecological opportunity are not independent: novel species traits might promote positive net diversification only in a new biotic (e.g., after antagonist extinction) or abiotic environment or in the presence of other traits (De Queiroz 2002). For this reason we believe that given a clade, it is more practical to consider two types of adaptive radiations: an increased diversification of all lineages of a clade and increased diversification of a single (or a few) lineages. We immediately admit that this dichotomy is a matter of scale, because a clade where all lineages are affected can be part of a larger clade where it is the only affected subclade, but the distinction is useful for a given clade.

While the initial, sudden positive shift in diversification rate is obviously evidence for adaptive radiations, the subsequent more gradual negative rate shift is just as indicative of adaptive radiation (Albertson et al. 1999; Lovette and Bermingham 1999; McPeek 2008; Phillimore and Price 2008). Because the negative rate shift is observed across clades of different ages and for different taxa, it seems unlikely to be caused solely by direct time dependence of diversification rates through external forcing (e.g., climate change reducing or increasing extinction rates), but it is the accumulation of species within the clade itself that limits further diversification. Therefore, this phenomenon is called diversity-dependent (or density-dependent) diversification (Valentine 1973; Sepkoski 1978, 1979, 1984; Walker and Valentine 1984; Schluter 2000; Rabosky 2009). The negative shift in diversification rate can, in principle, be due to a reduction in speciation rate or an increase in extinction rate as diversity increases (Sepkoski 1978), but an increase in extinction rate has been shown to be a very unlikely form of diversity dependence (Alroy 1996; Rabosky and Lovette 2008b; Burbrink and Pyron 2009; Quental and Marshall 2009). The reduction in speciation rate is considered to be a consequence of accessible niches becoming filled (Simpson 1953; Valentine 1980; Walker and Valentine 1984; Schluter 2000; Rabosky and Lovette 2008a; Gavrilets and Losos 2009) or a consequence of decreased range sizes following successive subdivisions of ancestral ranges (Sepkoski 1978; Rosenzweig 1996). Both mechanisms increasingly limit opportunity for ecological speciation and suggest that there is a clade-level carrying capacity for the number of species than can exist at any one time.

In this article we will present a framework for adaptive radiation that respects both the positive and negative rate shifts in a single coherent framework. In this framework, the negative shifts are due to diversity dependence and positive shifts are either due to (1) changes in the speciation and extinction rates or clade-wide shifts in carrying capacity or (2) decoupling of the diversity-dependent dynamics of a subclade from the main clade’s diversity-dependent dynamics, thereby affecting the subclade’s diversification rates. The decoupling occurs because of a key innovation or dispersal event that enables escape from diversity dependence regulated by competition for niche space. We emphasize that the key innovation could have occurred much earlier than the decoupling if the innovative trait required a new environment to cause the decoupling. We will compare the subclade-specific model of decoupling of diversity dependence with a model that assumes a clade-wide shift in speciation rate, extinction rate, or carrying capacity.

In the past decade molecular phylogenies have come to the fore as a valuable additional source of information on diversification patterns, especially where fossil evidence is lacking. Both temporal and topological methods have been developed to detect significant shifts in diversification rates (reviewed in Sanderson and Donoghue 1996; Mooers and Heard 1997; Barraclough and Nee 2001; Chan and Moore 2002; Ricklefs 2007). Temporal methods compare the observed distribution of speciation events through time, inferred from the lineage accumulation through time (the lineages-through-time plot; Paradis 1997), with the distribution expected by a null model of cladogenesis. These have the disadvantage that they identify only the timing of the change in the rate of lineage accumulation but not its location in the phylogenetic tree. Earlier, Farris (1976) had already argued that asymmetric phylogenies are more common than symmetric ones, and thus this asymmetry (or imbalance) conveys information. Topological methods build on this fact and focus on the relative diversity in the descendant subclades of a node (Vrba 1980). Single-node topological tests compare the observed difference in sister-group sizes to the expectation under a null model of cladogenesis (Slowinski and Guyer 1989a, 1989b, 1993; Goudet 1999). However, because these methods compare only present diversity to diversity at the clade’s origin, they cannot detect diversification rate increases and subsequent decreases (Hunter 1998). Whole-tree methods compute a single statistic of tree imbalance or asymmetry from the relative diversity of nodes throughout the entire tree (Shao and Sokal 1990; Kirkpatrick and Slatkin 1993; Chan and Moore 2002). These nonparametric approaches are robust in detecting rate shifts but do not lead us toward the underlying causes of these rate shifts. Alfaro et al. (2009) presented an algorithm to find significant rate shifts based on the likelihood of the phylogeny, but they used the constant-rate birth-death model and thus ignored the second phase of an adaptive radiation, despite the ample evidence for this second phase. This implies that this method may miss a shift in diversification due to a change in the carrying capacity or due to decoupling of diversity-dependent dynamics. Here we present not only the model of adaptive radiation incorporating its two phases in the light of diversity dependence as outlined above but also a likelihood-based inference method to detect decoupling of diversity dependence given the whole phylogenetic tree. This method accounts for asymmetry as well as the distribution of branching times and can handle incomplete phylogenies (if one knows the number of species that are not placed in the phylogeny but do have the same ancestral crown group).

We start by describing our general framework for studying adaptive radiations. We then describe our inference method for detecting decoupling of diversity-dependent dynamics, and we apply it to the radiation of the North American clade of the legume tribe Psoraleeae studied by Egan and Crandall (2008a). They found evidence for diversification rate shifts in this tribe, using both temporal and topological methods. We confirm this finding using our single unified approach. Moreover, not only do we estimate the diversification parameters, we also identify the timing of the decoupling.

A Coherent Framework for Adaptive Radiations

We view diversification as a birth-death process (Kendall 1948; Bailey 1964). This is the standard mathematical model of macroevolutionary diversification (Nee 2006) that has been frequently used in the study of adaptive radiations (e.g., Alfaro et al. 2009; Steeman et al. 2009). However, rather than assuming time-constant speciation and extinction rates or direct time dependence of these rates as has been done previously, we allow these rates to be diversity dependent. In general, it makes the most sense, as argued in the introduction, to assume that the diversity dependence is entirely in the speciation rate, but we note that our framework is easily amenable to account for diversity-dependent extinction. In the example we will study, we assume that the per-species speciation rate declines linearly with diversity (Etienne et al. 2012). The diversity at which the speciation rate is equal to the extinction rate is the carrying capacity of the clade, analogous to the definition of carrying capacity in population dynamics models. This diversity dependence describes phenomenologically the filling of accessible niches and hence the second phase of an adaptive radiation. Note that, at equilibrium, when speciation and extinction rates equal one another, not all niches are filled, unless the extinction rate equals 0.

In the first phase of an adaptive radiation, there are shifts in diversification parameters. These shifts can be clade-wide or subclade-specific (Table 1). We argue that diversity dependence again plays a key role, in both types of shift. The clade-wide rate shifts occur when external factors create new ecological opportunities for all lineages (McInnes et al. 2011), for instance, climate change (Ezard et al. 2011) or tectonic movement (e.g., cetaceans diversifying more rapidly due to restructuring of the oceans; Steeman et al. 2009). These external factors, in creating more niche space, most likely affect the clade-level carrying capacity rather than the intrinsic speciation rate or extinction rate.

Subclade-specific shifts result from either shifts in intrinsic speciation rate and extinction rate for the subclade, or (partial) decoupling of diversity-dependent dynamics of the subclade from the main clade. Shifts in intrinsic speciation rate and extinction rate crucially interact with diversity dependence: because all species compete for the same limited number of niches, those with elevated speciation rates or reduced extinction rates will eventually exclude the other species. This may explain evolutionary succession, where radiation of one taxon is accompanied by decline in another (Valentine 1973), for example, the replacement of therapsids by archosauromorphs during the Triassic facilitated by the latter’s higher growth rates (Sookias et al. 2012) or the replacement of Diplograptina by Neograptina during the late Ordovician (Bapst et al. 2012). Decoupling of diversity-dependent dynamics, complete or partial, can occur after key innovations or dispersal to a new environment, where the focal species’ diversification is released from diversity dependence. The focal species starts its own subclade that does not “feel” the presence of the other species, simply because the key innovation or dispersal event has allowed it to escape competition for niche space as it has moved far away in trait space or in geographical space, respectively. Note that dispersal to a new environment can be interpreted as a key innovation (De Queiroz 1998). We adopt this broad interpretation of key innovation from hereon (if only for notational convenience). In this article, we concentrate on the decoupling of diversity-dependent dynamics, as we believe it to be a major process in causing adaptive radiations and it increases our understanding of the role of key innovations. For instance, it is likely that evolutionary succession also involves (partial) decoupling triggered by innovative traits.

Extinction of antagonists can be positioned at several places in our framework. When an antagonist becomes extinct, the surviving species do(es) not experience the presence of the antagonist, simply because it is no longer there, but this does not necessarily constitute a decoupling of the dynamics. If it is merely a random extinction opening up the niche of the antagonist, then the extinction event is already part of the usual speciation-extinction dynamics. If the antagonist extinction results in a lower extinction rate of the surviving species, it can be seen as time dependence of the extinction rate and should thus be viewed as an example of a (temporary) rate shift. If the extinction of the antagonist creates more niches than just leaving its own, it can be viewed as a shift in carrying capacity. Only if these new niches are not available to all species, it constitutes decoupling of diversity-dependent dynamics.

The idea of decoupling of diversity-dependent dynamics can apply to a single species or to a group of species. In the case of a group of species where the decoupling is caused by an innovative trait, it is not necessary that this trait evolves multiple times; the trait may evolve in a single species and stay dormant until a change in (a)biotic environment activates it, by which time the species may have diversified. In our mathematical implementation of the decoupling model, we assume for simplicity that only a single species undergoes the decoupling event. Furthermore, to keep our focus, we do not look at shifts in subclade-specific intrinsic speciation and extinction rates without enabling escape from diversity dependence, that is, without decoupling of diversity-dependent dynamics. Nor do we consider partial decoupling.

Our model of key innovations causing decoupling of diversity-dependent dynamics is depicted in figure 1. The decoupling occurs somewhere along the lineage that will form a new subclade. We denote the time of this occurrence by . Time is a parameter in our model, so it can be estimated together with the other parameters of the model rather than set beforehand at fixed points (Alfaro et al. 2009) or allowed only at branching points (Rabosky 2006). Up to the lineage and its descendants contribute to the diversity dependence of the main clade, which we will call M. After the lineage no longer has an effect on the diversification in the main clade, but it now affects only the diversification in the subclade it initiated, which we will call S. We thus have a model with seven parameters: the decoupling time , three parameters for the diversification of the main clade (speciation rate, extinction rate, and carrying capacity), and similarly, three for the subclade. The number of parameters may be reduced by assuming that the intrinsic speciation rate and/or the extinction rate are the same for main clade and subclade. A model comparison based on an information criterion (we will use the Akaike information criterion [AIC]) will single out the most parsimonious model for the data.

Likelihood of a Phylogeny under the Key Innovation Model

The birth-death process (Kendall 1948; Bailey 1964) can be mathematically described with the master equation for the probability of n species at time t which we denote by :

2222

2222

We start with the algorithm to compute the likelihood of a phylogeny for the diversity-dependent diversification process (Etienne et al. 2012). It is based on the repeated integration of a system of ordinary differential equations. Consider a time t between two branching times where the phylogeny has k branches. Denote by the probability that a realization of the diversification process is consistent with the phylogeny up to time t and has n species at time t. The dynamical equation is (Etienne et al. 2012):

2222

At the branching times, we impose a branching event by multiplying the probabilities by the speciation rate , and increasing by 1 the number of branches k. This leads to the algorithm in box 1 to compute the probability of a phylogeny under the diversity-dependent birth-death model (Etienne et al. 2012). The resulting likelihood differs by a combinatorial factor , with q being the number of species, from the likelihood computed by Etienne et al. (2012). The latter allows direct comparison with the likelihood of Nee et al. (1994), but the combinatorial factor should be dropped if one wants to compute the likelihood of a phylogeny rather than the likelihood of a set of branching times.

Here we extend the algorithm of box 1 to account for the key innovation (KI) causing decoupling of dynamics at time . We separate the likelihood computation in a part for the main clade M and a part for the subclade S. For the main clade M, we impose the key innovation event at time by decreasing by one the number of branches k. For the subclade S, we use the algorithm in box 1, starting from the key innovation time with . We further define the branching times in the main clade before the key innovation, , and after the key innovation, , and in the subclade, . See also figure 1. Finally, if there are missing species in the phylogeny, they can be specified as and . Because we often do not know whether the missing species are in the main clade or the subclade, we can sum over all possible combinations of and such that , where m is the total number of missing species. Box 2 shows the algorithm to compute the likelihood of a phylogeny with main clade M and subclade S. This algorithm gives the likelihood without requiring that the two crown lineages survive. To condition on their survival, as advocated since Nee et al. (1994), we need to follow the number of descendants of the two crown species from time to time and compute the probability that the number of descendants of both are positive. We refer to appendix A for the technical details of the conditioning and to appendix B for details on how missing species are handled. We have made the likelihood computation available in the R package DDD (http://cran.r-project.org/web/packages/DDD/index.html), from version 1.0 onward. We also have Matlab code available on request.¹

Case Study

We illustrate our method by applying it to the North American clade of the legume tribe Psoraleeae (Egan and Crandall 2008a); see figure 2. The phylogenetic data for this clade are available in the Dryad Data Repository (Egan and Crandall 2008b). This clade contains 47 species (43 are in the phylogeny, 4 species are missing) within 5 genera. The most species-rich genus is Pediomelum (29 species) and contains 3 subgenera, Leucocraspedon, Disarticulatum, and Pediomelum, the first of which contains only 2 species. These asymmetries suggest an adaptive radiation in the latter two subgenera, or more precisely—assuming that antagonist extinction did not occur—a key innovation. Egan and Crandall (2008a) found evidence for significant rate shifts, using both temporal and topological approaches, but did not find support for the hypothesis that adaptation to xeric habitats (hot and dry environments) acted as key innovation. Egan and Crandall (2008a) dated the origin of this North American clade at 5.8 million years ago (mya), and of the genus Pediomelum at 3.2 mya. The major radiation of this genus did not start until around 2 mya, when the subclade containing subgenera Pediomelum and Disarticulatum diverged from the main clade (including subgenus Leucocraspedon). Egan and Crandall (2008a) detected a major shift in diversification rate around this time.

We tested whether our approach can indeed identify the seemingly evident key innovation around 2 mya or that the asymmetry occurs by chance and the increased branching is due only to a shift in diversification rates affecting the entire tree. We did this for two possible scenarios for the key innovation (see fig. 2). We determined the timing of this key innovation and whether there was only decoupling of diversity-dependent dynamics with main clade and subclade having their own carrying capacity or whether the intrinsic speciation rate and extinction rate also differed between the two. This resulted in a variety of models. The constant-rate (CR) models we tested were the standard (diversity-independent) birth-death model (CR0) and the diversity-dependent DDL+E model (CR1). The shifting rate (SR) models we tested include the Yule 2 rate model (SR0) that was also employed by Egan and Crandall (2008a) and models with a shift in the carrying capacity K (SR1), in K and the extinction rate μ (SR2), and in K and the intrinsic speciation rate (SR3). Our key innovation (KI) models include models where the subclade differs from the main clade only in K (KI1), in K and μ (KI2), in K and (KI3), and in K, , and μ (KI4). The models with a 0 denote models that we consider unrealistic because they either ignore extinction (SR0) or diversity dependence (CR0, SR0), but we analyzed them anyway as they are standard null models.

The key innovation models clearly stand out as the best description of the data regardless of the scenario (fig. 2) of the key innovation and regardless of whether the data set was truncated to compensate for lack of species recognition, as tables 2 and 3 show: the Akaike weights of the KI models were all greater than , while for the CR and SR models, they were all smaller than ). The scenario of figure 2B, where the key innovation occurs between the branching points at 2.02 and 1.73 mya seems more likely than the scenario of figure 2A, where it occurs between the branching points at 3.21 and 2.02 mya. Among the key innovation models, most support is lent to the model where the subclade not only differs in the carrying capacity K but also in the intrinsic speciation rate (KI3). This suggests that not only did the decoupled subclade experience a range of new niches, it was also able to colonize them more rapidly than would be expected from a reduction in the influence of diversity dependence alone. However, the Akaike weights of the KI models were similar and thus do not justify a strong preference for any particular model.

Egan and Crandall (2008a) suggested that the final negative rate shift may be due to lack of species recognition. This has been termed protracted speciation (Rosindell et al. 2010; Etienne and Rosindell 2012): incipient species go unrecognized because they have not completed speciation. While the most elegant way to account for this is to explicitly model protracted speciation (Etienne and Rosindell 2012), an expression or algorithm for the likelihood of this model given the phylogeny is not available when diversity dependence is acting. Hence, we followed Phillimore and Price (2008) by chopping off the last 0.22 million years (during which no branching occurs), and we repeated the analysis on this truncated phylogeny. Even though this is a rather crude method, we believe the obtained parameter estimates to be more realistic than the ones obtained without accounting for protracted speciation. For instance, applying our method to this truncated data set, that is, the data set that compensates for unrecognized species clearly yields higher extinction rates than when applied to the full data set where the extinction estimates are often close to 0 (Table 1).

We have to base our model comparison on the parameter estimates, because we cannot use AIC to compare the full and truncated data sets, as such a comparison is allowed only for different models on the same data set. To make a proper comparison of models on the same data set, one could interpret the analysis of the truncated data as applying to the full data set a model that assumes no speciation or extinction events occur in the last 0.22 million years. This model then performs better but has quite a complex structure, because it assumes one set of free parameters before 0.22 mya and another set of fixed parameters after 0.22 mya. The AIC is not well suited for penalizing this complexity for lack of parsimony.

Discussion

Central to the study of adaptive radiations is distinguishing chance variation from that which requires a more deterministic explanation (Chan and Moore 2002). Equal-rates models such as the constant-rate (diversity independent and diversity dependent) and shifting-rate birth-death models (i.e., models where the speciation and extinction rates are independent of the lineage at a given time; Rabosky 2006) assign identical probability to all topologies, including imbalanced ones. Observing a single imbalanced tree therefore does not constitute evidence against such a model (Slowinski and Guyer 1989a; McConway and Sims 2004). Still, a model that assigns a higher likelihood to imbalanced topologies will likely be preferred as indicated by likelihood ratio tests or information criteria after accounting for the number of parameters (McConway and Sims 2004). In this article, we have presented a model that does assign a higher likelihood to imbalanced trees at the cost of at least one but possibly four parameters. In our model, imbalance is caused by a key innovation (which we have broadly defined to include innovative traits as well as dispersal to a new environment) that decouples the diversity-dependent dynamics of a subclade from those of the main clade.

While we focused on decoupling of diversity-dependent dynamics in this article, we have presented a more general framework of adaptive radiations (Table 1) where diversity dependence plays a key role. First, diversity-dependent diversification explains the slowdown of diversification after an initial burst. Second, the initial burst can best be understood in the light of diversity dependence. When the burst is clade-wide, then a likely cause is a shift in carrying capacity, hence a temporary alleviation of the pressure of competing species. When the burst is subclade-specific, then without decoupling of diversity-dependent dynamics, the limit to diversity will crucially determine the outcome of competition for niches: species with elevated speciation rates or reduced extinction rates will be expected to dominate. With decoupling of diversity-dependent dynamics, a subset of species will be temporarily alleviated from competitive pressure. Thus, diversity dependence is a pervasive factor in adaptive radiation.

We have illustrated our method for detecting decoupling of diversity-dependent dynamics assuming that only a single key innovation occurs (i.e., only one subclade is treated differently). The method can be readily extended to detect more than one key innovation, similar to the way Alfaro et al. (2009) detect the locations of shifts in diversification rates, using the standard diversity-independent birth-death model (see also Sanderson and Donoghue 1996). However, the power of any method to find multiple key innovation events will be limited, unless the phylogeny is large, because each new key innovation event involves at least one (the key innovation time), most likely two (the key innovation time and a different K value for the subclade, e.g., KI1), and perhaps even four (e.g., KI4) additional parameters. Also, our results suggest that the data are not informative enough (e.g., the phylogeny is not large enough) to allow distinguishing clearly between the various key innovation models (KI1-KI4). Because the parameter estimates may vary substantially between these models, one should also be careful in interpreting the actual estimated parameter values (this can be more rigorously checked by computing confidence intervals or posterior probability distributions, which are both possible in our likelihood framework). Similarly, we expect that it will be difficult to distinguish between different models of diversity dependence (e.g., a power-law vs. linear decline of speciation rate with diversity) and between decoupling of diversity-dependent dynamics and subclade-specific shifts in speciation and extinction rates without decoupling. For example, successive radiations where one group of species replaces another can be explained by decoupling of diversity dependence due to antagonist extinction or by higher diversification rates in the new group due to a key innovation allowing them to overtake the old one without decoupling. An additional, more technical, caveat is that with many parameters, one may easily get stuck in a local likelihood optimum. Although we used various starting points, we do not rule out that one or more of the likelihood optima are not the global optima. Tables 2 and 3 show such local likelihood optima for the shifting-rate model.

Egan and Crandall (2008a) found no support for the hypothesis that a shift from mesic to xeric habitat between 3.21 and 2.02 mya acted as a key innovation. Our method does reveal a key innovation, but we find more support for a somewhat later key innovation, between 2.02 and 1.73 mya. The key innovation is therefore probably not associated with this shift from mesic to xeric conditions. Note that we have looked only at the phylogeny rather than at traits that are associated with the subclade that underwent the key innovation (Heard and Hauser 1995). Incorporating character evolution (see, e.g., Ree 2005; Maddison et al. 2007) into our model or, conversely, incorporating (decoupling of) diversity dependence in models of character evolution, would be an obvious, but mathematically challenging, next step (see Mahler et al. 2010 for an interesting attempt).

Rather than finding evidence for a key innovation due to a shift from mesic to xeric conditions, Egan and Crandall (2008a) detected an overall shift in diversification rates and attributed this to Quaternary climate changes, spurring the diversification of the subclade containing subgenera Disarticulatum and Pediomelum. They suggest that the glacial cycles coupled with the varied topography of the southwestern United States created many unique habitats that were ripe for the colonization of the ancestors of these two subgenera. Our analysis supports this suggestion.

It is still possible that a (climate-induced) clade-wide rate shift also occurred during the same time. Indeed, the phylogeny suggests that around 2 mya there is increased diversification in all lineages. Our method can be readily extended to detect clade-wide rate shifts and subclade-specific key innovations simultaneously.

Our model of adaptive radiation through key innovations does not require the assumption that speciation is intrinsically more rapid during adaptive radiation as proposed by Schluter (2000). Glor (2010) scrutinizes this assumption, because it lacks an established mechanistic explanation. Our model provides such a mechanistic explanation. Under most realistic parameter combinations, decoupling of diversity dependence automatically leads to exceptional diversification of the subclade. But if, for example, the carrying capacity of the subclade is substantially lower than that of the main clade, the intrinsic speciation rate and the extinction rate are similar to those of the main clade, and the main clade is still far from equilibrium, then we might not expect a strong asymmetry to arise due to the key innovation, and we would fail to detect it. Yet, a low carrying capacity of the subclade means that there is not much ecological space for the key innovation, so one may wonder whether this would still be considered a key innovation. The justification sought by Glor (2010) for extraordinary diversification to be a necessary ingredient for adaptive radiation comes from the implicit assumption that innovations should create opportunities. If they do not, one can perhaps hardly speak of ecological release. Thus, our model offers a plausible explanation of the problem identified by Losos (2010) that clades may fail to radiate although seemingly in the presence of ecological opportunity.

The property that speciation takes time, that is, that it is “protracted” (Rosindell et al. 2010), can in principle be incorporated in our model along the lines of Etienne and Rosindell (2012). As we explained above, we simply chopped of the last 0.22 million years in our analysis to mimic protracted speciation (Phillimore and Price 2008), because deriving a full likelihood of a phylogeny under both diversity dependence and protracted speciation remains a challenge. The analysis yielded more realistic (nonzero) extinction rates, suggesting that protracted speciation is an important factor to consider in future models of diversification.

There is evidence that, overall, diversity increases over time (see Wiens 2011 for some examples). This seems to conflict with evidence for ecological limits to diversification. Our model resolves this apparent conflict as follows: diversity increases by key innovations, but when key innovations are absent, diversification rates decrease due to diversity dependence (i.e., ecological limits). The end result is still a diversity increase but an irregular one. Of course, this picture is a bit of a caricature. In reality, innovations occur all the time, but some are (much) more influential than others and can be detected in phylogenies because they create many new opportunities causing a rapid radiation. Thus, our model reconciles the contradictory findings on the correlation between clade size and clade age. Clade size and clade age are expected to be correlated for clades where major key innovations are rampant. On average, older clades will have more such key innovations and thus are expected to hold more species. In contrast, clade size and clade age are largely uncorrelated for clades with no major key innovations, because the dominant process is diversity dependence, where ecological limits set the maximum clade size irrespective of clade age. Only for clades that are too young to experience diversity dependence, there may still be a noticeable correlation, even without major key innovations. In phylogenetic studies, often a group of relatively similar species is selected, for example, a genus, or perhaps Vincent and Brown’s (2005) G function (i.e., species with the same Bauplan), which causes a bias toward clades with few, if any, major key innovations. We might hypothesize that major key innovations occur between G functions but not within them.

In an attack on Rabosky’s (2009) ecological limits hypothesis, Wiens (2011) argues that ecological limits to diversification, if they exist at all, should not be seen as an alternative paradigm explaining diversity patterns but that the ecological limits hypothesis is simply part of the traditional idea that ecology influences evolution, and hence, ecological limits should affect diversification rates. This is exactly what our model does (as in fact does Rabosky’s model): the speciation rate is made dependent on clade diversity, but occasionally this dependence is avoided by key innovations.

In summary, we propose a new perspective on adaptive radiations in general and key innovations in particular, in the light of diversity-dependent diversification. Clade-specific key innovations result mostly from a decoupling of the diversity dependence of a focal subclade from the main clade, which may or may not be accompanied by subclade-specific shifts in intrinsic speciation and extinction rates. Clade-wide shifts in speciation or extinction rates, or in the clade-level carrying capacity are usually driven by external factors and offer new opportunities for all species alike, resulting in an adaptive radiation in the whole clade. We have shown how these scenarios can be translated mathematically to compute the likelihood, which can be used for statistical detection of clade-wide and subclade-specific parameter shifts and estimation of the time at which they took place. Diversity dependence is thus key to understanding adaptive radiations.

Enhancements