Abstract

Mesic forests in the North American Pacific Northwest occur in two disjunct areas: along the coastal and Cascade ranges of Oregon, Washington, and British Columbia as well as the Northern Rocky Mountains of Idaho, Montana, and British Columbia. Over 150 species or species complexes have disjunct populations in each area, and a priori hypotheses based on phytogeography and geology potentially explain the disjunction via either dispersal or vicariance. Here, we test these hypotheses in the disjunct salamander complex Plethodon vandykei and P. idahoensis by collecting genetic data (669 bp of Cyt b) from 262 individuals. Maximum likelihood analysis indicated reciprocal monophyly of these species, supporting the ancient vicariance hypothesis, whereas parametric bootstrap and Bayesian hypothesis testing allow rejection of the dispersal hypothesis. The coalescent estimate of the time since population divergence (estimated using MDIV) is 3.75 × 106 years, and the 95% credibility interval of this value overlaps with the geological estimate of vicariance, but not the hypothesized dispersal. These results are congruent with the pattern seen in other mesic forest amphibian lineages and suggest disjunction in amphibians may be a concerted response to a geological/climatological event. Within P. idahoensis, we tested the corollary hypothesis of an inland Pleistocene refugium in the Clearwater drainage with nested clade analysis and coalescent estimates of population growth rate (g). Both analyses support post-Pleistocene expansion from the Clearwater refugium. We corroborated this result by calculating Tajima's D and mismatch distribution within each drainage, showing strong evidence for recent population expansion within most drainages. This work demonstrates the utility of statistical phylogeography and contributes two novel analytical tools: tests of stationarity with respect to topology in the Bayesian estimation, and the use of coalescent simulations to test the significance of the population growth-rate parameter.

Biogeography and geology have been complimentary disciplines for most of the last two centuries, and comparative phylogeography is a recent extension of this tradition to a restricted taxonomic scale. Lineages with disjunct distributions have long attracted attention (e.g., Lyell, 1832), largely because they provide data critical to our understanding of evolutionary processes associated with lineage divergence. However, the study of disjunctions within species or between sister species is impeded by the continuum of evolutionary processes acting to structure genetic variation at these levels (Fig. 1). Researchers have typically approached phylogeographic data with methodological tools developed from either end of this continuum, with either a phylogenetic background (e.g., Sullivan et al., 1997) or a population genetic background (e.g., Gysels et al., 2004), and often limit the methods used to investigate lineage divergence to those consistent with their background. The integration of methods designed to detect signal at various positions of the continuum of genetic variation is critical to understanding the full spectrum of evolutionary and ecological forces that have structured genetic diversity. Furthermore, congruent results will reinforce confidence in their accuracy, whereas incongruent results may suggest ways that the data can be improved or instances where assumptions of a particular method may be violated. Although examples of the power of this integrated approach to phylogeography exist (e.g., Wares and Cunningham, 2001), the majority of phylogeographic studies employ a limited range of analytical methods. The importance of methods that attempt to infer signal from wide portions of this continuum is underscored by recent contention with respect nested clade analysis (Knowles and Maddison, 2003; Templeton, 2004), perhaps the only method that attempts to both account for geography explicitly and to utilize genetic information from across the continuum of genetic variation. Here, we apply a wide range of methods developed to infer signal at different levels of genetic diversity to the investigation of a lineage of salamanders with a disjunct distribution in the Pacific Northwest in order to maximize the inferences possible with the most basic phylogeographic data.

Figure 1

The continuum of genetic variation (top) and analytical methods used to interpret the genetic signal (bottom). The bars on either side of the method name show the approximate range to which we have applied these tests.

Figure 1

The continuum of genetic variation (top) and analytical methods used to interpret the genetic signal (bottom). The bars on either side of the method name show the approximate range to which we have applied these tests.

Mesic forests in the Pacific Northwest of North America (PNW) are characterized by both the late-successional dominance of western hemlock (Tsuga heterophyla) and western redcedar (Thuja plicata), and by mild, wet winters, and cool, relatively dry summers. The PNW can be divided longitudinally into three provinces: the Cascades and Coastal Mountains in the west, the Northern Rocky Mountains in the east, and the intervening xeric Columbia Basin. Mesic forest ecosystems occur as two geographically separate and ecologically isolated entities in the Cascades/Coastal Mountains and in the Northern Rocky Mountains. More than 150 plant, animal, and fungi lineages are known to have disjunct distributions, with populations in both coastal and inland mesic forests of the PNW (Brunsfeld et al., 2001; Neilson et al., 2001). The ecosystems of the PNW have been influenced by both ancient and recent geological events. The Northern Rocky Mountains were formed ca. 45 to 36 Mya during the Eocene and pollen records indicate the occurrence of coniferous forests on the western slopes of the Rocky Mountains since the mid-Eocene (Graham, 1993). The Cascades Range was formed during the early Pliocene (5 to 2 Mya), and this orogeny established a rain shadow that dried out the Columbia Basin, resulting in a well-developed steppe vegetation by the end of the Pliocene (Daubenmire, 1952; Graham, 1993). Subsequently, glaciers have repeatedly formed and retreated, especially the Cordilleran Ice Sheet in the Northern Rockies (Richmond et al., 1965; Waitt and Thorson, 1983; Delcourt and Delcourt, 1993).

This complex history and the broad taxonomic range of disjunct mesic forest taxa has led researchers to propose alternative hypotheses to explain the isolation of inland mesic forest lineages from their coastal counterparts. First, inland elements of the mesic forests could be the result of an ancient vicariance caused by the Pliocene xerification of the Columbia Basin, with subsequent persistence of inland taxa in refugia throughout the Pleistocene glacial cycles (Brunsfeld et al., 2001). This ‘ancient vicariance’ hypothesis predicts deep divergence and reciprocal monophyly between coastal and inland populations because it allows sufficient time for completion of lineage sorting. A corollary to the ancient vicariance hypothesis is the persistence of inland mesic forests in one or more refugia throughout the Pleistocene. The Clearwater Drainage in northcentral Idaho has been proposed as one such refugium, based on both geology and phytogeography (Daubenmire, 1975). This corollary hypothesis predicts a genetic signature of an expansion from the Clearwater Refugium. An alternate hypothesis explains the populations in the Northern Rockies as the result of recent (i.e., post-Pleistocene) dispersal from the Cascades (Mack et al., 1978; Brunsfeld et al., 2001). The ‘recent dispersal’ hypothesis predicts shallow (post-Pleistocene) divergence between northern coastal and inland populations, absence of reciprocal monophyly, and very little genetic diversity within and among inland populations. DNA evidence supporting ancient vicariance has been identified in the Tailed frog (Ascaphus truei/A. montanus) species complex (Nielson et al., 2001; mtDNA), whereas recent dispersal is supported by mtDNA in the water vole Microtus richardsoni (Demboski et al., unpublished; mtDNA), and cpDNA in the whitebark pine Pinus albicaulis (Richardson et al., 2002).

We focus here on the Plethodon vandykei species group, a salamander species complex that exhibits the characteristic disjunction and is most often found in affiliation with seeps or streams in areas with exposed and fractured bedrock (Wilson et al., 1995, 1997). It includes P. vandykei and P. idahoensis, species that were only recently recognized as specifically distinct (Howard et al., 1993). P. vandykei occurs in Washington in three isolated areas: the Olympic peninsula, the Mt. Rainier/Mt. St. Helens region, and Willapa Hills (Nussbaum et al., 1983). P. idahoensis is found in northern Idaho, extreme western Montana, and southern British Columbia (Wilson and Larson, 1998a; Wilson and Ohanjanian, 2002), and is isolated from P. vandykei by 300 km of desert sage/scrub and xeric forests (Fig. 2; Brunsfeld et al., 2001). Available data, including allozymes (Howard et al., 1993) and morphometrics (Wilson and Larsen, 1998b), indicate that each species is diagnosable. However, because each of those studies lacked a phylogenetic framework, explicit tests of alternative explanations for the disjunction have not been conducted.

