Rapid lineage accumulation in a non-adaptive radiation: phylogenetic analysis of diversification rates in eastern North American woodland salamanders (Plethodontidae: Plethodon)

Kenneth H Kozak, David W Weisrock, Allan Larson

Abstract

Adaptive radiations have served as model systems for quantifying the build-up of species richness. Few studies have quantified the tempo of diversification in species-rich clades that contain negligible adaptive disparity, making the macroevolutionary consequences of different modes of evolutionary radiation difficult to assess. We use mitochondrial-DNA sequence data and recently developed phylogenetic methodologies to explore the tempo of diversification of eastern North American Plethodon, a species-rich clade of woodland salamanders exhibiting only limited phenotypic disparity. Lineage-through-time analysis reveals a high rate of lineage accumulation, 0.8 species per million years, occurring 11–8 million years ago in the P. glutinosus species group, followed by decreasing rates. This high rate of lineage accumulation is exceptional, comparable to the most rapid of adaptive radiations. In contrast to classic models of adaptive radiation where ecological niche divergence is linked to the origin of species, we propose that phylogenetic niche conservatism contributes to the rapid accumulation of P. glutinosus-group lineages by promoting vicariant isolation and multiplication of species across a spatially and temporally fluctuating environment. These closely related and ecologically similar lineages persist through long-periods of evolutionary time and form strong barriers to the geographic spread of their neighbours, producing a subsequent decline in lineage accumulation. Rapid diversification among lineages exhibiting long-term maintenance of their bioclimatic niche requirements is an under-appreciated phenomenon driving the build-up of species richness.

Keywords:

1. Introduction

Molecular phylogenies and macroevolutionary models permit testing hypotheses of the tempo and mode of species diversification with increasing statistical rigor (Nee et al. 1992; Paradis 1997; Pybus & Harvey 2000; Harmon et al. 2003). Not surprisingly, adaptive radiations, clades that exhibit rapid divergence into an exceptional variety of adaptively disparate forms (Schluter 2000; Losos & Miles 2002), have served as model systems for quantifying rates of diversification. However, few studies have quantified the tempo of diversification in species-rich clades that contain negligible adaptive disparity (non-adaptive radiations; Gittenberger 1991). The importance of obtaining diversification-rate estimates for these systems is highlighted by recent arguments that allopatric speciation, considered the predominant mode of animal speciation (Barraclough & Vogler 2000; Coyne & Orr 2004), may occur rapidly in lineages multiplying through apparently non-adaptive, geographic processes (Near & Benard 2004). Without such empirical data, the macroevolutionary consequences of different modes of evolutionary radiation are difficult to assess.

Woodland salamanders of the genus Plethodon provide an excellent opportunity to study rates of diversification in an evolutionary radiation exhibiting little apparent adaptive phenotypic disparity. These lungless salamanders lack an aquatic larval stage, require moist terrestrial habitats for survival and reproduction, and occupy forested areas across portions of eastern and western North America (Highton 1995). They achieve their greatest diversity in eastern North America, where 46 of the 55 recognized species occur (Highton et al. 1989; Highton 1997, 1999, 2004; Lazell 1998; Highton & Peabody 2000; Mead et al. 2005). Many species are endemic to the Appalachian Highlands, where they are among the most abundant and ecologically important vertebrates in forest ecosystems (Davic & Welsh 2004). Despite ancient histories of evolutionary separation (Highton & Larson 1979; Maxson et al. 1979; Mahoney 2001), species of Plethodon exhibit extreme morphological stasis (Wake et al. 1983; Larson 1984; Carr 1996; Adams & Rohlf 2000), suggesting low rates of adaptive phenotypic evolution relative to other plethodontid taxa (García-París et al. 2000; Parra-Olea et al. 2004; Kozak et al. 2005). Geographic patterns of genetic variation are described for hundreds of natural populations and the geographic limits of many morphologically cryptic species are well known (Highton 1995). In addition, molecular-clock estimates of divergence times between major Plethodon lineages provide a foundation for quantifying the absolute time course of diversification (Highton & Larson 1979; Maxson et al. 1979; Hass et al. 1992; Larson et al. 2003).

