Replicate radiations provide powerful comparative systems to address questions about the interplay between opportunity and innovation in driving episodes of diversification and the factors limiting their subsequent progression. However, such systems have been rarely documented at intercontinental scales. Here, we evaluate the hypothesis of multiple radiations in the genus Lupinus (Leguminosae), which exhibits some of the highest known rates of net diversification in plants. Given that incomplete taxon sampling, background extinction, and lineage-specific variation in diversification rates can confound macroevolutionary inferences regarding the timing and mechanisms of cladogenesis, we used Bayesian relaxed clock phylogenetic analyses as well as MEDUSA and BiSSE birth–death likelihood models of diversification, to evaluate the evolutionary patterns of lineage accumulation in Lupinus. We identified 3 significant shifts to increased rates of net diversification (r) relative to background levels in the genus (r = 0.18–0.48 lineages/myr). The primary shift occurred approximately 4.6 Ma (r = 0.48–1.76) in the montane regions of western North America, followed by a secondary shift approximately 2.7 Ma (r = 0.89–3.33) associated with range expansion and diversification of allopatrically distributed sister clades in the Mexican highlands and Andes. We also recovered evidence for a third independent shift approximately 6.5 Ma at the base of a lower elevation eastern South American grassland and campo rupestre clade (r = 0.36–1.33). Bayesian ancestral state reconstructions and BiSSE likelihood analyses of correlated diversification indicated that increased rates of speciation are strongly associated with the derived evolution of perennial life history and invasion of montane ecosystems. Although we currently lack hard evidence for “replicate adaptive radiations” in the sense of convergent morphological and ecological trajectories among species in different clades, these results are consistent with the hypothesis that iteroparity functioned as an adaptive key innovation, providing a mechanism for range expansion and rapid divergence in upper elevation regions across much of the New World.
Ever since Darwin's seminal observations of Galá pagos finches (Grant 1986), radiations have offered some of the most spectacular and best-documented examples of diversification, providing uniquely powerful insights into how and why species diverge (Schluter 2000; Linder 2008; Gavrilets and Losos 2009; Losos 2010). This long interest is matched by an equally lengthy debate about what constitutes a radiation, with varying degrees of emphasis on species richness, rates of diversification, evolution of adaptive forms, and/or colonization of novel habitats (Givnish 1997; Sanderson 1998; Schluter 2009; Glor 2010). In many cases, radiations involve both high rates of lineage accumulation as well as ample circumstantial evidence that much of this diversity is of adaptive significance (Losos and Miles 2002; Gavrilets and Vose 2005; Gavrilets and Losos 2009; but see Rundell and Price 2009). However, aside from a few well-explored model groups from island and lake ecosystems (e.g., Hawaiian silverswords, Baldwin and Sanderson 1998; threespine sticklebacks, Rundle et al. 2000; east African cichlids, Seehausen 2006; and Anolis lizards, Losos 2009), mechanisms of diversification in most radiations remain poorly understood, especially in species-rich continental lineages, for which very little is known about the dynamics and likely causes of cladogenesis (Losos 2010).
Over the last decade, a growing number of recent and rapid plant radiations have been documented across widespread geographical areas and diverse habitats (e.g., Richardson et al. 2001a, 2001b; Klak et al. 2004; Hughes and Eastwood 2006; Moore and Donoghue 2007; Tank and Olmstead 2008; Guzmán et al. 2009; Marazzi and Sanderson 2010; Valente et al. 2010). These radiations have generally been attributed to ecological opportunity and/or evolutionary innovation, but relatively few of these studies have identified multiple diversification rate shifts that coincide with the evolution of derived traits and/or ecological or biogeographic transitions among closely related species. Such correlations assigned to specific nodes on a phylogeny offer the potential to pinpoint replicate radiations in time and space, reconstruct their ancestral characteristics and geographic locations, and hence infer the likely triggers and underlying causes of lineage accumulation (Ackerly et al. 2006; Moore and Donoghue 2007; Marazzi and Sanderson 2010). This is important if we are going to be able to discover why some lineages diversify and others do not and the extent to which this is attributable to synchrony—or lack of it—between trait evolution and opportunity.
There is also lack of consensus about whether radiations are contingent on geohistorical circumstances or follow deterministic trajectories (Losos et al. 1998; Taylor and McPhail 2000; Streelman and Danley 2003; Gavrilets and Vose 2005; Melville et al. 2006; Moore and Donoghue 2007; Linder 2008; Gavrilets and Losos 2009; Rundell and Price 2009; Young et al. 2009; Losos 2010), beyond the generally accepted idea that most radiations are characterized by an initial rapid burst of diversification (Gavrilets and Vose 2005; McPeek 2008; Phillimore and Price 2008; Rabosky and Lovette 2008a, 2008b; Mahler et al. 2010; but see Harmon et al. 2010). Fundamental questions remain about what factors dictate the outcomes of radiations in terms of possible limits and hence diversification rate slowdowns, imposed by carrying capacity as defined by geographic area and/or niche heterogeneity (Rabosky 2009a, 2009b, 2010) or even how relevant and appropriate such limits—and hence equilibrium models—are at all (Benton and Emerson 2007).
Ongoing development of methods for inferring the mode and tempo of speciation at macroevolutionary scales has provided new opportunities to address questions about the mechanisms and extent of among-lineage variation in diversification rates (e.g., Hunter 1998; Magallón and Sanderson 2001; Rabosky 2006a; Rabosky et al. 2007; Bokma 2008a; Alfaro et al. 2009; Wertheim and Sanderson 2010; Silvestro et al. 2011). Despite these advances, empirical studies may be confounded by incomplete sampling, nonrandom sampling, and/or background extinction, and uncertainty surrounding divergence time estimates (Nee et al. 1994; Pybus and Harvey 2000; Nee 2006; Ricklefs 2007; Bokma 2008b; Cusimano and Renner 2010; Wertheim and Sanderson 2010; Brock et al. 2011).
We investigate some of these questions via new phylogenetic and diversification rate analyses of the plant genus Lupinus (Leguminosae), a group for which there is evidence to suggest multiple radiations across broad areas of the New World. Lupinus comprises approximately 267 annual and perennial species distributed across lowland and montane environments (Fig. 1), with major concentrations of diversity in North America and South America, as well as a much smaller number of species in the Mediterranean region of Europe and northern Africa. The possibility that the genus has experienced multiple radiations was indicated by 2 studies describing exceptionally recent and rapid species diversification in the Andes (Hughes and Eastwood 2006) and western North America (Drummond 2008).
In the first of these studies, a time-calibrated phylogeny based on nuclear DNA (nDNA) data suggested that following recent divergence from Mexico and Central America, the Andean species of Lupinus diversified at net rates on the order of r = 2.5–3.7 lineages/ myr (Hughes and Eastwood 2006). This radiation is comparable to the fastest rates of net diversification that have been documented for any group of plants (e.g., Richardson et al. 2001a, 2001b; Klak et al. 2004; Scherson et al. 2008; Guzmán et al. 2009; Valente et al. 2010) and is especially remarkable for the continental scale over which it occurred. Estimates based on chloroplast DNA (cpDNA) found even higher rates of net diversification for the clade comprising the Andean and Mexican species (r = 1.4–5.7), with evidence for similar rates (r = 2.0–5.9) in an independently derived lineage of North American species (Drummond 2008).
Each of these rapidly diversifying clades is characterized by the prevalence of perennials that tend to occupy montane habitats, exhibit a diversity of phenotypes (prostrate/acaulescent perennial herbs to woody shrubs/trees), and occupy a very wide array of ecosystems (coastal dunes, chaparral, sagebrush steppe, and grasslands but especially open mountain forests, meadows, and disturbed slopes, extending to subalpine/alpine elevations) across approximately 100 of latitude from approximately 30S in the Andes to approximately 70N in Alaska (Figs. 1 and 2). This suggests that the derived evolution of iteroparity (i.e., reproducing more than once per lifetime) from semelparity (i.e., reproducing only once per lifetime) could have functioned as a key innovation for the invasion of higher elevations throughout the New World (Drummond 2008), providing an evolutionary pathway for adaptive diversification in novel and heterogeneous habitats that were previously inaccessible (Hughes and Eastwood 2006).
Here, we test the hypothesis of multiple radiations in Lupinus. Specifically, we propose that the derived evolution of perennial life history is correlated with the ability to colonize montane environments that were unsuitable for ancestral annual lineages, functioning as an adaptive “key innovation” that facilitated rapid range expansion and subsequent diversification. To address these questions, we assembled approximately 9 Kb of nDNA and cpDNA sequence data representing the most complete phylogenetic sampling of the genus to date (cf. online Appendices 1–4, doi:10.5061/dryad.17rc2f69). After conducting a series of Bayesian Markov chain Monte Carlo (MCMC) relaxed clock phylogenetic analyses to estimate the biogeographic history of diversification in Lupinus, we used 2 recently developed maximum- likelihood (ML) methods that accommodate nonrandom/incomplete sampling and extinction to: (i) infer diversification rate shifts (MEDUSA, Alfaro et al. 2009) and (ii) test for correlations between character states and diversification rates (BiSSE, Maddison et al. 2007; FitzJohn et al. 2009).
MATERIALS AND METHODS
Supermatrix Assembly and DNA Sequence Data
DNA sequences for 5 nuclear genes (GPAT1, GPAT2, ITS1+2, LEGCYC1A, and LEGCYC1B), 2 chloroplast genes (matK and rbcL), and 4 noncoding chloroplast regions (trnL intron, trnL–trnF, trnS–trnG, and trnT–trnL) were obtained from available GenBank accessions for Lupinus and previously unpublished sequences (online Appendix 1). The final data included 122 species of Lupinus (comprising 45.7% of ca. 267 species), augmented with 17 outgroup species from other genera of Genisteae. Each locus was aligned in MUSCLE 3.7 (Edgar 2004) and manually adjusted in Mesquite 2.72 (Maddison and Maddison 2009) to account for microstructural features, such as insertions, deletions, and inversions (Kelchner 2000). After removal of ambiguously aligned regions with uncertain homology, indels were scored in SeqState 1.4.1 (Müller 2005) using simple binary gap coding (Müller 2006). When the data included multiple sequences from different individuals of the same species, we created a single consensus sequence for that species using standard International Union of Pure and Applied Chemistry ambiguity coding to represent intraspecific nucleotide polymorphism (e.g., *bib63; McMahon and Sanderson 2006; Jones et al. 2007; Hughes and Page 2007). Polymorphic indels were coded as missing data. The advantages and limitations of this supermatrix approach using merged consensus sequences concatenated across loci are discussed in relation to alternative coalescent species-tree models in online Appendix 1. The data matrix itself is available online at TreeBASE accession S10651 (http://www.treebase.org) and as a BEAST XML input file in online Appendix 2.
Phylogenetic and Biogeographic Inference
The joint posterior distribution of topologies, divergence times, and ancestral geographic ranges for Lupinus were estimated using Bayesian MCMC searches in BEAST 1.5.4 (Drummond and Rambaut 2007) under a lognormal relaxed clock model (Drummond et al. 2002, 2006). We considered both pure-birth and birth–death coalescent tree priors (Gernhard 2008) to evaluate the potential effect of background extinction on tree shape, as well as incomplete sampling under the same models (Stadler 2009; but see Hartmann et al. 2010). These analyses incorporated a continuous-time Markov chain (CTMC) phylogeographic model with stochastic search variable selection (Lemey et al. 2009). Prior distributions for migration rates were defined by inverse great circle distances between the centroids of 6 primary areas of endemism (i.e., Old World, eastern North America, eastern South America, western North America, Mexico/Central America, and western South America, Fig. 2) based on distributional data assembled from approximately 15,000 georeferenced collections of Lupinus (online Appendix 4). The CTMC phylogeographic model explicitly integrates the process of phylogenetic inference with reconstruction of ancestral biogeography (Lemey et al. 2009) but assumes that ancestral ranges are limited to single regions (Sanmartín et al. 2008; Lamm and Redelings 2009; Ree and Sanmartín 2009). Since this assumption is true for all extant species of Lupinus (discounting anthropogenic invasives), we believe that the CTMC model represents a reasonable alternative to dispersal-vicariance parsimony (Ronquist 1997) or dispersal-extinction cladogenesis (Ree et al. 2005; Ree and Smith 2008).
Previous phylogenetic studies have recovered strong support for the placement of Lupinus among the genistoid legumes (i.e., Genisteae, cf. Wojciechowski et al. 2004; Lavin et al. 2005). Based on these data, rates of nucleotide substitution and indel evolution were calibrated to an absolute time scale by placing normal prior distributions on the mean and 95% intervals for the most recent common ancestor (MRCA) of the Genisteae minus Lupinus (12.8 Ma, CI = 6.3–18.0 Ma) and the stem node of Lupinus, that is the MRCA of the Genisteae (18.8 Ma, CI = 11.9–24.9 Ma), approximating estimates obtained from a dense fossil-calibrated phylogeny of the Leguminosae (Simon et al. 2009). Posterior sampling was conducted every 10,000 iterations, with at least 50,000,000 posterior iterations after 5,000,000 burn-in iterations from a random starting tree (An example XML input file for BEAST is included in online Appendix 2.). Convergence and stationarity were verified from the results of at least 3 MCMC searches, accompanied by monitoring of trace plots and effective sample sizes in Tracer 1.4.1 (Rambaut and Drummond 2008). The combined results from replicate runs were summarized in TreeAnnotator 1.4.8 (Drummond and Rambaut 2007) as the maximum clade credibility tree (MCCT) with median branch lengths derived from the posterior distribution.
Unlinked nucleotide substitution models for each locus were selected in MrModeltest 2.3 (Nylander 2008) based on optimal Akaike information criterion (AIC) scores (Posada and Buckley 2004) for each partition: HKY + Γ (LEGCYC1B); GTR + Γ (GPAT1, GPAT2, LEGCYC1A, matK, trnL intron, trnL–trnF, and trnT–trnL); GTR + I + Γ (ITS1 + 2, rbcL, and trnS–trnG). To avoid over-parameterization due to nonindependence of estimates for among-site rate variation and the proportion of invariable sites (Yang 2006), we favored the less complex GTR + Γ model over GTR + Γ + I for ITS1 + 2, rbcL, and trnS–trnG. Unlinked indel matrices were handled under a discrete Markov model for each locus, including a correction for ascertainment bias due to exclusive sampling of variable characters (Lewis 2001).
To assess the effect of birth–death and Yule coalescent tree priors on the posterior distribution of topologies and divergence times, we calculated Bayes factors (Newton and Raftery 1994) as the difference between the harmonic mean of marginal log-likelihood (lnL) scores from MCMC runs performed under each model. To evaluate the support for infrageneric groups used in the MEDUSA diversification rate analyses, we also considered the difference in Bayes factors between unconstrained and constrained MCMC searches where the monophyly and phylogenetic positions of terminal clades were restricted to the topology recovered in the MCCT summary tree. Because the raw harmonic mean can be an unstable estimator for Bayes factors, estimates were obtained using importance sampling in Tracer with 1000 bootstrap replicates (Suchard et al. 2005). Twice the difference in Bayes factors (BF01= 2(lnL0− 2lnL1)) on the order of BF01= 1–3 is considered negligible, BF01= 3–20 suggests positive support for B0, BF01= 20–150 indicates strong support for B0, and BF01> 150 shows very strong support for B0 (Kass and Raftery 1995; Raftery 1996).
Diversification Rate Analyses
Among-lineage variation in diversification rates was examined using MEDUSA (Alfaro et al. 2009), an extension to the birth–death likelihood model (Rabosky 2006a; Rabosky et al. 2007) that allows rate shifts to be inferred at one or more topological positions. To avoid the problem of unresolved incomplete and/or nonrandomly sampled lineages, MEDUSA requires a time- calibrated phylogeny in which extant taxonomic diversity can be assigned to monophyletic terminal clades, using a stepwise procedure to evaluate the support for increasingly complex models of diversification based on the difference in AIC scores (ΔAIC) as extra parameters are added to describe one or more rate shifts. Because the designation of an arbitrary cutoff for a “significant” change in ΔAIC can influence the inferred number of shifts (Alfaro et al. 2009), we instead used the lowest AIC score to select the most strongly weighted rate shift model. This approach is analogous to maximizing the likelihood while including a penalty for additional parameters (Burnham and Anderson 2002), and we report the marginal distribution of ΔAIC scores conditioned on 1000 subsampled posterior trees from the constrained pure-birth BEAST analysis.
To estimate levels of species richness in Lupinus clades, taxonomic data were compiled from a range of sources, including published floristic treatments, herbarium records, and a critical appraisal of the ILDIS legume database (Roskov et al. 2005). Based on the maximum clade credibility tree, we identified 19 well-supported infrageneric lineages for which extant species richness could be assigned to monophyletic terminal clades, corresponding to morphologically and/or geographically recognized groups (online Appendix 3). MEDUSA analyses were performed on a randomly selected subset of 1000 trees from the posterior distribution of the constrained pure-birth BEAST analysis using the R packages ape 2.3 (Paradis et al. 2004), Laser 2.2 (Rabosky 2006b), and Geiger 1.3-1 (Harmon et al. 2008b).
Subsampled posterior trees were pruned to a single branch for each terminal clade, and the outgroup was removed leaving Lupinus rooted. Per lineage rates of net diversification (r = birth - death) and relative extinction (ε= death/birth) were estimated under a birth–death and a pure-birth model. Following preliminary runs of MEDUSA, we adjusted the optimization interval in the source code for Geiger to avoid local optima defined by the default limits placed on the upper and lower bounds of the likelihood surface. The relative fit of alternative diversification rate models was evaluated based on the distribution of ΔAIC scores using the 1000 subsampled trees described above. Because MEDUSA uses stem ages for estimates of net diversification rates in unresolved terminal clades, we also calculated this parameter in Geiger under a pure-birth model with extant diversities calibrated to the crown node age of selected infrageneric groups (Magallón and Sanderson 2001).
Ancestral States and Correlates of Diversification
We examined 2 ecological characters hypothesized to be correlated with increased rates of diversification in Lupinus (Drummond 2008): (i) life history (semelparous vs. iteroparous) and (ii) habitat (lowland vs. montane). Data for life history were obtained from the taxonomic literature, georeferenced herbarium records, and personal field observations (online Appendices 3–4). Biennials and short-lived herbaceous perennials that reproduce only once or twice were coded as semelparous for the purpose of this study, although vegetative growth may persist across a limited number of growing seasons. Habitat was coded employing a compound metric incorporating both elevation and rugosity (a measure of topographic complexity such that = Ar/Ag, where Ar= true surface area and Ag= geometric or planar surface area). The motivation for using this metric was to capture the combined effects of environmental conditions associated with higher elevations, as well as the physiographic and ecological heterogeneity of montane habitats (cf. Körner 2007).
Digital elevation models (DEM) at 30 arc second resolution (ca. 1 km2) for each continent were obtained from GTOPO30 radar imaging data (United States Geological Survey 1996), imported into ArcMap 9.3 (ESRI 2008), and reprojected to Lambert azimuthal equal area format, allowing to be estimated using Benthic Terrain Modeler (Wright et al. 2005). After visual inspection of mapped onto each DEM, we used a map algebra filter to derive a new raster layer in which cells were coded as montane based on a threshold of 1.01, excluding all cells <500 m to remove lowland coastal foothills and desert canyons with topographically complex terrain, and including all cells >2500 m to retain alpine/ subalpine valleys and plateaus, yielding a well- characterized division between mountainous upper elevations and less rugged lower elevations. We then mapped latitude/longitude coordinates compiled for approximately 15,000 specimens of Lupinus to designate each species as lowland or montane based on the majority proportion of occupied habitat (online Appendix 4). For species with few georeferenced specimen data, we relied on data from the taxonomic literature, herbarium records, and personal field observations.
Ancestral state reconstructions were performed using a series of Bayesian reversible-jump hyperprior (RJHP) MCMC analyses in BayesTraits 1.0 (Pagel et al. 2004; Pagel and Meade 2006) on 1000 subsampled posterior trees from the unconstrained pure-birth BEAST analysis. This approach allowed us to obtain the posterior distribution of ancestral character states while accounting for phylogenetic uncertainty, with the results weighted among competing models for equal versus unequal transition rates between character states. Prior probabilities for discrete rate parameters were determ-ined by setting a hyperprior drawn from a uniform distribution of exponential means approximating the range of values obtained from preliminary ML estimates in BayesTraits.
To test for correlated evolution between life history and habitat, we used RJHP MCMC in BayesTraits to estimate posterior support for dependent versus independent models of state changes between characters. For each model, Bayes factors were calculated in Tracer (see above) to evaluate the relative fit of models, including an additional free parameter (κ) for rescaling branch lengths. When κ= 1.0, traits evolve linearly with branch lengths, whereas values of κ< 1.0 suggest changes are more likely to occur along shorter versus longer branches (Pagel 1999). For each BayesTraits analysis, at least 3 replicate MCMC runs were conducted for 100,000,000 iterations following a burn-in of 10,000,000 iterations, with posterior sampling every 10,000 iterations after adjusting the ratedev tuning parameter to ensure adequate mixing.
The BiSSE binary state speciation and extinction model (Maddison et al. 2007) implemented in diversitree 0.4-5 (FitzJohn et al. 2009) was used to examine whether life history and/or maximum elevation are directly correlated with differential rates of diversification. BiSSE employs ML optimization to estimate absolute rates of asymmetric character change (q), speciation (λ), and/or extinction (μ) by maximizing the likelihood of these parameters for a given topology with branch lengths (Maddison et al. 2007). Note that this parameteri- zation differs from MEDUSA, which uses likelihood equations based on net diversification (r) and relative extinction (ε) rather than absolute rates.
For each character, we compared 2 models in BiSSE, both times setting μ= 0 based on the weight of evidence from the BEAST and MEDUSA results, namely model (i) an unconstrained model in which q and λ were allowed to vary and model (ii) a constrained model in which λ for each character was set equal. If diversification rates are correlated with character states, then the unconstrained model should be favored over the constrained model. Because the likelihood ratio tests based on the approximation can have low power with sample sizes <500 species (Maddison et al. 2007), we instead calculated ΔAIC scores to evaluate the fit of alternative models. The distribution of ΔAIC scores between unconstrained and constrained BiSSE models was based on 1000 subsampled posterior trees from the unconstrained pure-birth BEAST analysis, with extant diversities of incompletely sampled clades assigned as in the MEDUSA analyses (Fig. 3) and the proportion of character states estimated from the taxonomic literature and herbarium records (online Appendices 3–4).
Phylogeny and Biogeography
We recovered a well-resolved phylogeny with strong support from posterior probabilities (PP > 0.90) along backbone nodes that define geographically structured relationships among key groups in Lupinus (Fig. 2; online Fig. S1). The genus is divided into 3 principal lineages descended from an Old World MRCA (95% highest posterior density = 7.2–17.7 Ma): (i) the Mediterranean and north African species (4.6–12.5 Ma) sister to the unifoliolate species from eastern North America (0.1–2.4 Ma); (ii) the eastern South American species (2.3–7.1 Ma) sister to the Texas bluebonnets from eastern North America (0.1–2.3 Ma); and (iii) the western New World species (5.0–13.2 Ma), comprising the Andean and Mexican species (1.19–3.50 Ma) derived from a paraphyletic assemblage of western North American species (2.1–5.5 Ma). The Old World/eastern North America unifoliolate clade (5.6–14.8 Ma) is sister to a large clade comprising the remaining New World species derived from a MRCA in eastern North America (6.0–15.3 Ma). The MRCA of the eastern South American/eastern North American Texas bluebonnet clade (2.8–9.6 Ma) was also placed in eastern North America, consistent with dispersal to eastern South America. Nonetheless, posterior support for ancestral ranges around the perimeter of the Atlantic basin was not decisive (Fig. 2). In contrast, ancestral area reconstruction across the western New World clade is largely consistent with a range expansion into contiguous geographic regions, progressing south from western North America into Mexico and the Andes (Fig. 2), although there is some uncertainty regarding the ancestral area of the Mexican and Andean clades.
We found strong support (PP > 0.90) for the phylogenetic placement of 19 terminal clades corresponding to morphologically and/or biogeographically distinct infrageneric groups (Fig. 3; cf. online Appendix 3). Although there was considerable uncertainty regarding the sister relationship between the Mexican and Andean clades (PP = 0.56), Bayes factors indicated that the constrained topology (Fig. 2) was favored over a model without constraints (BF01= 19.02), suggesting little (if any) loss of information regarding phylogenetic uncertainty. Likewise, Bayes factors also recovered strong support in favor of a pure-birth model with no extinction (BF01= 86.07). Accounting for incomplete sampling (i.e., assuming taxa have been randomly sampled from extant diversity—which is not the case here, cf. Fig. 3) did not appreciably change the posterior distribution of topologies, branching times, or support for clades under either model (results not shown). Accepting the clades recovered from the constrained pure-birth BEAST analysis as representative of the best-resolved terminal branches to which we could confidently assign extant diversities, we estimated the proportions of extant and sampled taxa for each of these lineages (Fig. 3; cf. online Appendix 3).
ΔAIC scores comparing the results of pure-birth versus birth-death models in MEDUSA favored the pure-birth model in 100% of the 1000 subsampled topologies (95% interval of ΔAIC values from ML estimates conditioned on the posterior distribution of trees = 4.81–10.08; cf. online Fig. S2). Under this model, MEDUSA consistently detected 3 shifts in net diversification rates (Fig. 3; online Fig. S2) relative to background levels in Lupinus (r = 0.18–0.48): (i) a primary shift in the montane regions of western North America (r = 0.48–1.76; ΔAIC= 43.68–116.12); (ii) a secondary shift associated with colonization of the Mexican highlands and Andes (r = 0.89–3.33; ΔAIC= 0.00–15.91); and (iii) a third shift in the lowlands of eastern South America (r = 0.36–1.33; ΔAIC= 0.09–18.18). Similar results were obtained under the birth–death model allowing for extinction (cf. online Fig. S3), although the eastern South American clade exhibited markedly lower rates of net diversification (r = 0.19–0.32) with high levels of background extinction (ε= 0.49–0.99), whereas the Mexican/Andean clade had slightly lower rates of net diversification (r = 0.26–3.03) with a bimodal distribution for relative extinction centered on the extremes of zero versus equal rates (ε= 0.00–0.99). In contrast, diversification rates estimated in Geiger using crown node ages under a pure-birth model were somewhat more rapid than those derived from stem node ages estimated in MEDUSA (Fig. 3): (i) western North America (r = 0.74–1.97); (ii) Mexico (r = 1.19–6.15); (iii) Andes (r = 1.56–5.21); and (iv) eastern South America (r = 0.39–1.22).
Correlates of Diversification
The major transition from lowland/annual to perennial/montane is strongly supported (PP > 0.90) along the same branch where the shift to higher rates of net diversification occurred in western North America (Fig. 4), including a secondary rate shift among the perennial/ montane species of Mexico and the Andes. In contrast, the rapidly diversifying eastern South American species are mostly perennials restricted to lower elevations (Fig. 4). RJHP MCMC results from BayesTraits showed very strong support (PP = 1.0) for correlated transitions between life history and montane habitat. These results also indicated positive support for κ< 1.0 (BF01= 8.86, κ= 0.00–0.39), consistent with higher probabilities of change along shorter branches for both characters. The distribution of ΔAIC scores from BiSSE (Fig. 4; cf. online Fig. S4) showed strong support for a correlation between higher rates of diversification for perennial life form (r = 0.93–2.77; ΔAIC= 62.78–125.38) and montane occurrence (r = 0.92–2.64; ΔAIC= 57.97–109.89) versus annual life form (r = 0.17–0.42) and lowland occurrence (r = 0.21–0.51).
Rapid Diversification and Multiple Radiations
Geotemporal patterns of species diversification remain poorly understood in species-rich continental clades, for which there is still limited insight into the dynamics and likely causes of diversification (Losos 2010). Here, we present evidence for a series of recent and rapid continental radiations in Lupinus, demonstrating a robust correlation between a derived life history trait, montane occurrence, and net diversification rates. Multiple continental radiations of this sort are likely to be much more common than currently documented for widely distributed species-rich plant genera (e.g. Pelargonium, Bakker et al. 2005; Astragalus, Scherson et al., 2008; Castilleja, Tank and Olmstead 2008, 2009; and Indigofera, Schrire et al. 2009), offering excellent scope to address fundamental questions about the patterns and processes underlying radiations.
The largest of the Lupinus radiations spans montane regions across much of the western New World (196 spp., r = 0.48–1.76). This “super radiation” encompasses a set of 3 separate nested radiations in the montane habitats of western North America (58 spp., r = 0.74–1.96), Mexico/Central America (46 spp., r = 1.19–6.15), and the Andes (81 spp., r = 1.56–5.21). Although rates of net diversification in these clades are comparable with estimates from previous studies (Hughes and Eastwood 2006; Drummond 2008; Silvestro et al. 2011), we also found evidence for accelerated diversification in the lower elevation tropical, subtropical, and subtemperate grassland and campo rupestre habitats of eastern South America (35 spp., r = 0.39–1.22). Crown node rates across all 4 of the most diverse terminal clades (Fig. 3) are approximately 2–12 times faster than background levels in Lupinus (r = 0.18–0.48), with the Mexican and Andean radiations exhibiting rates on a par with Dianthus (r = 2.2–7.6), which is perhaps the fastest documented radiation of plants to date (Valente et al. 2010).
Although rate shifts were detected at only 3 internal nodes (Fig. 3), the statistical power of the MEDUSA likelihood model to detect shifts among recently diverged lineages is limited by reliance on stem node dates (Alfaro et al. 2009). In the absence of reliable estimates of crown node dates for incompletely sampled or poorly resolved lineages, methods such as MEDUSA will be unable to diagnose parallel radiations in terminal sister clades with similar rates of lineage accumulation since there is no rate shift per se to be inferred beyond the node representing the MRCA of equally diverse clades. This highlights a fundamental problem with using rate shifts to identify parallel radiations that seem obvious when viewed in terms of concordance between phylogenetically and biogeographically structured levels of species richness (e.g., Mexican and Andean radiations of Lupinus; cf. Fig. 3).
We found no evidence for background extinction, with either the BEAST or the MEDUSA analyses indicating strong support for a pure-birth model. Although extinction undoubtedly plays an important role in shaping among-lineage variation in species richness (Stanley 1979; Sepkoski 1998; Benton and Emerson 2007; Alroy 2008; Rabosky 2009c; Quental and Marshall 2010), in the absence of data from the fossil record—as is the case for Lupinus—estimates of relative extinction rates from phylogenies based on extant taxa may be inherently problematic (Hunter 1998; Crisp and Cook 2009; Rabosky 2009c, 2010; Quental and Marshall 2010). Despite these limitations, support for a pure-birth model suggests that loss of lineages has been relatively inconsequential with respect to the observed phylogenetic patterns of extant species richness in Lupinus. This is not to say that extinction has not occurred but rather that the available data and models were inconsistent with extinction as a primary determinant of species richness over the relatively recent time scales examined here.
Conversely, estimates of net diversification rates could be confounded when speciation is density dependent (Rabosky 2009a, 2009b, 2010). The latter scenario represents an intriguing alternative hypothesis, and if true, would suggest that Lupinus has repeatedly experienced ecological release following colonization of new regions rather than undergoing changes in diversification rates per se. In this sense, the various New World radiations of Lupinus might simply represent different levels of steady-state diversity in response to carrying capacities at equilibrium between constant rates of speciation and extinction (Cracraft 1985; Losos and Schluter 2000; Rabosky 2009a, 2009b, 2010; Rabosky and Glor 2010).
The development of new models for densitydependent diversification provides promising directions for future research on these questions (Rabosky and Lovette 2008a, 2009; but see Bokma 2009; Quental and Marshall 2010). However, the patterns observed here do not exhibit the expected reduction in diversification rates commonly exhibited in adaptive radiations, whereby an “early burst” of cladogenesis is followed by stasis as ecological niches are filled (McPeek 2008; Phillimore and Price 2008; Rabosky and Lovette 2008a, 2008b; Mahler et al. 2010; but see Cusimano and Renner 2010; Harmon et al. 2010; Brock et al. 2011). In part, this may be a function of the differences between young versus mature radiations (Linder 2008), and we suggest that the striking absence of time dependency across the genus (i.e., older clades are not more diverse than younger clades) could be more consistent with the early stages of accelerated diversification in derived lineages (but see Rabosky 2009b).
Our results also demonstrate the importance of accounting for incomplete and/or biased taxon sampling in macroevolutionary studies of diversification (Marazzi and Sanderson 2010; Brock et al. 2011; Ryberg and Matheny 2011). Because previous comparative phylogenetic analyses of Lupinus (Hughes and Eastwood 2006; Drummond 2008; Moore and Donoghue 2009; Silvestro et al. 2011) were restricted by more limited availability of DNA sequence data and appropriate methodologies, key shifts, and correlates of diversification were not necessarily recognized. For example, in contrast to the results presented here, Moore and Donoghue (2009) found no evidence for a shift to increased rates of diversification following the colonization of eastern South America by Lupinus, we believe in large part because their Bayesian cross- validation predictive test was applied to a biased selection of extant diversity. Similarly, Silvestro et al. (2011) postulated only 2 rate shifts located on different branches of the tree from the 3 rate shifts found in this study, in part because their analysis was based on more limited data that reconstructs a different tree topology from that presented here. The issue of nonrandom phylogenetic sampling is distinct from the “pull of the present” when sampling is incomplete (but effectively random) or when background extinction leaves an apparent surplus of recently diverged lineages (Pybus and Harvey 2000; Crisp and Cook 2009; Brock et al. 2011). It is increasingly clear for species-rich taxa spanning broad geographical areas that dense taxon sampling will be required to accurately detect rate shifts and assign them to particular nodes, with implications for the inferred locations, divergence times, and ancestral trait combinations that are needed to provide insights into the factors driving species radiations (Ackerly 2000; Heath et al. 2008; Cusimano and Renner 2010; Marazzi and Sanderson 2010; Ryberg and Matheny 2011).
Mechanisms of Diversification
The biogeographic and ecological history of radiations in Lupinus coincides with several major physiographic changes in the New World. Expansion into western North America was initially limited to annual lineages at lower elevations centered in California, followed by rapid colonization of higher elevations by perennial lineages (Fig. 5), contemporaneous with Pliocene uplift of the Sierra Nevada, Coastal Ranges, Transverse Ranges, and Cascades in the past approximately 5 myr (Graham 1999; Harden 2004) and invasion of the neighboring Rocky Mountains, which had already experienced final uplift during the Laramide orogeny ending approximately 40 Ma (Dickinson et al. 1988; McMillan et al. 2006). The inferred timing of this transition is also coincident with the separate expansion of a small clade of perennial species into boreal regions of Alaska/Canada and lowland open forests of eastern North America, likely proceeding into the highlands of Mexico and Central America (Fig. 2, online Appendix 3 and Fig. S1), which were largely in place as early as approximately 15 Ma (Ferrari et al. 1999; Ferrusquía-Villafranca and González-Guzmán 2005). The most recent burst of diversification followed the uplift of high mountains in Costa Rica and Panama associated with the formation of the Isthmus of Panama (Coates and Obando 1996) and especially final uplift of the northern Andes and emergence of significant areas of upland habitats suitable for the establishment of Lupinus in the last approximately 4 myr (Gregory-Wodzicki 2000; Hughes and Eastwood 2006; Garzione et al. 2008; Ehlers and Poulsen 2009).
Life history theory predicts that semel- parous species are expected to perform better than iteroparous species in drier environments with greater seasonality and interannual fluctuation in water availability (Cole 1954; Charlesworth 1994; Roff 2002). Likewise, numerous studies have shown that annual species of plants are generally favored in warmer and drier climates, whereas perennial species are favored in cooler and wetter environments (e.g., Schaffer and Gadgil 1975; Mulroy and Rundel 1977; Evans et al. 2005; Smith and Beaulieu 2009). The results from BayesTraits and BiSSE indicate that iteroparity and montane habitat are strongly correlated in Lupinus, and both traits are associated with increased rates of diversification (Fig. 4; online Fig. S4). Given that most annual species of Lupinus occupy lower elevation xeric environments with extreme temporal/spatial fluctuations in rainfall and temperature (Fig. 5), the derived evolution of iteroparity may have functioned as a key innovation for the invasion of upper elevation ecosystems where rainfall regimes are moderated by orographic effects and thus more constant (Drummond 2008; cf. Roe 2005; Ehlers and Poulsen 2009). Conversely, studies of other plant genera have shown that transitions from perennial to annual are often accompanied by colonization of xeric lowland environments (e.g., Oenothera, Evans et al. 2005, 2009 and Scorzoneroides, Cruz-Mazo et al. 2009), suggesting that transitions between semelparity and iteroparity may be associated with colonization and/or diversification in novel habitats.
The finding of increased species richness among montane clades of Lupinus mirrors a broader trend across many other taxonomic groups in which higher levels of species richness have been documented at middle to upper elevations (e.g., Rahbek 1995; Fjeldså and Lovett 1997; Körner 2000, 2002; Lomolino 2001; Badgley 2010). Several ideas have been proposed to explain this pattern, including “ecological opportunity” or release from prior environmental constraints (Simpson 1944; Hughes and Eastwood 2006; Moore and Donoghue 2007; Harmon et al. 2008a; Yoder et al. 2010), habitat heterogeneity or the idea of a “montane mosaic” (Janzen 1967; Hughes and Eastwood 2006; *bib14; Kozak and Wiens 2007), “ecological niche conservatism” in which vicariant habitat and nonadaptive isolation by drift lead to divergence within prior adaptive tolerances (Wiens and Graham 2005; Kozak et al. 2006; Kozak and Wiens 2006, 2010), “time- to-speciation” or “montane museum” with earlier colonization allowing longer times for lineage accumulation (Stephens and Wiens 2003; Smith et al. 2007; Wiens et al. 2007), a “species–area relationship” where diversity is controlled by spatial effects such as total available geographic area (MacArthur 1969; Diamond 1977; Lomolino 2000; Losos and Schluter 2000), and/or a “mid-domain effect” in which diversity is determined by the bounds of minimum to maximum accessible elevations (Colwell and Lees 2000; Colwell et al. 2004; Smith et al. 2007).
Given that the greatest levels of species richness are found in younger rather than older clades of Lupinus (Fig. 3), the montane museum hypothesis is clearly inapplicable in this case. A general species–area relationship also seems unlikely since the most diverse clades in Lupinus occupy geographic areas that are not apparently larger than those occupied by the much less diverse eastern North American and Mediterranean/north African lineages (Fig. 2, also see Fig. 1 in Plitmann and Heyn 1984). Aside from debate regarding appropriate null models for examining biogeographically and phylogenetically structured diversity (Colwell et al. 2005; Davies et al. 2005; Zapata et al. 2005), evidence for a mid-domain effect in Lupinus is equivocal. In western North America, species can be found from sea level to subalpine/alpine zones, although primarily concentrated at middle to upper elevations (Fig. 5), whereas species richness in the Andean perennial lineage is concentrated at higher elevations (>2500 m, see Hughes and Eastwood 2006). We are also doubtful about the ecological niche conservatism hypothesis with respect to Lupinus, given the wide diversity of life history strategies and habitats occupied by very closely related species (Fig. 1). Instead, we favor a combination of ecological opportunity and montane mosaics as the main drivers of diversity, with sharp elevational gradients yielding heterogeneous and isolated habitats that allow speciation to proceed rapidly (Hughes and Eastwood 2006).
Contingency, opportunity, and key innovation.—
Given the evidence for multiple radiations in Lupinus, an important question is whether these patterns support the hypothesis that shifts in life history and habitat functioned as a “key innovation” (e.g., Simpson 1953; Hunter 1998; Schluter 2000; Galis 2001; Marazzi and Sanderson 2010) that led to accelerated diversification in Lupinus. Considering the biogeographic scenarios outlined above as well as the strong correlation between montane habitat, iteroparity, and increased rates of lineage accumulation, we propose a multi-tiered scenario that invokes geohistorical contingency, key innovation, dispersal to novel ecosystems, and ecological opportunity as the principal factors leading to the observed patterns of accelerated diversification in Lupinus (Yoder et al. 2010; Moore and Donoghue 2007). Under this scenario, orographic uplift created new island-like environments that could not be fully exploited by ancestral annual lineages in North America. Following the derived evolution of perennial life history (potentially in situ as lowland areas were transformed by tectonic uplift), these species were increasingly able to invade upper elevations (Fig. 5), encountering a heterogeneous array of novel environments across the western New World (Hughes and Eastwood 2006; Drummond 2008). In these highly fragmented montane regions, diversification is expected to proceed rapidly due to the combined influences of geographic isolation and environmental variation, allowing local selection and genetic drift to drive differentiation among populations in the absence of significant gene flow between incipient species (Endler 1977; Barton 1996; Schluter 2000; Emerson and Gillespie 2008; Gavrilets and Losos 2009).
Although the proposed combination of key innovation, ecological opportunity, and “montane mosaic” scenario is compelling, perennial life history might not be causally related to the colonization of montane ecosystems since codistributed characters may have influenced patterns of diversification (Maddison et al. 2007; Moore and Donoghue 2007; Losos 2011) or iteroparity may have evolved secondarily to environmental factors other than montane habitat. Furthermore, the categorical definition of “montane” used here is at best a proxy for other environmental factors that directly shape the ecology of organisms (Körner 2007). Likewise, our definition of “perennial” is predicated on the distinction between patterns of reproductive life history (i.e., iteroparity vs. semelparity—see Materials and Methods section) and ignores covariance in growth form, body mass, and/or life span (Enquist et al. 1999; Niklas and Enquist 2003), any of which might also be related to success in montane environments. Although more sophisticated multivariate quantitative models (FitzJohn 2010) could circumvent these simplifications, simple binary character state definitions were used to analyze correlations between character states and diversification rates because, unlike the multivariate models, BiSSE allows for missing taxa and/or unresolved terminal clades (FitzJohn et al. 2009).
The eastern South American radiation of Lupinus represents a predominantly perennial derived clade largely restricted to mid to low elevation habitats that share few environmental characteristics with the montane regions of western North America and South America (Overbeck et al. 2007; cf. Monteiro and Gibbs 1986; Simon et al. 2009), apparently contrary to the key innovation hypothesis. This may reflect the extent to which our definition of montane is a proxy for environments and overall habitat favoring perennials over annuals. Within this clade, there are still strong indications that annuals occupy lower and drier habitats than perennials, which grow in cooler and more mesic conditions, albeit only at middle elevations here classified as lowland. Likewise, lowland perennial species from western North America tend to occupy relatively mesic coastal or riverine environments. These patterns suggest that the correlation between iteroparity and montane habitat represent a more general phenomenon of ecological adaptation, even though the overwhelming majority of diversification has taken place at upper elevations. In contrast, the relatively depauperate Old World clade of lowland annuals has remained species poor, despite the proximity of montane environments and temperate habitats, where introduced North American perennial lupines are naturalized and invading, further supporting the idea that evolutionary contingency and key innovation have been important factors underlying diversification of Lupinus.
The diversity of morphological and ecological niches occupied by Lupinus (Fig. 1) suggests that adaptive pathways for speciation may have been recapitulated in the multiple radiations outlined here. We currently lack hard evidence for examples of “species-for-species” trait matching based on fine-scale phenotype/environment correlations across independently diversifying clades (i.e., “replicate adaptive radiations” sensuSchluter 1990, 2000; Losos and Miles 2002; Glor 2010; Losos 2010). However, we expect that such patterns will be detected as additional data are assembled even if they will have been tempered by inevitable geohistorical circumstances that define the array of ecological and evolutionary opportunities available at any particular point in time and space (Taylor and McPhail 2000; Moore and Donoghue 2007; Young et al. 2009; but see Losos et al. 1998; Melville et al. 2006), as well as differences in the adaptive landscapes over which species have evolved (Gavrilets and Vose 2005; Gavrilets and Losos 2009).
Support for this project was provided by a postdoctoral fellowship from the Center for Research on Invasive Species and Small Populations at the University of Idaho [C.S.D.]; a Royal Society University Research Fellowship at the University of Oxford [C.E.H.]; a Biotechnology and Biological Sciences Research Council studentship [R.J.E.]; the Genetics Society; and the Stanley Smith Horticultural Trust. Computational support was provided by the Initiative for Bioinformatics and Evolutionary Studies at the University of Idaho [NIH/ NCRR P20RR16448 and P20RR016454].
We thank Luke Harmon, Richard FitzJohn, Peter Linder, Daniel Rabosky, Stephen Smith, Jack Sullivan, and David Tank for critical discussion and/or access to source code; Donovan Bailey, Stephan Beck, Tim Budden, Ruth Clark, Aniceto Daza, Alfonso Delgado, Chris Fagg, Rob Forrest, Martin Gardner, Greg Kenicer, Bente Klitgaard, Gwil Lewis, Helga Ochoterena, David Neill, Arely Palabral, John Pannell, Terry Pennington, Toby Pennington, Carolyn Proenca, Carlos Reynel, Mario Sousa, Richard Spellenberg, Maria Teresa Schifino- Wittmann, Salvador Talavera, John Wood, and Marty Wojciechowski for assistance with fieldwork and provision of plant samples; and the authorities in Ecuador, Peru, Bolivia, and Brasil for permission to collect material. We thank Alexandre Antonelli, an anonymous reviewer, and especially Susanne Renner, for extensive critical feedback that greatly improved the manuscript.