Figure 2

(A) Map of North America. (B) Collection localities for P. idahoensis and P. vandykei. Locality numbers correspond to those in Appendix 1. The approximate ranges of P. vandykei and P. idahoensis are shown with dotted white lines. The maxima of the Cordilleran Ice Sheet is shown with a black line. States or provinces are abbreviated as follows: Montana (MT), Idaho (ID), Oregon (OR), Washington (WA), and British Columbia (BC). (C) Details of localities in central Idaho. The approximate boundary of the Clearwater River drainage is shown with a dotted line.

Figure 2

(A) Map of North America. (B) Collection localities for P. idahoensis and P. vandykei. Locality numbers correspond to those in Appendix 1. The approximate ranges of P. vandykei and P. idahoensis are shown with dotted white lines. The maxima of the Cordilleran Ice Sheet is shown with a black line. States or provinces are abbreviated as follows: Montana (MT), Idaho (ID), Oregon (OR), Washington (WA), and British Columbia (BC). (C) Details of localities in central Idaho. The approximate boundary of the Clearwater River drainage is shown with a dotted line.

Here, we combine statistical phylogenetics with coalescent simulations to assess the geographic structure of mtDNA variation in the P. vandykei species group. Specifically, we test the dispersal versus vicariance hypotheses, and test the Clearwater refugium hypothesis. The nested nature of the ancient vicariance/Clearwater refugium hypothesis makes the application of a range of methods that span the continuum of genetic variation particularly appropriate here. The results of this study are then interpreted in light of phylogeographic studies of other mesic forest disjuncts in an ongoing synthesis of the ecosystem genetics of mesic forests in the PNW.

Materials and Methods

Sample Collection