This study provides a comprehensive assessment of the tempo of lineage diversification in eastern North American Plethodon using an mtDNA phylogeny. This analysis goes beyond earlier phylogenetic studies in using representatives of nearly all of the 46 extant species of eastern Plethodon, including intraspecific geographic sampling for many species. We use a recently developed hypothesis-testing framework to investigate whether rates of lineage accumulation (speciation minus extinction) vary among major clades of eastern Plethodon and throughout the phylogenetic history of the group. Our study demonstrates that peak rates of diversification in groups exhibiting extreme phenotypic stasis can be comparable to those reported for classic examples of adaptive radiation.

2. Material and methods

(a) Taxon sampling

In phylogenetic studies of diversification, it is important to use a sampling regime that adequately reflects the extant variation in the study system. Knowledge of species limits in Plethodon is remarkable in that geographic limits of many morphologically cryptic species are well known as a result of comprehensive population-level studies of geographic genetic variation (Highton et al. 1989; Highton 1995, 1997, 1999; Highton & Peabody 2000). This study includes 44 of the 46 recognized species of eastern Plethodon (excluding the recently described P. ainsworthi; Lazell 1998 and P. sherando; Highton 2004), many of which have geographically restricted ranges in the southern Appalachian Highlands and do not harbour additional cryptic, independent evolutionary lineages (Weisrock and Larson, in press; Weisrock et al. 2005). For species with more widespread distributions, we included multiple geographic samples to prevent exclusion of phylogenetically distinct population-level lineages not detected in earlier studies. We include only two samples of the widespread species, Plethodon cinereus; however, a recent mtDNA study using greater geographic sampling indicates that intraspecific genetic divergence is not as great as divergence between species of the P. cinereus group (Sites et al. 2004). Two species of western North American Plethodon, two additional plethodonines, a bolitoglossine, and a spelerpine were used as out-group taxa. Geographic sampling locations for all samples are shown in electronic supplementary material.

DNA extraction, PCR, and sequencing were performed as in Weisrock et al. (2005). Amplification and sequencing of approximately 1100 nucleotide bases of mtDNA sequence spanning the complete ND2 and tRNATrp genes employed primers L4437 (Macey et al. 1997), L5195 (Weisrock et al. 2005), and a new reverse primer anchored in tRNAAla (5′-AAAGTGTTTGAGTTGCATTCA-3′). Sequence alignment was straightforward and unambiguous.

(b) Phylogenetic reconstruction

Maximum-parsimony trees were constructed in PAUP* 4.0b10 (Swofford 2002) using a heuristic search with 100 random-addition replicates and TBR branch swapping. To assess support for individual nodes, 100 non-parametric bootstrap replicates were performed under similar conditions, except that 10 random additions were performed per bootstrap replicate.

