The disparity in species richness among groups of organisms is one of the most pervasive features of life on earth. A number of studies have addressed this pattern across higher taxa (e.g. ‘beetles’), but we know much less about the generality and causal basis of the variation in diversity within evolutionary radiations at lower taxonomic scales. Here, we address the causes of variation in species richness among major lineages of Australia's most diverse vertebrate radiation, a clade of at least 232 species of scincid lizards. We use new mitochondrial and nuclear intron DNA sequences to test the extent of diversification rate variation in this group. We present an improved likelihood-based method for estimating per-lineage diversification rates from combined phylogenetic and taxonomic (species richness) data, and use the method in a hypothesis-testing framework to localize diversification rate shifts on phylogenetic trees. We soundly reject homogeneity of diversification rates among members of this radiation, and find evidence for a dramatic rate increase in the common ancestor of the genera Ctenotus and Lerista. Our results suggest that the evolution of traits associated with climate tolerance may have had a role in shaping patterns of diversity in this group.
Why do some groups of organisms contain vastly greater numbers of species than other groups? A prominent explanation holds that variation in species richness among groups can be explained in part by differences in per-lineage rates of species origination and extinction (Slowinski & Guyer 1989; Mooers & Heard 1997). Many previous studies have tested for variation in diversification rates across large phylogenetic and taxonomic scales, such as angiosperms (Magallon & Sanderson 2001; Sims & McConway 2003; Davies et al. 2004) or passerine birds (Phillimore et al. 2006; Ricklefs 2006). By comparison, few studies have addressed the extent to which diversification rates vary within species-level radiations, and particularly within groups that do not represent classic adaptive radiations (Ruber & Zardoya 2005; Kozak et al. 2006; McKenna & Farrell 2006). Extending the analysis of diversification rates to species-level radiations of non-model groups is important, because recent work suggests that clade age and not diversification rate might be the dominant signal influencing the distribution of diversity among major animal clades (McPeek & Brown 2007). Understanding the contribution of diversification rate variation to differences in species richness requires more data from evolutionary radiations at lower taxonomic scales.
Australian lizards known as sphenomorphine skinks are a particularly appropriate system in which to test the extent to which diversification rates vary among lineages during continental evolutionary radiations. This group comprises at least 232 species that appear to have diverged in situ within mainland Australia (Cogger 2000; Reeder 2003; Skinner in press), making it by far the most diverse vertebrate radiation in Australia, and one of the most diverse continental radiations across all amniotes. Sphenomorphine skinks occur in virtually all terrestrial habitats in Australia, and body elongation and limb reduction have evolved repeatedly within the group (Greer 1989). Of the 15 Australian sphenomorphine genera, the genus Ctenotus alone includes nearly 100 described species (Cogger 2000). Although it has been suggested previously that Ctenotus may have experienced elevated diversification rates (Reeder 2003), there have been no formal tests for heterogeneity of diversification rates within the sphenomorphine skink radiation.
Here, we present the first molecular phylogenetic analysis of major lineages within Ctenotus and use these historical data to test whether it is necessary to invoke among-lineage variation in diversification rates to explain the high disparities in extant diversity among major groups of this radiation. We also test the monophyly of the genus Ctenotus to exclude the possibility that the high species diversity of this genus is simply a taxonomic artefact. Finally, we use a likelihood framework based on the birth–death process to test whether high species diversity within Ctenotus is the outcome of increased diversification rates within the genus, or whether it reflects a decline in diversification elsewhere in the radiation. Our results indicate a dramatic increase in diversification occurring in the lineage leading to Ctenotus and its sister taxon Lerista and suggest that arid-adapted lineages may have diversified explosively in response to the expansion of the Australian arid zone over the past 20 Myr.
2. Material and methods
(a) Taxon sampling
We obtained tissue samples for 34 species of Ctenotus, including multiple representatives of the 12 recognized ‘species groups’ (Storr et al. 1999) and most of the phenotypic diversity present in this high-diversity genus (Greer 1989; Cogger 2000). A summary of major morphotypes within Ctenotus is provided as electronic supplementary material (electronic supplementary material, table s1). Genomic DNA was extracted using DNeasy Tissue kits (Qiagen). We sequenced 718 bp of mitochondrial ND4 and 181 bp of three flanking tRNAs (tRNA-ser, tRNA-his, tRNA-leu) using primers and protocols described by Reeder (2003). In addition, we sequenced 671 bp of the nuclear ATP synthetase β-subunit intron, after Skinner (in press). We aligned new DNA sequences for Ctenotus to those used in a previous analysis of sphenomorphine skink relationships (Reeder 2003; Skinner in press). The full data matrix contained 81 species, including representatives of all Australian sphenomorphine genera, and 6 outgroup taxa. In addition to ND4, tRNA and ATP synthetase intron sequences, we included 1442 bp of mitochondrial 12S and 16S rRNA sequences for 48 of these species that were available from previous studies (Reeder 2003; Skinner in press). Alignment of ND4 did not require gap insertion and was straightforward. Ribosomal 12S and 16S alignments used secondary structure models from the European ribosomal database (Wuyts et al. 2001), tRNAs were aligned after Macey & Verma (1997) and ATP synthetase intron sequences were aligned using ClustalX (Thompson et al. 1997) with default parameters. GenBank accession numbers, museum voucher numbers and sampling localities for all sequences included in this study are available as electronic supplementary material (electronic supplementary material, table s2).
(b) Phylogenetic analysis
Bayesian phylogenetic analyses were performed on the combined dataset using separate partitions for intron (ATP synthetase), protein coding (ND4) and structural (tRNAs, rRNAs). Appropriate models of molecular evolution for the three partitions were determined using MrModeltest v. 2.2 (Nylander 2004). Recognition of additional data partitions (Brandley et al. 2005) did not change results, and we present analyses based on only these three partitions. We performed two simultaneous runs of Metropolis-coupled Markov chain Monte Carlo (MCMCMC) sampling using the parallel processor version of MrBayes v. 3.1.2 (Altekar et al. 2004), with four incrementally heated chains run simultaneously for 10 million generations and trees sampled every 1000th generation. We discarded the first 7.5 million generations as burn-in, and checked for convergence by (i) testing for stationarity of log L values, (ii) by plotting the posterior probabilities for nodes as a function of the number of generations, and (iii) examining standard deviations of split frequencies for independent runs (Huelsenbeck et al. 2001). We performed two additional runs using identical conditions to verify that results converged on the same posterior distribution of trees.
(c) Variation in diversification rates
We obtained an ultrametric tree for diversification analyses by applying penalized likelihood (PL; Sanderson 2002) to the Bayesian consensus phylogram using the software package r8s (Sanderson 2003). Prior to rate smoothing, we removed all non-Australian sphenomorphine taxa from the tree and applied cross-validation procedure implemented in r8s with an additive penalty to select an optimal smoothing parameter. The diversification rate analyses discussed below are robust to absolute age estimates and require only relative divergence times. However, to facilitate comparison with other Australian radiations for which approximate time scales are available (e.g. Jennings et al. 2003; Crisp et al. 2004), we inferred the age of the basal split among Australian sphenomorphine species. No fossil calibration points exist within either the sphenomorphine radiation or the closely related outgroup taxa (Martin et al. 2004), so we used published rates of squamate mtDNA sequence evolution (Zamudio & Greene 1997; Calsbeek et al. 2003; Richmond & Jokusch 2007). We determined mean pairwise maximum-likelihood (ML) distances between taxa in the two basal clades under the best-fit model and parameters inferred by MrModeltest for this data partition and estimated the timing of the split under an mtDNA divergence rate of 0.895% corrected sequence divergence per Myr (range: 0.47–1.32%; Zamudio & Greene 1997).
Diversification rate analyses were conducted on the genus-level tree created by pruning all but one lineage per genus from the PL tree; species totals were assigned to each terminal following Wilson & Swan (2003). We tested whether some groups of Australian skinks are characterized by exceptionally high or low diversification rates using the approach of Magallon & Sanderson (2001). We plotted the number of species within Australian sphenomorphine genera as a function of stem clade age, as inferred from the PL chronogram. The stem clade age is simply the age of the divergence between the crown group and its extant sister group, corresponding in this case to genus age; this assumes that genera are monophyletic. We do not have complete sampling of all sphenomorphine species, therefore we cannot distinguish between crown- and stem-group ages, and these divergence times thus represent the transition time from a single ancestral species to a clade with some level of standing diversity in the present.
We then considered the Australian sphenomorphine radiation as a whole and estimated the net diversification rate r or the difference between the speciation rate (λ) and the extinction rate (μ). We compared the extant diversity of sphenomorphine genera with the 95% CIs around expected species diversity for clades of similar age that have diversified with the net diversification rate observed for the overall radiation (Magallon & Sanderson 2001, eqn (10)).
Net diversification rates were estimated for the entire sphenomorphine radiation using two methods. First, we used the whole-clade estimator described by Magallon & Sanderson (2001, eqn (7)). This estimator is a function of three parameters: the number of species in the clade; the age of the clade; and the ratio of the extinction rate to the speciation rate μ/λ (denoted by a). Even with a constant net diversification rate, the expected outcome of a stochastic diversification process can vary with the extinction fraction a (e.g. Raup 1985). We therefore computed CIs around expected diversity under a=0 and 0.99; these values represent extremes on a continuum of relative extinction rates. We computed one-tailed probabilities of the observed diversity of each clade, given the estimated overall diversification rate and clade age, to test whether sphenomorphine skinks are characterized by an excess of species-rich or species-poor lineages; we combined individual p-values using the Z-transform method (Whitlock 2005; additional information in the electronic supplementary material, table s4).
We repeated these analyses using an estimator of r that uses both taxonomic and phylogenetic data, because our data consist of terminal taxa of known diversity (e.g. numbers of species within genera), as well as the phylogenetic relationships of these taxa. Our estimator generalizes the Yule process estimator of Sanderson & Wojciechowski (1996) to the birth–death process. The log likelihood of the observed diversity of T terminal taxa is given by(2.1)where(2.2)ni is the diversity of the ith terminal; ti is the stem age of the terminal; and a is the extinction fraction. The log likelihood that each of the N internal branches in the phylogeny will result in exactly two descendant species after some time t is given by(2.3)where xi is the birth time of the ith lineage in units of time before present. The log likelihood of the combined taxonomic and phylogenetic data is simply log L=log LP+log LT. ML estimates of r were found using the one-dimensional ‘optimize’ routine in R (http://cran.r-project.org/). A derivation of this estimator and its relationship to previous taxonomic/phylogenetic estimators is given in appendix A.
(d) Shifts in diversification rate
We tested for shifts in diversification rate within the Australian sphenomorphine radiation by contrasting the likelihood of the data under a model with equal diversification rates for all lineages to the likelihood under a model where an ancestral diversification rate r1 shifts to a new rate r2 along some branch in the tree (Sanderson 1994; Sanderson & Wojciechowski 1996). To compute likelihoods and estimate r, we used the birth–death estimator based on phylogenetic and taxonomic data described in this study. For this model, referred as the flexible-rate model, we sequentially split the tree at each branch and optimized r onto the resulting pair of subtrees. The node resulting in the maximum combined likelihood for the bipartite tree was the ML estimate of the shift point. Analyses were conducted under extinction fractions of a=0 and 0.99.
We predicted that the clade containing Ctenotus and Lerista would show an increase in diversification rates relative to other lineages simply owing to the greater diversity of these genera relative to other sphenomorphine genera. However, an alternative hypothesis is that this clade retained an ancestral but elevated diversification rate, while another clade or clades exhibited a decline in diversification. To distinguish between these possibilities, we repeated our flexible-rate analysis with the constraint that the subtree partition containing Ctenotus must include the root node. Thus, no rate shift was permitted along the path between the root node and Ctenotus. If the likelihood of the data under this constrained model is similar to that under the flexible-rate model, we cannot distinguish between an increase in diversification rates for Ctenotus and a decrease in diversification rates for some other clade. To avoid conditioning results on any particular tree topology and branch lengths, we sampled 1000 phylogenies from the set of post-burn-in topologies (1 tree every 5000 generations) and applied PL to each. We repeated our analyses on this set of trees to generate the posterior distribution of likelihood differences between these two variable-rate diversification models, with the expectation that the rate-decrease model would consistently provide a poorer fit to the data than the flexible-rate model. Functions for all diversification analyses described in this study have been included in the LASER library for the R programming language (Rabosky 2006a).
3. Results and discussion
(a) Phylogenetic relationships
Bayesian phylogenetic analysis supports Ctenotus as monophyletic with a high level of confidence (figure 1), implying that the exceptional species diversity within Ctenotus cannot be explained as a taxonomic artefact. In addition, our results strongly support sister clade relationship between Ctenotus and Lerista, the other highly diverse genus of sphenomorphine skink (79 species). Although taxon sampling of Lerista is limited, morphological and molecular data (Greer 1989; A. Skinner 2007, personal communication) strongly suggest that the genus is monophyletic.
Our results corroborate those of several previous studies (Reeder 2003; Skinner in press) and provide further support for the monophyly of Australian sphenomorphine skinks, as well as a basal split between Notoscincus and all other Australian sphenomorphine lineages. These results were also supported by parsimony and likelihood bootstrap analyses (electronic supplementary material, table s3). As found previously (Reeder 2003), several sphenomorphine genera are not monophyletic (e.g. Eulamprus, Glaphyromorphus). Taken together, monophyly of Ctenotus and Ctenotus plus Lerista suggests higher diversification rates for these lineages relative to the remainder of the Australian sphenomorphine radiation, because 75% of the species in the radiation fall within these genera.
(b) Variation in diversification rates
Cross-validation analysis of PL trees with different smoothing parameters and of a tree with branch lengths estimated under a molecular clock constraint indicated optimal performance of PL with a smoothing parameter of 100. We pruned the PL tree to include 23 lineages for which we could assign levels of taxonomic diversity (figure 2). We inferred a basal divergence time of 28.2 Myr BP (Myr ago; range: 19.1–53.6 Ma) using ML distances between Notoscincus and other Australian sphenomorphine skinks. In the absence of robust fossil calibration points, we prefer to treat this age estimate provisionally, but note that a recently published substitution rate for squamate reptile ND4 (1.6% Myr−1) would imply that the radiation has occurred even more recently (Feldman & Spicer 2006). With the exception of the absolute magnitude of inferred diversification rates, all results are independent of the proposed age of the sphenomorphine radiation.
We were able to assign taxonomic diversity levels to the four paraphyletic genera, because taxon sampling included most or all species within three of those genera (2/2 Coeranoscincus, 2/3 Ophioscincus, 11/13 Glaphyromorphus) and because major lineages within Eulamprus correspond to phenotypically and ecologically distinctive categories (O'Connor & Moritz 2003). To verify that results were robust to alternative species richness assignments, we generated 100 datasets by randomly partitioning generic diversity into paraphyletic lineages (e.g. the 15 recognized species of Eulamprus were divided at random among Eulamprus lineages A, B, C, D and E in each replicate dataset) and repeated the full complement of diversification analyses on each. Shuffling diversity levels within paraphyletic genera in this fashion did not change results described below (electronic supplementary material, figure s1).
Comparisons of observed species diversity within terminal taxa relative to those expected under a homogeneous net diversification rate across all Australian sphenomorphine skinks reject the null hypothesis that net diversification rates have been constant among lineages (p<0.013, a=0; p<0.031, a=0.99; figure 3). Regardless of whether sphenomorphine diversification rates are estimated using the combined taxonomic/phylogenetic approach described in this study (figure 3) or the whole-clade estimator (p<0.001; electronic supplementary material, table s4), there is clearly an excess of species-poor or species-rich lineages. Although this test is conservative (see electronic supplementary material, table s4e), these results are robust to assumptions about the underlying model of extinction. Therefore, clade age alone does not explain the striking disparity in species richness among sphenomorphine skinks.
(c) Shifts in diversification rate
The data clearly reject the constant diversification model in favour of the flexible-rate model (p<0.001; table 1). The ML estimate of the shift point in the flexible-rate model is the node corresponding to the most recent common ancestor (MRCA) of Ctenotus and Lerista (figure 2) under extinction fractions of both a=0 and 0.99; ML estimates of r for Ctenotus/Lerista and all other Australian sphenomorphine lineages suggest that the net diversification rate within this clade has increased approximately 2.5-fold (a=0) to 15-fold (a=0.99) relative to that observed across the remainder of the tree. While palaeontological evidence supports the view that relative extinction rates have generally been high (e.g. Stanley 1979; Gilinsky 1994; Alroy 1996), we are aware of no evidence suggesting that extant clades have diversified in the absence of extinction. It is thus probable that the true magnitude of the rate increase for Ctenotus and Lerista exceeded that inferred under a=0. Rate shifts at other internal nodes are far less likely than those at the MRCA of Ctenotus/Lerista, as assessed by the difference in likelihood scores (ΔL) between the best-fit location of the rate shift and alternative nodes (figure 2): MRCACtenotus–Hemiergis, ΔL=4.06; MRCACtenotus–Anomalopus, ΔL=9.99; shift on Ctenotus branch only, ΔL=7.05; and shift on Lerista branch only, ΔL=8.36 (ΔL given for a=0, but results for a=0.99 are qualitatively similar).
The hypothesis that Lerista/Ctenotus have merely retained high diversification rates from the ancestral sphenomorphine lineage is not supported, as this scenario would require that diversification rates have decreased multiple times throughout the tree. The model specifying increased diversification at the Lerista/Ctenotus MRCA performs much better than a model where no rate increase is permitted between the root node and Lerista/Ctenotus (rate-decrease model; table 1), as assessed by the AIC. This result is not conditional on the topology and branch lengths shown in figure 2, because the posterior distribution of the difference in likelihood scores between the flexible-rate and rate-decrease models strongly favours the former model (electronic supplementary material, figure s2).
Because our analyses are partially dependent on taxonomy, our results could have been influenced by a failure to adequately account for true species diversity within genera. Very few intraspecific phylogeographic studies have been conducted on sphenomorphine skinks, and we predict that future investigations will reveal morphologically cryptic species within many taxa. However, this issue would only influence our results if Lerista and Ctenotus harbour proportionately fewer cryptic species than other taxa.
(d) Aridification of Australia and diversification rates
The progressive expansion of the Australian arid zone during the past 25 Myr appears to have catalysed increased diversification in a number of sclerophyllous plant clades, including eucalypts and casuarinas (Crisp et al. 2004). Likewise, some lineages of agamid lizards experienced rapid diversification during the aridification of Australia (Harmon et al. 2003; Hugall & Lee 2004; Rabosky 2006b). The majority of species within both Lerista and Ctenotus occur in the arid and semi-arid regions of Australia, suggesting a possible link between aridification and diversification rates in this group. Temporal calibration of the sphenomorphine tree (figure 2) yields dates for the radiations of Ctenotus and Lerista that accord well with the known chronology for the aridification of Australia and that show broad congruence with the timing of diversification in these other arid-adapted groups.
Given the diversity of these two genera in the Australian arid zone, it is striking that only a few other lineages of sphenomorphine skinks (figure 2) can be considered arid zone taxa—two species each of Notoscincus and Eremiascincus. Thus, while nearly half of all the species of sphenomorphine skinks occur in the arid and semi-arid regions of the continent, this diversity is concentrated within just 4 of 23 major lineages (figure 2). Yet the arid and semi-arid climatic regions of Australia are defined by a climatic zone that occupies more than three-fourths of the continent's surface area (James & Shine 2000). That so few major lineages have succeeded in entering a geographically vast region, which nonetheless comprises a small fraction of the total climatic diversity in Australia (James & Shine 2000), suggests the possibility that phylogenetic conservatism of climatic tolerances (Ricklefs & Latham 1992; Wiens & Graham 2005) may underlie the exceptional diversification of a few sphenomorphine lineages. Such phylogenetic conservatism of traits is hypothesized to limit the dispersal of lineages between tropical and temperate regions (Ricklefs & Latham 1992; Wiens et al. 2006) and among elevational zones in the montane tropics (Wiens et al. 2007); we similarly hypothesize that such constraints have influenced the ability of lineages to shift from mesic-to-arid environments.
Can these diversification patterns be explained by niche conservatism in conjunction with the larger area of the arid zone, relative to the more mesic regions of continental Australia? There is considerable evidence that the geographical area occupied by a clade exerts a strong effect on diversification rates (Losos & Schluter 2000; Ricklefs 2003; Davies et al. 2005). The mesic-to-arid gradient might act as an environmental filter, limiting the dispersal of lineages between those environments, and the larger area of the arid zone would thus provide an expanded theatre of diversification for clades that successfully made this transition. If arid Australia harbours a diversity of sphenomorphine species proportional to its geographical area, and if a restricted number of clades have been able to enter this climatic zone, we would expect to observe increased diversification rates in those clades as a consequence of the interaction between the geographical area and phylogenetic conservatism of climatic tolerance. Indeed, Ctenotus species richness within major climatic regions of Australia is proportional to the geographical area occupied by those zones (James & Shine 2000). This suggests that part of the answer to ‘why do so many Ctenotus species occur in the arid zone’ (Pianka 1972; Morton & James 1988; James & Shine 2000) may come from understanding why other sphenomorphine clades fail to show this pattern.
The phylogenetic conservatism/geographical area hypothesis might only be part of the explanation for the increased diversification of Ctenotus and Lerista relative to other sphenomorphine skinks. Assuming that the ancestral Australian sphenomorphine was not arid adapted, which is likely based on the general climatic tolerances of sphenomorphine outgroup taxa, only four major lineages have successfully made the transition to arid environments (if the putative shift in climatic tolerance occurred in the ancestor of Ctenotus and Lerista, then only three lineages have made this transition). Yet two of these lineages failed to radiate (Notoscincus and Eremiascincus). Previous studies have found that the Ctenotus lineage is characterized by a substantial increase in both critical thermal maximum temperatures and preferred active temperatures relative to other sphenomorphine species, including Eremiascincus (Benett & John-Alder 1986; Huey & Bennett 1987; Garland et al. 1991). Although data are needed for Lerista and Notoscincus, these results suggest the possibility that traits related to thermal physiology might underlie the dramatic radiation of these groups in the arid zone. However, despite a close phylogenetic relationship, Ctenotus and Lerista have followed radically different evolutionary paths: Ctenotus is diurnal, surface active and shows minimal interspecific variation in body shape, whereas Lerista is fossorial to cryptozoic, frequently nocturnal and shows a tremendous range of interspecific variation in limb reduction and body elongation (Greer 1989).
This study demonstrates exceptional heterogeneity in diversification rates within a major continental radiation. The high species richness of Ctenotus and Lerista relative to other genera of Australian skinks cannot be explained by clade age or as a taxonomic artefact.
Tremendous variation in diversification rates appears to have characterized several continental plant radiations (Hodges & Arnold 1995; Klak et al. 2004; Hughes & Eastwood 2006), but we are unaware of a comparable, sustained contrast in diversification rates within animal radiations restricted to a geographically contiguous region. The diversification of Australian sphenomorphine skinks has occurred against a background of apparently strong phylogenetic niche conservatism, and we speculate that a key physiological innovation may have catalysed diversification in some, but not all, arid-adapted lineages.
This research was supported by NSF-OSIE-0612855, the Mario Einaudi Center at Cornell University, and an EAPSI Fellowship to D. Rabosky administered jointly by NSF and the Australian Academy of Science. All research was conducted in accordance with the legal requirements of the USA and Australia. We thank A. Agrawal, M. Hutchinson, A. McCune and the CLO Evolutionary Biology Lab discussion group for comments on the manuscript and/or valuable discussion of these topics. We are especially grateful to A. Skinner for providing sequences and alignments from Skinner (in press).