We collected tail clips from 243 P. idahoensis from 49 localities and 18 P. vandykei from 3 localities (see Appendix 1 for detail, available at the Society of Systematic Biologists website, http://systematicbiology.org). A single E. eschscholtzii was collected and used as an outgroup (Mahoney, 2001). The P. idahoensis samples are from seven river drainages and one lake basin whereas the P. vandykei samples are from three river drainages (Fig. 2). Sample sizes ranged from 1 to 13 per locality, with an average n = 5. When P. idahoensis samples were pooled by drainage, sample sizes ranged from 9 to 60 (Table 1), with an average n = 30.4. Therefore, following Saunders et al. (1984; P = [n − 1]/[n + 1]), there is a high probability that we have sampled the deepest coalescent event in each drainage (P = 0.8–0.97). However, these estimates should be conservative because ancestral haplotypes are expected to be most frequent in the population (Waterson and Guess, 1977; Crandall and Templeton, 1996).

Table 1.

Diversity statistics by river drainage for P. idahoensis. River drainage, the number of samples collected in that drainage (n), the probability of the samples from the drainage capturing the deepest coalescent event (Prob), the number of haplotypes in the drainage (No. hap), the number of polymorphic sites (s), the nucleotide diversity (π) ± 1 standard deviation, Tajima's D (with significant values in bold), and the P-value for the mismatch distributions are shown

Drainage n Prob No. hap s π D Mismatch 
Clearwatera 89 0.98 29 38 0.0104 ± 0.005 n/a n/a 
 Selway River 15 0.88 0.0015 ± 0.001 −1.89 0.58 
 Lochsa River 23 0.92 10 0.0024 ± 0.002 −1.37 0.35 
 NFCb River 51 0.96 18 23 0.0073 ± 0.004 −0.13 0.65 
St. Joe River 60 0.97 13 28 0.0031 ± 0.002 −2.10 0.80 
Coeur d'Alene River 19 0.9 0.0064 ± 0.004 −1.81 0.69 
Kootenai River 38 0.95 16 37 0.0017 ± 0.001 −1.97 0.71 
Columbia River 28 0.93 11 0.0002 ± 0.000 −1.16 0.73 
Lake Coeur d'Alene 0.8 n/a 
P. idahoensis (total) 243 0.99 63 97 0.0069 ± 0.004 n/a n/a 
Drainage n Prob No. hap s π D Mismatch 
Clearwatera 89 0.98 29 38 0.0104 ± 0.005 n/a n/a 
 Selway River 15 0.88 0.0015 ± 0.001 −1.89 0.58 
 Lochsa River 23 0.92 10 0.0024 ± 0.002 −1.37 0.35 
 NFCb River 51 0.96 18 23 0.0073 ± 0.004 −0.13 0.65 
St. Joe River 60 0.97 13 28 0.0031 ± 0.002 −2.10 0.80 
Coeur d'Alene River 19 0.9 0.0064 ± 0.004 −1.81 0.69 
Kootenai River 38 0.95 16 37 0.0017 ± 0.001 −1.97 0.71 
Columbia River 28 0.93 11 0.0002 ± 0.000 −1.16 0.73 
Lake Coeur d'Alene 0.8 n/a 
P. idahoensis (total) 243 0.99 63 97 0.0069 ± 0.004 n/a n/a 
a

NFC, Lochsa, and Selway are the major tributaries of the Clearwater drainage.

b

NFC = North Fork of the Clearwater.

DNA was extracted from 10 to 20 mg tail clips, which had been stored in 90% EtOH, with either the DNeasy Tissue kit (Qiagen, Inc.; Valencia, CA), following the manufacturer's instructions for rodent tails, or the method described in Nielson et al. (2001). A 669-bp fragment of the mitochondrial cytochrome b gene (Cyt b) was amplified with the forward primer MVZ 15 (5′-GAACTAATGGCCCACACWWTACGNAA-3′; Moritz et al., 1992) and with H-14963 (5′-GGCAAA-TAGAARTATCATT-3′; Sullivan et al., 2000) as the initial reverse primer. For subsequent polymerase chain reaction (PCR) amplification, we designed a lineage-specific reverse primer (5′-AGGAGTGAGAGTAGAGTAAGTA-3′). Amplicons were purified using polyethylene glycol precipitation, and sequencing reactions were performed with the BigDye Kit (Applied Biosystems, Inc., Foster City, CA) with 10 to 20 ng of PCR product in 10 μ L volumes. CentriSep columns (Princeton Separations, Adelphia, NJ) were used to filter sequencing reactions, and samples were run on an ABI 377 automated sequencer using 5% Long Ranger polyacrylamide gels. Cyt b was sequenced in both the 5′ and 3′ directions, and edited and aligned with Sequencer 3.0 (GeneCodes, Ann Arbor, MI). Sequences were deposited in GenBank under accession numbers AY572046 to AY572108.

We used two methods of preliminary data exploration: descriptive diversity statistics and AMOVA. ML distances were computed between E. eschscholtzii and P. vandykei, between P. vandykei and P. idahoensis, within P. idahoensis, and within the Clearwater P. idahoensis. In addition, within P. idahoensis, we calculated several diversity statistics (nucleotide diversity [π], the number of pairwise differences [k], and the number of polymorphic sites [s]) using Arlequin 2.0 (Schneider et al., 2000). Analysis of molecular variance (AMOVA; Cockerham, 1969; Excoffier et al., 1992) was conducted using Arlequin to investigate population structure by testing the global hypothesis that haplotype variance is partitioned among drainages. We further grouped the seven drainages in which P. idahoensis occur into two groups, the Clearwater drainages (Selway, Lochsa, North Fork Clearwater) and the northern drainages (St. Joe, Coeur d'Alene, Kootenai, Columbia).

Phylogenetic Analyses

Phylogeny estimation

We generated maximum likelihood (ML) and maximum parsimony (MP) estimates of phylogeny to identify the major clades and to test phylogenetic hypotheses in the P. vandykei species group. We pruned all redundant haplotypes, and performed searches with Paup* 4.0 (Swofford, 2002). Under parsimony, heuristic searches included 150 random addition-sequence replicates and TBR branch swapping. Characters were equally weighted, and there was no limit to the number of trees retained. Nodal support was assessed using bootstrap analysis with 1000 replicates (Felsenstein, 1985). We used Dt-ModSel (Minin et al., 2003) to select a model; this method incorporates fit, a penalty for overparameterization, and performance into model selection, and is expected to select simpler models that nevertheless provide better branch-length estimates than other alternative model-selection methods. Because model-selection methods cannot assess the adequacy of chosen models, we evaluated the absolute goodness of fit of the chosen model chosen using the parametric bootstrap approach (e.g., Goldman, 1993; Sullivan et al., 2000). We then conducted heuristic searches under ML with the chosen model, and TBR branch swapping, 10 random-addition replicates, and a random starting tree. Nodal support for the ML tree was assessed using nonparametric bootstrap values (Felsenstein, 1985), computed from 100 replicates.

Frequentist hypothesis testing

We used parametric bootstrap analyses (Huelsenbeck et al., 1996) to test the ancient vicariance and recent dispersal hypotheses. Because ancient vicariance predicts reciprocal monophyly of P. idahoensis and P. vandykei, we searched for the best ML topology that refutes reciprocal monophyly of P. idahoensis and P. vandykei (i.e., a converse constraint). We then used Seq-Gen (Rambaut and Grassley, 1997) to simulate 100 data sets on the best constrained tree, with model parameters and branch lengths estimated using the ML model selected above. Each simulated data set was then searched with Paup* for both the best tree and the best tree constrained for nonreciprocal monophyly. The distribution across all replicates of the difference in log-likelihood scores (lnLunconstrainedlnLconstrained) served as the null distribution, which was used to assess the difference calculated from the real data. The recent dispersal hypothesis was tested in a similar fashion. This hypothesis proposes that inland populations are the result of dispersal from the Cascades through the Okanogan Highlands (Fig. 1), and predicts that haplotypes from the Columbia River Drainage will be basal to other P. idahoensis. Therefore the tree used to test the recent dispersal hypothesis was found by constraining the haplotypes from the Columbian drainage to be basal within the P. idahoensis clade.

Bayesian hypothesis testing

We used MrBayes (Huelsenbeck and Ronquist, 2001,) to assess the posterior probability of each of the two a priori hypotheses (i.e., ancient vicariance and inland dispersal). This approach attempts to determine the posterior probability that some hypothesis is correct, given the model and the data, by using MCMC to estimate the joint posterior probability distribution. It further assumes that the Markov chain is sampling from the true posterior distribution (i.e., the chain has reached stationarity with respect to the parameter of interest; Huelsenbeck et al., 2002) and, because we are assessing conformation to topological predictions in these tests, stationarity with respect to topology is critical. Because the topology parameter may be particularly susceptible to non-stationarity (Huelsenbeck et al., 2002), we employed a stationarity test similar to one proposed by Huelsenbeck et al. (2002). We conducted four independent heated runs (each composed of four Metropolis-coupled chains), and started each run with a different random tree. We ran the chains for 107 generations and sampled trees every 1000 generations. If the four independent runs have each converged on the true joint posterior probability distribution, the four samples of trees should represent independent samples drawn from that distribution. In order to assess this statistically, we saved the last 5000 trees from each run and computed the symmetric-difference distance between each tree in the sample and the ML tree (found above) using Paup*. We then conducted a standard analysis of variance (ANOVA) on the four groups of tree-to-tree distances to assess whether the four chains sampled independently from the same underlying (and hopefully true) joint posterior probability distribution. Although a nonsignificant result for this test would not guarantee that the runs have reached stationarity with respect to topology, it would provide much stronger evidence of such than would the standard examination of lnL plots. To complete the hypothesis test, we then imported the sample of trees from the Bayesian analysis into Paup* and filtered them with the same constraint trees used to test the ancient vicariance and recent dispersal hypotheses. The proportion of trees in the sample consistent with the topology predicted by each hypothesis is the Bayesian posterior probability of that hypothesis.

Population Genetic Analyses

MDIV analysis

The ancient vicariance hypothesis predicts that the coalescent depth of divergence between P. vandykei and P. idahoensis will be great due to their isolation since the xerification of the Columbia basin began some 5 Mya. To investigate the depth of divergence, we analyzed the data with MDIV (Neilson and Wakeley, 2001), an MCMC-based program designed to distinguish migration from isolation. We performed several iterations of the analysis before we were satisfied that we were analyzing the data with appropriate priors on M and T. Our initial analysis of P. vandykei and P. idahoensis, with M = 10 and T = 5, indicated that the migration rate was close to zero. This was congruent with the results of the phylogenetic hypothesis testing, so we reanalyzed the data with prior values of M = 0 and T = 6. We conducted 3 × 106 generations of the Markov chain, and then repeated the analysis using a different seed value to assess stationarity. We assumed a mutation rate of 10− 7 when converting the estimate of θ to Ne, which was then used to convert the posterior distribution of TMRCA to years.

The a priori corollary hypothesis of ancient vicariance, the Clearwater refugium hypothesis (Daubenmire, 1975), predicts a post-Pleistocene population expansion from a refugium located in the Clearwater Drainage. We used two methods, nested clade analysis (NCA) and coalescent estimates of the rate of population growth (g) to investigate P. idahoensis for the genetic signature of recent population expansion.

Nested clade analysis

Of all available population genetic approaches, NCA is the method that most explicitly incorporates geography. We first used NCA to assess the significance of the correlation of haplotypes and geography. Unweighted geographic distances were calculated from latitude and longitude, and TCS (Clement et al., 2001) was used to estimate a haplotype network for P. idahoensis. Groupings in this network were nested following Templeton (1998) and ambiguities were resolved following the guidelines in Templeton et al. (1987). GeoDis (Posada et al., 2000) was used to calculate the average geographical distance of the haplotypes to the center of each clade (Dc), the average distance of the haplotypes in each clade to the center of the clade at the next nesting level (Dn), and the difference between tip clades and interior clades for Dc and Dn. We tested for the significant association of haplotype with geography using a categorical permutation contingency analysis, and used the inference key (Templeton, 2004) to infer the likely population process that contributed to the pattern of genetic diversity observed within P. idahoensis.

Under the Clearwater refugium hypothesis, population expansion should result in little association of haplotype with geography because there has been insufficient time for lineage sorting to occur among populations derived from the same source, and the NCA should fail to reject the null hypothesis of no association between haplotype and geography. By contrast, upon inclusion only of samples from the putative refugium, the hypothesis predicts a significant association haplotype and geography because these samples represent populations descended from long-isolated populations. As a test of the Clearwater refugium hypothesis, we performed a second NCA using only the samples from the Clearwater drainage to test for a significant association of haplotype and geography.

Coalescent simulations

Population expansion can also be assessed using methods that directly estimate parameters of an expanding population by including the growth rate (g) in a coalescent model. Such genealogical methods are more efficient than methods based on diversity statistics (Kuhner et al., 1998), particularly when the g is low or negative. Expanding populations produce genealogies that are dominated by the tips of the tree, because most of the coalescent events occur in the recent past after the population expansion. The coalescent model for expanding populations of Kuhner et al. (1998) rescales the total time of the coalescent process so that g = (1 −egt)/T, where t is the original time to the common ancestor and T is the rescaled time. A significantly large g indicates population expansion. We used the program Fluctuate (Kuhner et al., 1998), which implements a Metropolis-Hastings geneology sampler to estimate g for the entire sample of P. idahoensis (gT) and the Clearwater samples alone (gC). The hypothesis predicts that gT will be significantly greater than zero, because P. idahoensis has expanded into recently deglaciated areas. Genealogical samplers rely on the shape of the likelihood surface to generate confidence intervals (CIs), and under some conditions these CIs have been shown to be too narrow for parameters of structured coalescent models (Abdo et al., 2004). In order to generate a null distribution of gT, we conducted coalescent simulations under a constant population size model (e.g., gi = 0), with the effective population size (Ne) and mutation rate (μ) parameters estimated from the original data. Because these parameters are interdependent, we used the compound parameter θ, which is equal to 2Neμ for mitochondrial data, in the simulations. We used Migrate (Beerli, 2002) to estimate θ for the Clearwater samples (θC) and for all the samples (θT). We then simulated 100 data sets with θC and 100 data sets with θT using Treevolve (Grassly et al., 1999) and a model of constant population size. All Markov chains in Fluctuate were initialized on the empirically estimated values of θ, but were allowed to move into regions of higher log-likelihood, and we used the estimates of gi from each of these simulated data sets to construct the distribution of gT and gC under the null hypothesis of stable populations.

Analyses within drainages

We explored genetic structure within different populations of P. idahoensis by calculating several diversity statistics for samples collected within a single river drainage. These drainages likely contain isolated populations of P. idahoensis because most are separated by ridge lines that exceed the recorded elevation limits of P. idahoensis (Nussbaum et al., 1983). Within each drainage, genetic diversity was estimated by calculating nucleotide diversity (π), the number of pairwise differences (k), and the number of polymorphic sites (s) using Arlequin 2.0 (Schneider et al., 2000). These parameters were used in the calculation of Tajima's D (Tajima, 1989), which can be used to test for population expansion (Rogers, 1995). We also tested for population expansion using the demographic expansion model of Rogers and Harpending (1992) to analyze the pairwise mismatch distributions. Populations that have been stable historically are predicted to have multimodal mismatch distributions, whereas those that have undergone a recent expansion are predicted to be unimodal. We calculated the mismatch distribution with Arlequin (Schneider et al., 2000) for samples collected within each of the eight river drainages in which P. idahoensis occurs, and used the Raggedness Index of Harpending (1994) to identify populations in which the demographic model of a unimodal mismatch distribution can not be rejected.

Results

We sequenced 669 nucleotides, corresponding to positions 16294 to 16963 of the Xenopus laevis mitochondrial gene Cyt b (Roe et al., 1985). The pattern of molecular evolution was consistent with mitochondrial DNA of plethodontids and inconsistent with the presence of nuclear pseudogenes. Seven haplotypes were sampled in P. vandykei, and 62 haplotypes were sampled in P. idahoensis (Table 2). One of the haplotypes (A) was found in 45% of all P. idahoensis individuals sampled (109/243), and was sampled in 31 of the 49 populations.

Figure 3

ML phylogeny estimate for the P. idahoensis/P. vandykei complex under a HKY+Γ model of sequence evolution. This is one of two ML trees with a log-likelihood score of −2511.685. ML bootstrap values from 100 replicates are shown for nodes recovered in > 50% of the bootstrap data sets. Haplotypes are grouped by drainage, and shown in shaded boxes, but haplotype ‘A’ is found in the North Fork of the Clearwater, the St. Joe, the Coeur d'Alene, the Kootenai, and the Columbia drainages as well as the Lake Coeur d'Alene population. Haplotype ‘N’ is found on both the Kootenai and the Columbia Rivers, but is sister to a haplotype (haplotype R) found only on the Kootenai River. Haplotype ‘O’ is found on the St. Joe, but is sister to haplotypes found on the Columbia River. Haplotype ‘bcc255’ is found on the Kootenai River, but is sister to a clade found on the North Fork of the Clearwater River.

Figure 3

ML phylogeny estimate for the P. idahoensis/P. vandykei complex under a HKY+Γ model of sequence evolution. This is one of two ML trees with a log-likelihood score of −2511.685. ML bootstrap values from 100 replicates are shown for nodes recovered in > 50% of the bootstrap data sets. Haplotypes are grouped by drainage, and shown in shaded boxes, but haplotype ‘A’ is found in the North Fork of the Clearwater, the St. Joe, the Coeur d'Alene, the Kootenai, and the Columbia drainages as well as the Lake Coeur d'Alene population. Haplotype ‘N’ is found on both the Kootenai and the Columbia Rivers, but is sister to a haplotype (haplotype R) found only on the Kootenai River. Haplotype ‘O’ is found on the St. Joe, but is sister to haplotypes found on the Columbia River. Haplotype ‘bcc255’ is found on the Kootenai River, but is sister to a clade found on the North Fork of the Clearwater River.

Table 2.

Number of haplotypes by drainage for Plethodon idahoensis. Haplotypes sampled in multiple individuals are labeled A to V (M and U were found in P. vandykei). The final column, denoted with the *, lists the number of individuals from each drainage with a haplotype not sampled in any other P. idahoensis individual. Haplotype letters correspond to those in Figures 2 and 3

 Haplotype 
 
 
Drainage 
Selwaya · · 12 · · · · · · · · · · · · · · · · · 
Lochsaa · · · 10 · · · · · · · · · · · · · · · 
NFCa,b 11 · · · · · · · · · · · · · 11 
St Joe 40 · · · · · · · · · · · · · · · · 
Cd'Ac 18 · · · · · · · · · · · · · · · · · · · 
LkCd'A · · · · · · · · · · · · · · · · · · · 
Kootenai 16 · · · · · · · · · · · · · 10 
Columbia 15 · · · · · · · · · · · · · · · · 
 Haplotype 
 
 
Drainage 
Selwaya · · 12 · · · · · · · · · · · · · · · · · 
Lochsaa · · · 10 · · · · · · · · · · · · · · · 
NFCa,b 11 · · · · · · · · · · · · · 11 
St Joe 40 · · · · · · · · · · · · · · · · 
Cd'Ac 18 · · · · · · · · · · · · · · · · · · · 
LkCd'A · · · · · · · · · · · · · · · · · · · 
Kootenai 16 · · · · · · · · · · · · · 10 
Columbia 15 · · · · · · · · · · · · · · · · 
a

NFC, Lochsa, and Selway are the major tributaries of the Clearwater drainage.

b

NFC = North Fork of the Clearwater.

c

Coeur d'Alene River.

Minimum and maximum ML distances (HKY+Γ) were as follows: E. eschscholtziiP. vandykei (min = 0.26411 substitutions/site; max = 0.27034), P. vandykeiP. idahoensis (min = 0.08529; max = 0.10604), within P. idahoensis (max = 0.02776), and within Clearwater P. idahoensis (max = 0.02329). With the exception of the Lake Coeur d'Alene sample (which was fixed for haplotype A), all populations of P. idahoensis were polymorphic for Cyt b, with the largest number of nucleotide polymorphism (s) occurring on the Kootenai River Drainage (Table 1). The AMOVA suggested that most of the molecular variance occurred within drainages (51.69%; P < 0.001) and between the southern and northern drainages (37.17%; P < 0.001), with only little variation among drainages (11.14%; P = 0.084).

Phylogenetic Analyses

Phylogeny estimation

We selected the HKY+Γ model of sequence evolution using Dt-ModSel (Minin et al., 2003) with equilibrium base frequencies of πA = 0.305; πC = 0.235; πG = 0.14; πT = 0.32; transition-transversion ratio = 1.336); and Γ -distribution shape parameter (α = 0.716. Furthermore, the absolute goodness-of-fit test indicated that the HKY+Γ model could not be rejected (P = 0.64), indicating that the parametric bootstraps and Bayesian hypothesis tests should not be compromised by reliance on a single substitution model.

There were 91 parsimony-informative sites in the complete data set. We found 2160 MP trees (not shown), with a tree length of 278. The topology of the MP trees was similar to that of two ML trees with log-likelihood scores − lnL = 2511.685. The ML topologies differed from one another in only minor ways, and one of these ML trees is shown (Fig. 3). In the ML estimate of phylogeny, P. idahoensis and P. vandykei form reciprocally monophyletic clades that are well supported (100% bsML; 99% bsMP; 1.0 Bpp) in both ML trees and MP trees. Within P. idahoensis, haplotypes from the Selway River (a tributary of the Clearwater) were monophyletic in the optimal tree, but this clade was weakly supported (61% bsML; 54% bsMP; 0.429 Bpp). Most of the haplotypes from the Lochsa River (another Clearwater tributary) formed a sister clade to the Selway haplotypes in the optimal topology, but this relationship was also poorly supported (60% bsML; 46% bsMP; 0.261 Bpp). One haplotype (bcc66), collected two miles downstream from the confluence of the Selway and Lochsa Rivers, was basal to a large and poorly resolved clade containing all the haplotypes sampled in drainages to the north.

Frequentist hypothesis testing

The best topology on which the two species are not reciprocally monophyletic (thereby refuting the ancient vicariance hypothesis) was 9.33 −lnL units worse than the ML tree, and when we constrained Paup* to search for the best topology consistent with the inland dispersal hypothesis, there was a deterioration of 14.323 −lnL units. The parametric bootstrap tests indicated that both values are higher than can be attributed to phylogenetic uncertainty (PNRM < 0.01; PRD < 0.01), supporting the establishment of the disjunction by an ancient vicariance event

Bayesian hypothesis testing

The four independent Bayesian runs had average tree-to-tree distances from the ML tree of 54.23, 53.90, 53.87, and 53.87. The ANOVA indicated that these values were not significantly different (FOBS = 2.28; 0.1 > P > 0.05), a result which we interpreted as evidence that the four MrBayes runs were sampling topologies from the same joint posterior probability distribution. Bayesian posterior probabilities supported the ancient vicariance hypothesis (P > 0.99995), and rejected the recent dispersal hypothesis (P < 0.00005).

Population Genetic Analyses

Estimation of population divergence

Use of MDIV (Nielson and Wakely, 2001) to estimate the joint distribution of theta (θ) and the time since population divergence (Tdiv) resulted in estimates of θ = 0.03589 and Tdiv = 2.088 * 2Ne generations. Assuming μ = 10−7, the effective population size Ne = 179450. Converting the estimate of Tdiv into years requires a value for the generation length, and no good estimate is available for P. vandykei/P. idahoensis. Nussbaum et al. (1983) report that P. idahoensis do not begin breeding behavior until their fifth year, so if we consider 5 years to be a lower bound for the generation length, our MDIV-derived estimate of TMRCA = 2.088 * 2Ne generations corresponds to a divergence time of 3.75× 106 years with a 95% credibility interval of 1.31 × 106 to 8.2 × 106 years. The maximum predicted time of lineage divergence under the ancient vicariance hypothesis (5.0 × 106 years) falls within the 95% credibility interval of our Tdiv estimate, and the confidence interval on Tdiv excludes the post Pleistocene divergence expected under the recent dispersal hypothesis.

Nested clade analysis

In the NCA of the entire P idahoensis sample, we were not able to reject a random association of haplotypes with geography (P = 0.284), although several subclades (3–X, 2–10, 1–31, 1–25, 1–15, 1–12, and 1–10; Fig. 4) did exhibit geographically structured variation. However, we were able to detect significant geographic structure using NCA when we examined just the samples from the Clearwater drainage (P < 0.001), a result consistent with the Clearwater refugium hypothesis (a corollary to the ancient vicariance hypothesis accepted above). In each of the nested clades where there was a significant association of haplotype and geography, the inference that we reached by following the inference key (Templeton, 2004) was one of population expansion.

Figure 4

Haplotype network among all P. idahoensis haplotypes. Haplotypes are connected by lines that represent a single substitution, and missing haplotypes are labeled with a black dot. Nesting levels for the NCA are shown, with the one-step clades in white, the two-step clades in light gray, and the three-step clades in dark gray. The clade numbers are also shown, with clades that have a significant association of haplotype and geography denoted with an asterisk.

Figure 4

Haplotype network among all P. idahoensis haplotypes. Haplotypes are connected by lines that represent a single substitution, and missing haplotypes are labeled with a black dot. Nesting levels for the NCA are shown, with the one-step clades in white, the two-step clades in light gray, and the three-step clades in dark gray. The clade numbers are also shown, with clades that have a significant association of haplotype and geography denoted with an asterisk.

Coalescent simulations

Theta estimates, calculated with the assumption of a stable population, were θT = 0.0473 for all P. idahoensis and θC = 0.0145 for the Clearwater samples. When the coalescent model allowed for population expansion, both the total sample and the Clearwater sample exhibited high growth rate estimates (gC = 435.02; gT = 429.02). Simulations indicated that both gC and gT were significantly greater than zero (PC < 0.01; PT < 0.01). These results are consistent with the hypothesis of a population expansion both within the Clearwater (the proposed refugium) as well across the entire species.

Analyses within drainages

When we pooled samples by drainage and used Tajima's D and mismatch distributions to infer population expansion within the drainages, we found that there is evidence in each of the drainages for a recent population expansion (Table 1), although in two drainages the value of Tajima's D was not significant.

Discussion

Methodological Issues

Phylogeographic data exist at the nexus of phylogenetics and population genetics, and analyses of such data are best conducted using techniques developed for each. This study demonstrates the utility of combining statistical phylogeographic and population genetic approaches in assessing nested a priori hypotheses derived using data from geology, paleoecology, and paleobotany. Such hypotheses predict specific patterns of genetic variation that can be tested while simultaneously accounting for the stochastic nature of the coalescent process and the amplifications of these effects over geological time. Phylogenetic trees are extremely useful for testing the predictions of hypotheses that relate to the deepest genetic structure that may be present in a phylogeographic data set. For example, our ML estimate of phylogeny supports both the ancient vicariance and Clearwater refugium hypotheses and the statistical tests such as the parametric bootstraps or Bayesian posterior tests eliminate phylogenetic uncertainty as an explanation. The congruence between the parametric bootstrap tests and the Bayesian hypothesis tests, coupled with the greater speed of the latter, suggest that Bayesian hypothesis testing may be a more efficient use of computational resources, although power tests are needed to explore the issue. Similarly, both have been criticized as relying too heavily on the adequacy of a particular model of sequence evolution (e.g., Felsenstein et al., 2003; Buckley, 2001). This emphasizes the need to both choose a model rationally and to assess the absolute fit of the data to the model. In our case, the chosen model adequately approximates the unknown process that generated the data.

Due primarily to the stochasticity of the coalescent process, it is considerably more difficult to design phylogenetic tests of recent events; it is not possible to construct appropriate predicted topologies in order to test hypotheses regarding events in the recent past such as post-Pleistocene expansion. Such hypotheses are best evaluated with predictions derived from population genetic theory and test statistics that account for coalescent stochasticity (Knowles and Maddison, 2002). We have investigated post-Pleistocene expansion using several complementary methods. On the among-populations scale, NCA was used to test the prediction that haplotypes within the Clearwater drainage (the location of the putative refugium) would have a significant association with geography. Because NCA cannot be used to test an a priori null hypothesis of a rapid expansion (but can identify a signature of population expansion at specific nested clades), the significance of our empirical estimate of the growth rate parameter (g) was assessed by simulating a null distribution. On smaller geographic scales, we validated the expectation of population expansion using Tajima's D and analysis of mismatch distributions.

Perhaps the greatest difficulty of the statistical phylogeographic approach is choosing which test(s) are most appropriate to the question being asked; each of the tests used here has shortcomings that may lead to misinterpretation. One advantage to utilizing a variety of tests designed to detect signal across a wide range of genetic divergence is that the tests can be cross validated. For example, the standard interpretation of a significantly negative Tajima's D is that selection and/or population expansion are present (Tajima, 1989). Similarly, both sampling and unaccounted for genetic structure can lead to spurious inferences of Tajima's D (Hammer et al., 2003). Nevertheless, we interpret the signal present in our data as population expansion because this result was also found when we estimated the population growth rate and calculated the mismatch distributions. In two drainages, we had non-negative values of Tajima's D but could not reject the unimodal mismatch distribution. We suspect that this discrepancy is related to the low values of π for these drainages; π is incorporated into Tajima's D but not the distribution of mismatches. The NCA results were consistent with our prediction that extant populations within the refugium would retain significant geographic structure; there was a significant association of geography and haplotypes in the Clearwater, and this association degraded when all samples were analyzed. This apparently results from insufficient time for lineage sorting to have occurred in the northern portion of the distribution. In each nested clade that had a significant association of haplotype and geography, the inference derived by use of Templeton's (2004) inference key was one of population expansion (results not shown). This may be one population process that NCA is particularly capable of inferring (Templeton, 1998). The appeal of NCA is that it represents perhaps the only analytical framework that attempts to infer genetic processes ranging from small scale population processes to those acting across a complex geography and on a longer time scale. However, we agree with Knowles and Maddison (2002) that further study is needed on the circumstances under which NCA works well and where it fails. Pending the results of that research, there is a need to cross-validate inferences made with NCA using other data (Templeton, 2004) and other methods.

Genetics of the Plethodon vandykei Species Group

We found deep divergence between P. vandykei and P. idahoensis with mitochondrial data. The depth of this divergence (corrected distances of 0.0853 to 0.106 substitutions per site), as well as the results of the parametric bootstrap Bayesian hypothesis tests, and coalescent-based estimates of TMRCA all support the ancient vicariance hypothesis. The corollary of the ancient vicariance hypothesis is that P. idahoensis persisted throughout the Pleistocene glaciations in a mesic forest refugium, putatively situated in the Clearwater drainage, and subsequently expanded into deglaciated drainages at the close of the Pleistocene. The results of the NCA, as well as our coalescent-based estimates of the growth rate, support recent population expansion. Furthermore, our data suggest that extant populations within the refugial Clearwater drainage have undergone significant post-Pleistocene growth. One possible explanation is that refugial mesic forests within the Clearwater may have been elevationally compressed during Pleistocene glacial maxima, and the retreat of glaciation at the close of the Pleistocene led to an increase of available habitat within the drainage. The three major tributaries that comprise the Clearwater drainage (the North Fork, Lochsa, and Selway Rivers) may also have functioned as independent refugia during Pleistocene cold periods. These results suggest that our understanding of the Clearwater refugium is preliminary, and that more refined hypotheses are needed to understand the dynamics of this complex refugium. We are currently screening multilocus SNP data to further test the hypotheses that we have refined here. Several authors have suggested that these data offer significant advantages for phylogeographic studies (Arbogast et al., 2002; Brumfield et al., 2003; Funk and Omland, 2003).

Ecosystem Phylogeography

These results increase our understanding of the assembly of the inland mesic forest ecosystem in the Pacific Northwest. Three other vertebrate lineages exhibit disjunct mesic forest distributions: tailed frogs (Ascaphustruei/A. montanus), giant salamanders (Dicamptodon spp.), and water voles (Microtus richardsoni). Studies of these taxa suggest that the vertebrate components of this ecosystem were not assembled in a concerted manner. Good (1989) studied genetic diversity in Dicamptodon using allozymes and found that the deepest divergence was between coastal and inland populations. He considered this divergence deep enough to recommend recognition of the inland forms as a distinct species (Good, 1989). Nielsen et al. (2001) found a similar pattern in Ascaphus using data from two mitochondrial genes, and also argued for the recognition of the inland Ascaphus montanus as a distinct species. Although the distribution of these lineages, as well as that of the P. vandykei species group, is best explained by an ancient vicariance and persistence through the Pleistocene glaciations, the pattern found in Microtus richardsoni suggests a northern dispersal (J. R. Demboski and J. Sullivan, unpublished). Demboski et al. (in preparation) have found shallow genetic divergence in the water vole, and haplotypes that are shared between populations from the northern Cascades and the northern Rocky Mountains. This northern dispersal also appears to explain the distribution of white-bark pine (Pinus albicaulis; Richardson et al., 2002), a subalpine disjunct.

Based on certain common life history traits, the similarity in the genetic pattern between tailed frogs and Pacific giant salamanders might be expected. Adults of each have lungs and live in the wet habitat surrounding high gradient mountain streams and both groups have larval stages that utilize the same streams. However, the life history of Plethodon vandykei and P. idahoensis is quite different; the genus Plethodon is characterized by direct development (no free-living larval stage), and is lungless and dependent on ambient moisture for cutaneous respiration. Although P. vandykei and P. idahoensis occupy the same broad ecosystem as Ascaphus and Dicamptodon, they primarily exploit rocky seeps and talus slopes rather than high gradient streams. They also have much lower reproductive capacities, laying 4 to 12 eggs per year, relative to Dicamptodon, which lays 135 to 200 eggs per year, or Ascaphus, which lays 44 to 75 eggs per year (Nussbaum et al., 1983). In light of these life histories, the differences among the northern distributional limits of the three inland species is striking. Dicamptodon aterrimus is unknown from north of the St. Joe River Drainage (see Fig. 1), whereas Ascaphus montanus is known from the extreme southeastern portion of British Columbia, apparently reaching no more than 50 km into Canada (Dupruis, 2002). By contrast, Plethodon idahoensis has recently been found to extend 95 km north of Revelstoke, British Columbia (our locality no. 49; Wilson and Ohanjanian, 2003), which is several hundred kilometers north of the international border and 300 km farther north than Ascaphus montanus reaches. The northern extent of our sampling indicates that the range of P. idahoensis has expanded over 600 km since the close of the Pleistocene. This rate of expansion, > 50 m annually, is striking in comparison to the slower apparent rates of A. montanus and D. aterrimus (assuming both were present in the Clearwater Refugium; genetic data for each suggest this; Nielson et al., in preparation; Carstens et al., in preparation). The best explanation for this rapid expansion may be the ability of P. idahoensis to exploit ephemeral habitat created in the wake of retreating glaciers. The wet moraines deposited by retreating glaciers are in many ways ideal habitat for P. idahoensis, and deglaciation may have created temporary colonization routes that allowed them to expand their range to a greater degree than the other codistributed amphibians.

Conclusions

The advantages of statistical phylogeography may be most pronounced in the comparison across all the taxa that comprise a given ecosystem. Organisms as diverse as bryophytes and conifers, millipedes and salamanders, chipmunks and fungi have vastly different life histories and will necessarily be studied with data from divergent genomes. Such data can best be synthesized into ecosystem phylogeography through rigorous statistical tests of independent a priori hypotheses. It is this commonality of such testable hypotheses that will lead to a significant improvement in our understanding of regional phylogeography and the historical assembly of ecosystems.

Acknowledgments

We received critical advice from A. Wilson concerning locations to search for P. idahoensis and P. vandykei. We would also like to thank P. Ohanjanian, M.-A.Beauche, C. Burch, C. Steele, K. Villa, and T. Antifeau for assistance in collecting samples. K. Segraves, C. Simon, J. Sites, J. Good, and two anonymous reviewers provided valuable comments on earlier versions of this manuscript. Additionally, S. Brunsfeld, L. Forney, O. Pellmyr, D. Althoff, C. Smith, K. Berger, J. Good, B. Hyde, and D. Rokyta participated in important discussions pertaining to various portions of this research. Funding for B. Carstens was provided by the University of Idaho Presidential Fellowship. Funding for computer hardware used for much of the computational aspect of this work was provided with funding from NIH National Center for Research Resources grants P20 RR16454 for the BRIN Program, and P20 RR16448 for the COBRE program. Both these programs are part of the University of Idaho's Initiative in Bioinformatics and Evolutionary studies (IBEST). Funding for J. Degenhardt was provided by an NSF REU. Funding used for laboratory reagents and supplies was provided by NSF EPSCoR 809935 (to JS). JS was supported by NSF DEB-9974124 while this research was being conducted.

References

Abdo
Z.
Crandall
K. A.
Joyce
P.
Evaluating the performance of likelihood methods for detecting population structure and migration
Mol. Ecol.
 , 
2004
, vol. 
4
 (pg. 
837
-
851
)
Arbogast
B. S.
Edwards
S. V.
Wakeley
J.
Beerli
P.
Slowinski
J. B.
Estimating divergence times from molecular data on phylogenetic and population genetic timescales
Ann. Rev. Ecol. Syst.
 , 
2002
, vol. 
33
 (pg. 
707
-
740
)
Beerli
P.
Migrate: Documentation and program, part of Lamarck. Version. 1.5
2002
 
Brumfield
R. T.
Beerli
P.
Nickerson
D. A.
Edwards
S. V.
The utility of single nucleotide polymorphisms in inferences of population history
Tr. Ecol. Evol.
 , 
2003
, vol. 
18
 (pg. 
249
-
256
)
Brunsfeld
S. J.
Sullivan
J.
Soltis
D. E.
Soltis
P. S.
Silvertown
J.
Antonovics
J.
Comparative phylogeography of northwestern North America: A synthesis
Integrating ecology and evolution in a spatial context
 , 
2001
Williston, Vermont
Blackwell Publishing
(pg. 
319
-
339
Pages
Buckley
T. R.
Model misspecification and probabilistic tests of topology: Evidence from empirical data sets
Syst. Bio.
 , 
2002
, vol. 
51
 (pg. 
509
-
523
)
Clement
M.
Posada
D.
Crandall
K.
TCS: A computer program to estimate gene genealogies
Mol. Ecol.
 , 
2001
, vol. 
9
 (pg. 
1657
-
1660
)
Cockerham
C. C.
Variance of gene frequencies
Evolution
 , 
1969
, vol. 
23
 (pg. 
72
-
83
)
Crandall
K. A.
Templeton
A. R.
Harvey
P. H.
Brown
A. J. L.
Maynard Smith
J.
Nee
S.
Applications of intraspecific phylogenetics
New uses for new phylogenies
 , 
1996
New York, New York
Oxford University Press
(pg. 
81
-
102
Pages
Daubenmire
R.
Davis
R. J.
Plant geography of Idaho
Flora of Idaho
 , 
1952
Provo, Utah
Brigham Young University Press
(pg. 
1
-
17
Pages
Daubenmire
R.
Floristic plant geography of eastern Washington and northern Idaho
1975
Provo, Utah
Brigham Young University Press
Delcourt
P. A.
Delcourt
H. R.
Paleoclimates, paleovegetation, and paleofloras during the late Quaternary
Flora of North America Editorial Committee
 , 
1993
, vol. 
volume 1
 
first edn.
New York, New York
Oxford University Press
(pg. 
71
-
94
Pages
Dupruis
L.
Distribution of Ascaphus montanus in the Yahk River and neighboring watersheds
2002
Box 612, Squamish
British Columbia
 
Report to Tenbec Industries and the Columbia Basin Fish and Wildlife Compensation Program. Ascaphus Consulting, VON 3GO
Excoffier
L.
Smouse
P.
Quattro
J.
Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data
Genetics
 , 
1992
, vol. 
136
 (pg. 
343
-
359
)
Felsenstein
J.
Confidence limits on phylogenies: An approach using the bootstrap
Evolution
 , 
1985
, vol. 
39
 (pg. 
783
-
791
)
Felsenstein
J.
Inferring phylogenies
2003
Sunderland, Massachusetts
Sinauer
Funk
D. J.
Omland
K. E.
Species-level paraphylly and polyphyly: Frequency, causes, and consequences, with insights from animal mitochondrial DNA
Ann. Rev. Ecol. Syst.
 , 
2003.
, vol. 
34
 (pg. 
397
-
423
)
Goldman
N. J.
Statistical tests of models of DNA substitution
J. Mol. Evol
 , 
1993
, vol. 
36
 (pg. 
182
-
198
)
Good
D. A.
Hybridization and cryptic species in Dicamptodon (Caudata: Dicamptodontidae)
Evolution
 , 
1989
, vol. 
43
 (pg. 
728
-
744
)
Graham
A
History of the vegetation, Cretaceous—Tertiary
Flora of North America Editorial Committee
 , 
1993
, vol. 
volume 1
 
first edn.
New York, New York
Oxford University Press
(pg. 
57
-
70
Pages
Grassley
N. C.
Harvey
P. H.
Holmes
E. C.
Population dynamics of HIV-1 inferred from gene sequences
Genetics
 , 
1999
, vol. 
151
 (pg. 
427
-
438
)
Gysels
E. S.
Hellemans
B.
Papoulie
C.
Volckaert
F. A. M.
Phylogeography of the common goby, Pomatoschistus microps, with particular emphasis on the colonization of the Mediterranean and the North Sea
Mol. Ecol.
 , 
2004
, vol. 
13
 (pg. 
403
-
417
)
Harpending
R. C.
Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution
Hum. Biol.
 , 
1994
, vol. 
66
 (pg. 
591
-
600
)
Howard
J. H. L.
Seeb
W.
Wallace
R.
Genetic variation and population divergence in the Plethodon vandykei species group (Caudata: Plethodontidae)
Herptologica
 , 
1993
, vol. 
49
 (pg. 
238
-
247
)
Huelsenbeck
J. P.
Hillis
D. M.
Jones
R.
Ferraris
J. D.
Palumbi
S. R.
Parametric bootstrapping in molecular phylogenetics: Applications and performance
Molecular zoology: Advances, strategies, and protocols
 , 
1996
New York, New York
Wiley-Liss
(pg. 
19
-
45
Pages
Huelsenbeck
J. P.
Lagret
B.
Miller
R. E.
Ronquist
F.
Potential applications and pitfalls of Bayesian inference of phylogeny
Syst. Biol.
 , 
2002
, vol. 
51
 (pg. 
673
-
688
)
Huelsenbeck
J. P.
Ronquist
F.
MrBayes: Bayesian inference of phylogenetic trees
Bioinformatics
 , 
2001
, vol. 
17
 (pg. 
754
-
755
)
Knowles
L. L.
Maddison
W. P.
Statistical phylogeography
Mol. Ecol.
 , 
2002
, vol. 
11
 (pg. 
2623
-
2635
)
Kuhner
M. K.
Yamato
J.
Felsenstein
J.
Maximum likelihood estimation of population growth rates based on the coalescent
Genetics
 , 
1998
, vol. 
149
 (pg. 
429
-
434
)
Lyell
C.
Principles of geology, volume 2
1832
London
John Murray
 
(Facsimile edition published by University of Chicago Press, 1991.)
Mack
R. N.
Rutter
N. W.
Bryant
V. M.
Jr.
Valastro
S.
Reexamination of post-glacial vegetation history in northern Idaho: Hager Pond, Bonner Co
Quat. Res.
 , 
1978
, vol. 
10
 (pg. 
241
-
255
)
Mahoney
M. J.
Molecular systematics of Plethodon and Aneides (Caudata: Plethodontidae: Plethodontini): Phylogenetic analysis of an old and rapid radiation
Mol. Phyl. Evol.
 , 
2001
, vol. 
18
 (pg. 
174
-
188
)
Mead
L. S.
Tilley
S. G.
Katz
L. A.
Genetic structure of the Blue Ridge dusky salamander (Desmognathus orestres): inferences from allozymes, mitochondrial DNA, and behavior
Evolution
 , 
2001
, vol. 
55
 (pg. 
2287
-
2302
)
Minin
V.
Abdo
Z.
Joyce
P.
Sullivan
J.
Performance-based selection of likelihood models for phylogeny estimation
Syst. Biol.
 , 
2003
, vol. 
52
 (pg. 
674
-
683
)
Moritz
C.
Schneider
C. J.
Wake
D. B.
Evolutionary relationships within the Ensatina eschsholtzii complex confirm the ring species interpretation
Syst. Biol.
 , 
1992
, vol. 
41
 (pg. 
273
-
291
)
Neilson
M.
Lohman
K.
Sullivan
J.
Phylogeography of the tailed frog (Ascaphus truei): Implications for the biogeography of the Pacific Northwest
Evolution
 , 
2001
, vol. 
55
 (pg. 
147
-
160
)
Nielsen
R.
Wakeley
J. W.
Distinguishing migration from isolation: An MCMC approach
Genetics
 , 
2001
, vol. 
158
 (pg. 
885
-
896
)
Nussbaum
R. A.
Brodie
E. D.
Jr.
Storm
R. M.
Amphibians and reptiles of the Pacific Northwest
1983
Moscow, Idaho
University of Idaho Press
Posada
D.
Crandall
K. A.
Templeton
A. R.
GeoDis: A program for the nested cladistic analysis of the geographical distribution of genetic haplotypes
Mol. Ecol.
 , 
2000
, vol. 
9
 (pg. 
487
-
488
)
Rambaut
A. E.
Grassly
N. C.
Seq-Gen: An application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees
Comp. Appl. Biosci.
 , 
1997
, vol. 
13
 (pg. 
235
-
238
)
Richardson
B. A.
Brunsfeld
S. J.
Klopfenstein
N. B.
DNA from bird-dispersed seed and wind-disseminated pollen provides insights into postglacial colonization and population genetic structure of whitebark pine (Pinus albicaulis)
Mol. Ecol.
 , 
2002
, vol. 
11
 (pg. 
215
-
227
)
Richmond
G. M.
Fryxell
R.
Neff
G. E.
Weis
P. L.
Wright
H. E.
Jr.
Frey
D. G.
The Cordilleran ice sheet of the northern Rocky Mountains and related Quaternary history
The Quarternary of the United States
 , 
1965
Princton, New Jersey
Princton University Press
(pg. 
231
-
242
Pages
Roe
B. A.
Ma
D. P.
Wilson
R. K.
Wong
J. F.
The complete nucleotide sequence of the Xenopus mitochondrial genome
J. Biol. Chem.
 , 
1985
, vol. 
260
 (pg. 
9759
-
9774
)
Rogers
A.
Genetic evidence for a Pleistocene population explosion
Evolution
 , 
1995
, vol. 
49
 (pg. 
608
-
615
)
Rogers
A. R.
Harpending
H.
Population growth makes waves in the distribution of pairwise genetic differences
Mol. Biol. Evol.
 , 
1992
, vol. 
9
 (pg. 
552
-
569
)
Saunders
I. W.
Tavare
S.
Watterson
G. A.
On the geneology of nested sub-samples from a haploid population
Adv. Appl. Prob.
 , 
1984
, vol. 
16
 (pg. 
471
-
491
)
Schneider
S.
Roessli
D.
Excoffier
L.
Arlequin v2.0.: Documentation and program
2000
 
Sullivan
J.
Arellano
E.
Rogers
D. S.
Comparative phylogeography of Mesoamerican highland rodents: Concerted versus independent responses to past climatic fluctuations
Am. Nat.
 , 
2000
, vol. 
155
 (pg. 
755
-
768
)
Sullivan
J.
Markert
J. A.
Kilpatrick
C. W.
Phylogeography and molecular sytematics of the Peromyscus aztecus species group (Rodentia: Muridae) inferred using parsimony and likelihood
Syst. Biol.
 , 
1997
, vol. 
46
 (pg. 
426
-
440
)
Swofford
D. L.
Paup*. Phylogenetic Analysis Using Parsimony (and Other Methods). Version 4
2002
Sunderland, Massachusetts
Sinauer Associates
Taberlet
P. L.
Fumagalli
A.
Wust-Saucy
G.
Cossons
J.-F.
Comparative phylogeography and postglacial colonization routes in Europe
Mol. Ecol
 , 
1998
, vol. 
7
 (pg. 
453
-
464
)
Tajima
F.
The effect of change in population size on DNA polymorphism
Genetics
 , 
1989
, vol. 
105
 (pg. 
437
-
460
)
Templeton
A. R.
Nested clade analysis of phylogeographic data: Testing hypotheses about gene flow and population history
Mol. Ecol
 , 
1998
, vol. 
7
 (pg. 
381
-
397
)
Templeton
A. R.
Statistical phylogeography: Methods of evaluating and minimizing inference errors
Mol. Ecol.
 , 
2004
, vol. 
4
 (pg. 
789
-
809
)
Templeton
A. R.
Boerwinkle
E.
Sing
C. F.
A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. I. Basic theory and analysis of alcohal dehydrogenase activity in Drosophila
Genetics
 , 
1987
, vol. 
117
 (pg. 
343
-
351
)
Waitt
R. B.
Jr.
Thorson
R. M.
Wright
H. E.
Jr.
Porter
S. C.
The Cordilleran ice sheet in Washington, Idaho, and Montana
Late Quaternary Environments of the United States, the Late Pleistocene
 , 
1983
Minneapolis, Minnesota
University of Minnesota Press
(pg. 
53
-
70
Pages
Wares
J. P.
Community genetics in the northwestern Atlantic intertidal
Mol. Ecol.
 , 
2002
, vol. 
11
 (pg. 
1131
-
1144
)
Wares
J. P.
Cunningham
C. W.
Phylogeography and historical ecology of the North Atlantic intertidal
Evolution.
 , 
2001
, vol. 
55
 (pg. 
2455
-
2469
)
Waterson
G. A.
Guess
H. A.
Is the most frequent allele the oldest? Theor
Pop. Biol.
 , 
1977
, vol. 
11
 (pg. 
141
-
160
)
Wilson
A. G.
Jr.
Larsen
J. H.
Jr.
Biogeographic analysis of the Coeur d'Alene Salamander (Plethodon idahoensis)
N.W. Sci.
 , 
1998a
, vol. 
72
 (pg. 
111
-
115
)
Wilson
A. G.
Jr.
Larsen
J. H.
Jr.
Morphometric analysis of salamanders of the Plethodon vandykei species group
Am. Mid. Nat.
 , 
1998b
, vol. 
141
 (pg. 
266
-
276
)
Wilson
A. G.
Jr.
Larsen
J. H.
Jr.
McAllister
K. R.
Distribution of Van Dyke's salamander (Plethodon vandykei Van Denburgh)
Am. Mid. Nat
 , 
1995
, vol. 
134
 (pg. 
388
-
393
)
Wilson
A. G.
Jr.
Ohanjanian
P.
Plethodon idahoensis Slater and Slipp. Coeur d'Alene Salamander
Catalogue of American amphibians and reptiles
 , 
2002
St. Louis, Missouri
Society for the Study of Amphibians and Reptiles
 
Pages 741.1–741.3
Wilson
A. G.
Jr.
Wilson
E. M.
Groves
C. R.
Wallace
R. L.
U.S. distribution of the Coeur d'Alene Salamander (Plethodon idahoensis) Slater and Slipp
Gr. Bas. Nat.
 , 
1997
, vol. 
57
 (pg. 
359
-
362
)