The appropriate model of nucleotide substitution was determined through likelihood-ratio tests implemented in MrModeltest v.2.2 (written by J. A. A. Nylander: http://www.ebc.uu.se/systzoo/staff/nylander.html). Bayesian phylogenetic analysis was performed using the parallel-processor version of MrBayes v.3.04 (Altekar et al. 2004). Five Markov chains were used with the temperature profile at the default setting of 0.2. Random starting trees and flat priors for all parameters were used to begin each Markov chain. Two million generations were run with tree samplings taken every 1000 generations. Sampled trees were parsed with MrBayes to construct a phylogram based upon mean branch lengths and to calculate posterior probabilities of all branches using a majority-rule consensus. All saved trees from generations prior to a stationary ln-likelihood (lnL) value were discarded. To account for the possibility that an analysis may not approach the optimal posterior distribution, three additional runs were performed using identical conditions.

(c) Estimates of diversification rates

To obtain ultrametric trees for diversification analyses, the Bayesian consensus phylogram was subjected to rate smoothing by penalized likelihood (PL; Sanderson 2002a) using the software package r8s (Sanderson 2002b). All out-group taxa were pruned prior to rate smoothing. PL requires fixation of the root node to estimate dates of divergences. We set the divergence time between the P. cinereus group and the remaining species of eastern Plethodon to 27 million years, a calibration point derived from allozymic and immunology-based studies of divergence times (Highton & Larson 1979; Maxson et al. 1979). An optimal smoothing value was obtained using the cross-validation procedure implemented in r8s.

To visualize the tempo of lineage accumulation (speciation minus extinction), we plotted the natural logarithm of the number of lineages against the branch-length distance from the root node on the rate-smoothed Bayesian consensus phylogram. We tested for a significant departure from the null hypothesis of a constant rate of diversification using the constant-rate (CR) test (Pybus & Harvey 2000; Pybus & Rambaut 2002), which uses the γ-statistic to compare the relative positions of nodes in the phylogeny to those expected under a CR model of diversification. Negative values of γ indicate that nodes are closer to the root than expected under the null model and imply a deceleration in the accumulation of lineages. Positive values of γ indicate that nodes are closer to the tips and imply an acceleration of the accumulation of lineages. A CR model of diversification can be rejected at the 95% level of significance if γ<−1.645 (Pybus & Harvey 2000).

We also analysed the distribution of relative divergence times among species using survivorship models described by Paradis (1997). Diversi v.0.20.0 (Paradis 2000) was used to test three alternative models of lineage accumulation in a maximum-likelihood framework. Model A specifies a constant rate of lineage accumulation (δ) and represents the null hypothesis of no heterogeneity in rates through time. Model B specifies a gradual change in the rate of lineage accumulation through time by estimating an additional parameter, β. Values of β>1 indicate that δ has decreased through time; values of β<1 indicate that δ has increased through time. Model C specifies two different rates of lineage accumulation before (δ1) and after (δ2) a defined breakpoint in time, TC. Because model A is a special case of models B and C, a likelihood-ratio test (LRT) is appropriate to determine if A is favoured over B or C. For comparing models B and C, which are not nested, the Akaike Information Criterion (AIC; Akaike 1973) is used to determine the favoured model. The model with the lowest AIC is the one that best describes temporal variation in relative divergence times.

An interpretation that the CR test and survivorship analyses reflect instantaneous lineage-accumulation rates at past times assumes that speciation and extinction occur with equal probability among lineages during a given time interval and that all extant species have been sampled. To test the former assumption, we used the relative-cladogenesis statistic, which uses a broken-stick distribution to identify branches ancestral to a greater proportion of the extant species than expected by chance (Nee et al. 1992; Rambaut et al. 1997). This test calculates the probability (Pk) that a particular lineage at time t will have k tips given the total number of tips at time 0 (the present).

Incomplete taxon sampling at the species level tends to underestimate the number of nodes toward the present and may cause an apparent slowdown in diversification rate (Pybus & Harvey 2000). To explore the potential effects of unidentified missing lineages, we used Monte Carlo simulations to generate γ-statistic values under the assumption that the inferred number of lineages (n) in a clade represents only a proportion ( f) of the actual clade total (Pybus et al. 2002). This procedure permits an assessment of whether the γ-statistic for the actual dataset could be an artefact of incomplete sampling. For values of f ranging from 0.05 to 1.0, we simulated 1000 pure-birth phylogenies with n/f tips using PhyloGen (Rambaut 2002) and then calculated their corresponding γ-values with Genie (Pybus & Rambaut 2002). For each of the simulated datasets, the mean γ-value and its 0.975 and 0.025 percentiles were plotted to estimate the 95% confidence interval for the number of lineages that must be missing from the sample to get a γ-statistic as extreme as the measured one if γ actually is zero.

3. Results

(a) Phylogenetic relationships

Eighty sequences representing 1108 aligned bases are reported for 44 species of Plethodon and four out-group taxa. Absence of premature stop codons in the ND2 protein-coding region, functional stability of the tRNATrp gene, and strong bias against guanine on the light strand indicate that we have collected authentic mitochondrial DNA sequences. Including the out-group taxa, the dataset contains 758 variable sites of which 651 are parsimony informative (584 within eastern Plethodon). Likelihood-ratio tests favour the Gtr+I+Γ model of nucleotide substitution.

Bayesian phylogenetic analysis estimates a posterior distribution with a mean lnL of −18 700.09 (s.d.=9.91) following a ‘burn in’ of 100 000 generations (figure 1). Parsimony analysis produces a single most-parsimonious tree of 3842 steps. Because both analyses produce highly congruent estimates of phylogenetic relationship among the major haplotype clades, only the Bayesian consensus phylogram is presented with bootstrap values from the parsimony analyses included for shared branches.

Figure 1

Bayesian 50% majority-rule consensus phylogram. Posterior probabilities greater than or equal to 0.95 are denoted with an asterisk above the branches and non-parametric bootstrap proportions for the parsimony analysis below. Branch lengths represent means from the posterior distribution of trees.

Phylogenetic relationships among major lineages of Plethodon are congruent with those of earlier studies based on nuclear-encoded allozymes (Highton & Larson 1979), albumins (Hass et al. 1992) and a different mtDNA region (Mahoney 2001). Eastern North American Plethodon form a strongly supported clade. In addition, the P. cinereus, P. wehrlei and P. glutinosus groups are each resolved as strongly supported monophyletic groups. Species of the P. welleri group are resolved as a strongly supported clade, except for P. websteri, which is deeply divergent from any of the recognized species groups. With the exception of the sister-group relationship between the P. cinereus group and a clade containing all other species of eastern Plethodon, phylogenetic relationships among the species groups are not strongly supported.

(b) Diversification rates

The PL tree for 48 putatively independent population-level lineages of eastern Plethodon is shown in figure 2. Survivorship analyses reject a CR model of lineage accumulation in favour of declining rates of diversification (table 1). The likelihood value for model B, which specifies a gradual decline in rate of lineage accumulation, is nearly indistinguishable from model C, which assumes an abrupt change in the diversification rate at 2 million years before present (MYBP). The CR test does not reject a constant-rates model of diversification. However, the negative γ-statistic is consistent with the survivorship analyses, indicating that the origins of extant lineages are clustered disproportionately early in the history of the group.

Figure 2

Ultrametric Bayesian consensus phylogram for 48 lineages of eastern North American Plethodon. Branch lengths are proportional to absolute time based on a calibration point of 27 MYBP for the phylogenetic split between the P. cinereus group and all remaining species of eastern Plethodon (Highton & Larson 1979; Maxson et al. 1979). Dashed branches denote lineages that experienced significantly higher rates of diversification according to the relative-cladogenesis statistic (Nee et al. 1992; Rambaut et al. 1997).

View this table:
Table 1

Tests of the null hypothesis that the rate of lineage accumulation (speciation minus extinction) was constant during the phylogenetic history of eastern North American Plethodon. (Model A specifies a constant lineage-accumulation rate. Model B specifies a gradual decrease in the lineage-accumulation rate. Model C specifies different lineage-accumulation rates before and after a breakpoint in time, TC, which is the value of TC that maximized the ln likelihood for Model C. The three models were fit to the entire eastern Plethodon clade, and to the P. glutinosus group separately because the latter clade diversified at a significantly higher rate than the other lineages of eastern Plethodon according to the relative-cladogenesis statistic (Nee et al. 1992; Rambaut et al. 1997).)

The relative-cladogenesis statistic indicates significant differences in diversification rates among lineages within eastern Plethodon. In particular, the P. glutinosus group has experienced a higher rate of diversification relative to all other nodes in the phylogeny (Pk=0.004; figure 2). The lineages-through-time (LTT) plot for this clade is shown in figure 3. Both the CR test and survivorship analyses strongly reject the hypothesis that rates of lineage accumulation have remained constant over time in the P. glutinosus group (table 1). The AIC favours model C with an abrupt shift from a higher to lower diversification rate at 8 MYBP. Model C estimates a lineage-accumulation rate of 0.80 species per million years (s myr−1) (95% CI: 0.22–1.38 s myr−1) between 11.3 and 8 MYBP and 0.16 s myr−1 between 8 and 0 MYBP (95% CI: 0.09–0.21 s myr−1). Monte Carlo simulations of γ-statistic values indicate that 134 P. glutinosus-group lineages (95% CI: 15–670) would need to be added to measure a γ-statistic of −2.30 when the actual value is not significantly different from zero.

Figure 3

Semi-logarithmic plot of lineages-through-time for the P. glutinosus group based on the calibration shown in this figure.

4. Discussion

Major clades of eastern North American woodland salamanders exhibit significant differences in rates of lineage accumulation. In particular, the relative-cladogenesis statistic indicates that the ancestral lineage leading to the P. glutinosus group has generated a greater proportion of the extant species diversity of eastern Plethodon than expected by chance (figure 2). Furthermore, the tempo of lineage accumulation within this exceptionally diverse group has not been constant throughout its phylogenetic history. Both the CR test and survivorship analyses of branching times in the mtDNA phylogeny indicate that the P. glutinosus group underwent an early burst of lineage splitting that produced a significantly larger number of surviving lineages than expected under a CR model of diversification (table 1). In comparison to the P. glutinosus group, species in the P. cinereus, P. welleri, and P. wehrlei groups have more northern distributions with large portions of their current distributions in areas unsuited for occupation during Pleistocene glacial maxima (Highton 1995). Although some overlap occurs in the geographic distributions of these major clades, the P. glutinosus species group has a larger and more southern geographic distribution than the other species groups (Highton 1995), which may have provided greater opportunities for speciation and long-term persistence of lineages (the time-for-speciation effect; Stephens & Wiens 2003).

The P. glutinosus group exhibits remarkable rates of lineage accumulation early in its evolutionary history. The estimated rate of lineage accumulation between 11.3 and 8 MYBP is 0.80 s myr−1. This rate is comparable to those documented for examples of adaptive radiation including Hawaiian silverswords (0.56 s myr−1; Baldwin & Sanderson 1998), Neogene horses (0.50–1.4 s myr−1; Hulbert 1993), Lake Tanganyika cichlids (0.75–1.49 s myr−1; McCune 1997), and south African ice plants (0.77–1.75 s myr−1; Klak et al. 2004). Because this is a net diversification rate that incorporates both speciation and extinction, it is likely that rates of lineage formation during this time-period were higher. This high rate of diversification is particularly noteworthy for its temporal depth. Some of the most striking cases of rapid lineage accumulation are more recent (e.g. Darwin's finches, Sato et al. 1999; Lake Malawi cichlids, Albertson et al. 1999), with less time for extinction to lower overall species diversity.

Use of molecular phylogenies to assess the tempo of lineage accumulation relies critically on adequate representation of extant independent evolutionary lineages (Pybus & Harvey 2000; Pybus et al. 2002). Incomplete taxon sampling underestimates the number of lineages toward the present and may produce an artefactual slowdown in lineage diversification. Knowledge of geographic patterns of lineage diversity in eastern Plethodon is unparalleled. Much of the species-level diversity in the group is morphologically cryptic, and most species have been subjected to detailed studies of geographic genetic variation (Highton et al. 1989; Highton 1995, 1997, 1999; Highton & Peabody 2000; Weisrock and Larson, in press; Weisrock et al. 2005). Therefore, it is unlikely that our diversification-rate conclusions are biased by a failure to account for cryptic lineage diversity. Furthermore, our phylogenetic simulations demonstrate that many additional cryptic species would need to be present to make the measured decrease in diversification an analytical artefact.

Two lines of evidence suggest that the rapid early diversification of the P. glutinosus species group was not driven by adaptive differentiation. First, the vast majority of species pairs in the P. glutinosus group have allopatric distributions, and cases of sympatry involve distantly related species; this spatial and temporal pattern indicates that speciation occurred through vicariant geographic processes followed by range shifts (Barraclough & Vogler 2000; Losos & Glor 2003) without obvious adaptive divergence. Second, morphological disparity within all major species groups of eastern Plethodon is minimal and consists primarily of minor variation in coloration and body size (Highton et al. 1989; Highton 1995; Carr 1996; Highton & Peabody 2000). Character displacement occurs in sympatry at the borders of some species' ranges (Adams & Rohlf 2000; Adams 2004), but these differences are not characteristic of the species as a whole.

What then has driven the exceptional rates of diversification that characterize the early history of the P. glutinosus species group? We support recent arguments that the long-term stability of an ecological niche has profound consequences for the build-up of species diversity because it causes species to track their preferred habitats in the face of environmental change, promoting vicariant fragmentation and formation of new evolutionary lineages (e.g. Levinton 2001; Wiens 2004a,b). We hypothesize that the rate of vicariant lineage splitting is accelerated in clades that exhibit long-term maintenance of an ecological niche and whose geographic distributions experience temporal climatic fluctuations producing alternating periods of expansion and contraction of favourable habitat. The evolutionary history of woodland salamanders exhibits both of these characteristics; Plethodon has occupied a forest-floor niche since the Eocene (Maxson et al. 1979; Larson et al. 1981) and is endemic to the North American temperate zone, a region with a dramatic history of climate-driven forest contraction and expansion (Webb & Bartlein 1992; Jansson & Dynesius 2002). Recent declines of some Plethodon populations illustrate their susceptibility to fragmentation by habitat alteration (Highton, in press).

We propose that ecological and genetic interactions (Hairston 1980; Nishikawa 1985; Highton 1995) limit the distributional overlap of species and account for the shift to a lower diversification rate in the P. glutinosus group; once the major geographic lineages were established early in the history of this clade, this pre-emptive occupation of geographic space (Jockusch & Wake 2002) has limited opportunities for subsequent lineage expansion and accumulation. In support of this hypothesis, Plethodon populations can expand rapidly into favourable habitats not occupied by close relatives, as evidenced by the rapid geographic expansion of populations into a vast, formerly glaciated territory within the past 10 000 years (Highton & Webster 1976). Hybridization is common but restricted to the geographic borders of P. glutinosus-group species (Highton & Peabody 2000), thereby limiting their geographic expansion and accumulation of new lineages.

The relationship between niche conservatism and rates of diversification proposed here has received little attention compared to models that link a build-up of species richness to adaptive ecological divergence. Nevertheless, it is potentially a general phenomenon as evidenced by strikingly similar temporal patterns of species accumulation in eastern North American forest-dwelling birds. Extant lineages of Dendroica warblers display limited disparity in traits that characterize other avian adaptive radiations and originate disproportionately early in the evolutionary history of that group near the Miocene–Pliocene boundary (Lovette & Bermingham 1999). An abrupt increase in temperature and aridity spanning this time-period contracted and fragmented North America's forested habitats and increased grassland habitats (Webb et al. 1995), promoting rapid vicariant speciation of co-distributed woodland populations that require mesic, forested habitats for survival and reproduction (Highton 1995; Lovette & Bermingham 1999).

5. Conclusions

We demonstrate that allopatrically distributed lineages exhibiting extreme phenotypic stasis can experience peak rates of lineage formation that match the tempo of diversification in adaptive radiations. Where a lineage's ecological niche is evolutionarily stable, and where its favoured habitat alternates temporally between geographically continuous versus fragmented conditions, lineages accumulate rapidly early in the group's history and form a vicariant geographical pattern that persists through long-periods of evolutionary time. Rapid accumulation of vicariant lineages exhibiting long-term stasis of their bioclimatic niche is an under-appreciated phenomenon in evolutionary biology.

Acknowledgements

We are especially grateful to R. Highton for his gracious contribution of samples for phylogenetic analysis. P. Chippindale and J. Wiens also provided key DNA samples. Funding was provided by Highlands Biological Station's Bruce Family and Charles Ash Scholarships (KHK), a Highlands Biological Station Grant-in-Aid (DWW), American Society of Ichthyologist's and Herpetologists Gaige Award (KHK), and NSF Doctoral Dissertation Improvement Grants (KHK & AL and DWW & AL).

Footnotes

  • Present address: Department of Ecology and Evolution, Stony Brook University, Stony Brook, NY 11794, USA.

  • Present address: Department of Biology, University of Kentucky, 101 Morgan Building, Lexington, KY 40506-0225, USA.

  • These authors contributed equally to the development of this study.

  • The electronic supplementary material is available at http://dx.doi.org/10.1098/rspb.2005.3326 or via http://www.journals.royalsoc.ac.uk.

    • Received July 15, 2005.
    • Accepted September 2, 2005.

References

View Abstract