Historical biogeography of western heather voles (Phenacomys intermedius) in montane systems of the Pacific Northwest

Abstract Quaternary climate fluctuations and topographical variation in the Pacific Northwest region of North America have interacted to affect the historical biogeography of biota in this region. High-elevation mammals have unique diversification patterns due to their isolation on mountaintops and potential for population growth and range expansion in lowland refugia that are available during glacial periods. We examined the phylogeographic structure, dates of lineage diversification, and historical demography of western heather vole (Phenacomys intermedius) populations across several mountain ranges in the Pacific Northwest. Our analysis of sequence variation in the mitochondrial control region using both maximum-likelihood and Bayesian methods identified 3 major geographically distinct lineages: an Oregon and California lineage, Washington lineage, and Northern and Interior lineage. Our estimate of divergence times using a Bayesian relaxed molecular-clock method revealed that diversification among these major lineages began ∼1.8 million years ago (mya) in the early Pleistocene with the split of the Oregon and California lineage followed by the split of the Washington lineage and the Northern and Interior lineage ∼1.5 mya. All 3 clades remain allopatric, suggesting that they did not share a common refugium during cold climatic intervals of the Pleistocene. Further diversification within each major clade occurred in the middle Pleistocene when populations in isolated mountain ranges became distinct. Demographic estimates from Bayesian skyline plots indicate that each of the 3 major clades has experienced population decline since the early Holocene, possibly due to the redistribution of populations into higher-elevation habitats that became restricted to mountaintops following continental and alpine deglaciation.


E 2010 American Society of Mammalogists
The interaction of historical climate changes with topographical variation can produce both intermittent lowland biogeographic refugia and high-altitude isolation of montane biota. The Pacific Northwest region of North America possesses a heterogeneous landscape that transects major mountain ranges and intervening lowland valleys and plateaus. Throughout the glacial cycles of the Pleistocene, between 2.58 million years ago (mya- Gibbard et al. 2010) and 11,700 years ago (Walker et al. 2009), the Pacific Northwest was exposed to a succession of about 20 glacial periods characterized by expansion and retraction of continental ice sheets and alpine glaciers. However, the timing of peak expansion of alpine glaciers varied across major mountain ranges (Porter et al. 1983;Thackray 2008) resulting in dynamic changes in the location of lowland ice-free zones that served as biotic refugia. Furthermore, periods of climate warming forced the redistribution of lowland biota into alpine zones that were isolated from one another by the topography of the Pacific Northwest.
Our understanding of the responses of North American plants and animals to Quaternary climate history has been enhanced by studies of genetic variation of populations across geographic space (Avise 2009). The genetic signatures of populations that resulted from historical vicariance and dispersal events developed through the interplay of 2 important evolutionary forces: genetic drift and mutation. When populations become isolated and small, genetic drift acts to reduce genetic diversity and thereby to promote differentiation among populations. By comparing both the genetic diversity of contemporary populations and the genealogical relationships among populations over geographic space, investigations of phylogeography can provide clues as to where and when isolation events and dispersal routes occurred.
Climate-mediated impacts on the structuring of genetic variation of mammal populations in the Pacific Northwest have been investigated in several low-elevation species (Arbogast and Kenagy 2001;Miller et al. 2006;Yang and Kenagy 2009;Zheng et al. 2003). One common pattern that has emerged from these studies is the apparent diversification of lineages during full glacial periods due to climate-driven shifts of populations into disjunct refugia. Another pattern is the major latitudinal shift in distribution and accompanying population expansion that follows the retraction of continental ice sheets (Lessa et al. 2003). However, mammals currently restricted to high-altitude environments experienced major elevation shifts due to these climate changes, which has led to different predictions on the outcomes of their genetic structuring and demographic changes. The expansion of alpine glaciers during cold climatic periods will force downslope shifts in populations from disjunct mountain ranges and possibly allow them to merge in common lowland refugia with greater habitat availability, which will lead to increases in population size. Moreover, upslope shift of populations into restricted mountaintops during warm periods will lead to smaller population sizes and facilitate significant lineage divergences due to founder effects (Wright 1942). Recent examples of the response of alpine specialists to these environmental scenarios demonstrate patterns that are consistent with these predictions (Carstens et al. 2005;Galbreath et al. 2009;Knowles 2001;Spaeth et al. 2009).
The western heather vole (Phenacomys intermedius) is distributed sparsely in alpine and subalpine habitat across several mountain ranges in the Pacific Northwest, including the Coastal Mountains of Canada, Rockies, Cascades, Olympics, and the Sierra Nevada (McAllister and Hoffmann 1988). The earliest known fossil records of this genus date back to the Pliocene/Pleistocene boundary ,2.5 (mya) and were discovered in both Siberia and Alaska (Repenning et al. 1987). Subsequently the purported ancestors of modern P. intermedius spread southward into Canada and the conterminous United States in the early Pleistocene (Repenning 2001). Several late-Pleistocene and early-Holocene fossils of Phenacomys have been found in wide-ranging lowland areas outside their contemporary distribution (Diveley 1999;Emslie 1986;Graham and Mead 1987;Grayson 1981;Karrow et al. 1995;Lundelius et al. 1983;Lyman 2008;Rensberger and Barnosky 1993;Rogers et al. 2000;Walker 1987; Fig. 1). Based on palynological studies, habitat likely to be suitable for P. intermedius was widespread throughout the Pacific Northwest lowlands during full glacial conditions (Barnosky 1985;Daniels et al. 2005;Thompson and Anderson 2000;Worona and Whitlock 1995). However, it is uncertain whether refugial locations of P. intermedius represented 1 large common refugium or disjunct refugia separated by physical barriers. Furthermore, a large portion of the contemporary distribution of P. intermedius occurs in high northern latitude areas that were glaciated by continental ice sheets during the last glacial maximum. The Cascades and Rockies represent the 2 major mountain axes that span nonglaciated and glaciated areas and might have served as important northward dispersal routes for P. intermedius following the retraction of continental ice sheets.
To assess the impacts of late-Quaternary climate change and geological features on the historical biogeography of P. intermedius across montane systems of the Pacific Northwest, we obtained samples throughout several montane regions from northern British Columbia, Canada, southward to California. Our objective was to incorporate phylogenetic inference analyses, relaxed molecular-clock dating computations, and historical demographic estimations into a phylogeographic framework to test the hypotheses that contemporary isolated populations were forced into a series of disjunct refugia separated by topographical barriers during full glacial conditions in the Pleistocene rather than into a common refugium; dispersal into northerly areas that were previously occupied by continental ice sheets occurred from multiple refugial sources; postglacial shifts of elevation into isolated Northwest. Also identified is the geographic extent of each of the 3 major clades as defined in Fig. 2. The inset shows the entire contemporary distribution of P. intermedius in dark outline (redrawn from locality data in Hall [1981], Nagorsen [2005], and Verts and Carraway [1998]), approximate margins of the Cordilleran and Laurentide Ice Sheets (,14,000 years before present) as the shaded area north of the dotted line (Siegert 2001), and major alpine glaciers as shaded areas south of the ice sheets (Porter et al. 1983). Late-Pleistocene and Holocene fossil localities outside the modern heather vole range are indicated by black circles. mountain ranges resulted in lineage divergence; and major lineages are experiencing declines in population since the end of the Pleistocene due to their upslope range shift and reoccupation of restricted montane regions.

MATERIALS AND METHODS
Sampling.-Our analysis is based on 89 samples of P. intermedius from 18 localities representing the northern and western portions of the geographic range, with a special emphasis on isolated montane populations of the Pacific Northwest (Appendix I; Fig. 1). The samples consisted of 36 frozen tissues from the liver of specimens collected between 1999 and 2008, which were deposited in the Burke Museum, University of Washington (UWBM), and the Museum of Vertebrate Zoology, University of California (MVZ). The remaining 53 samples were obtained from skins of museum study specimens collected between 1925 and 2000, which were deposited in the UWBM; the Slater Museum of Natural History, University of Puget Sound (PSM); and the Department of Fisheries and Wildlife, Oregon State University (OSUFW; Appendix I). Research was conducted in adherence with guidelines of the American Society of Mammalogists (Gannon et al. 2007).
Laboratory methods.-We extracted whole genomic DNA from fresh liver or spleen tissue using the prescribed protocol of DNeasy Tissue Kit (Qiagen, Valencia, California). For samples from museum study skins genomic DNA was extracted following a protocol developed by Mullen and Hoekstra (2008), which included an ethanol wash every 3 h for 24 h to remove salts and polymerase chain reaction (PCR) inhibitors that may have been used in preserving museum skins. For these samples we also undertook several steps to avoid contamination. First, we used a new sterilized razor blade for the removal of snippets from each specimen. Second, we performed ethanol washes in a separate room from where PCR amplifications were performed to avoid contamination from PCR amplicons. Next, we divided our samples into smaller subsets and staggered the dates of extractions and PCR amplifications for each subset. We also checked for contamination by checking for PCR products from negative-control extractions and negative-control PCR master mixes with gel electrophoresis. Finally, we repeated DNA extractions for 5 random individuals using skin snippets from another part of the study specimen and compared their sequences with original samples.
We amplified a 325-base pair (bp) sequence of the mitochondrial control region using the PCR with primers RTVL16000 (59-GTCAACACCCAAAGCTGACA-39) and RtvInt1r (59-GTTGGTTTCACGGAGGATGG-39- Miller et al. 2006). PCR mixtures totaled 25 ml and followed slightly different protocols for frozen tissue and museum study skin samples. For the frozen tissue samples each reaction mixture contained 14.9 ml of nuclease-free H 2 O, 2.5 ml of 10X PCR buffer, 1 ml of 25 mM MgCl 2 , 0.2 mM of each deoxynucleoside triphosphate, 2.5 ml of each 1 X primer, 0.125 ml of JumpStart Taq DNA polymerase (Sigma, St. Louis, Missouri), and 1 ml of genomic DNA. We performed PCR amplifications in a GeneAmp PCR System 9700 thermocycler (Applied Biosystems, Foster City, California). The cycle conditions included a denaturing step at 94uC (10 min), followed by 30 cycles (1 min at 94uC, 30 s at 50uC, and 1 min at 72uC) with a final extension period of 5 min at 72uC. Because DNA from study skins was presumed to be somewhat degraded, we modified the PCR master mix to improve PCR amplification by adding 2.5 ml of 10X bovine serum albumin, 2 ml of 25 mM MgCl 2 , 0.5 ml of JumpStart Taq DNA polymerase, and 9.9 ml of genomic DNA, and by not adding the 14.9 ml of nucleasefree H 2 O. We also increased the number of amplification cycles to 40 during PCR. We treated all PCR products with ExoSapIT (USB Corp., Cleveland, Ohio) to remove unincorporated nucleotides and primers. Sequencing reactions for PCR products of each sample were performed separately for each forward and reverse primer. Sequencing reaction mixtures of 10 ml for each sample included 1 ml of PCR range. The phylogeny inferred with maximum-likelihood methods yielded an essentially identical topology. Numbers above nodes represent maximum-likelihood bootstrap percentages . 80 and those below represent Bayesian posterior probabilities . 0.85. Sampling localities (A-R) and haplotype numbers (1-58), as shown in Appendix I, are indicated at branch tips. product, 7.5 ml of nuclease-free H 2 O, 0.5 ml of the Big Dye Terminator Cycle Sequencing Ready Reaction Mix (Applied Biosystems), 0.5 ml of the Big Dye 5X Sequencing Buffer, and 0.5 ml of primer. Samples were run on an ABI 3100 Genetic Analyzer (Applied Biosystems) for sequencing of nucleotides. Sequences were aligned using SEQUENCHER 4.6 (Gene Codes Corp., Ann Arbor, Michigan) and deposited in GenBank under accession numbers GQ903588-GQ903676.
Phylogenetic analyses.-We used both maximum-likelihood and Bayesian inference methods to reconstruct phylogenetic relationships of haplotypes. We determined that the Hasegawa-Kishino-Yano (Hasegawa et al. 1985) model of nucleotide substitution was the best-fitting model for our phylogenetic analyses by examining Akaike's information criteria corrected for small sample size differences (DAIC c ) and Akaike weights in MODELTEST (version 3.7-Posada and Crandall 1998). Other parameters estimated were as follows: base frequencies (A 5 0.3547, C 5 0.0722, G 5 0.2312, T 5 0.3419); proportion of invariable sites (I) 5 0.6089; and gamma shape parameter 5 0.5954. We reconstructed the maximum-likelihood phylogeny using PHYML 2.4 (Guindon and Gascuel 2003) with 1,000 bootstrap replicates to evaluate nodal support for phylogenetic clades (Felsenstein 1985). The Bayesian inference phylogenetic reconstructions were performed with MRBAYES 3.1 (Ronquist and Huelsenbeck 2003) using Markov chain Monte Carlo sampling. Runs were set for 5,000,000 generations, with 1 cold and 3 heated chains and with trees sampled every 100 generations. The 1st 25% of sampled trees were discarded as burn-in after visual inspection revealed that these chain samples had not converged on stable likelihoods (stationarity). Outgroup sequences were obtained from GenBank from 1 Arborimus albipes (accession DQ198850) and 1 Dicrostonyx richardsoni (AF192739).
Divergence dating.-We estimated the time to the most recent common ancestor of all haplotypes and several monophyletic clades using BEAST 1.5.2 . BEAST brings together several complementary evolutionary models (i.e., substitution models, insertiondeletion models, demographic models, tree-shape priors, relaxed clock models, and node-calibration models) into a single coherent framework for evolutionary inference. Based on our MODELTEST results, we used an Hasegawa-Kishino-Yano substitution model with a gamma rate distribution and proportion of invariant sites. We also determined from a likelihood ratio test (Felsenstein 1981) that a relaxed molecular-clock model with an uncorrelated lognormal tree prior best accounted for the varying rates of evolution across lineages in our data set. To calibrate the tree we added several outgroup taxa belonging to higher taxonomical levels and used estimated fossil dates to constrain node heights and topological features. These node constraints included (F1) the earliest fossil representing the subfamily Arvicolinae, which was found 5.5 mya (95% confidence interval [95% CI] 5 5.0-6.0 mya) in East Asia from the early Pliocene (Chaline et al. 1999;Jin and Zhang 2005); and (F2) the earliest fossil representing the tribe Phenacomyini, found in the early Pleistocene 2.4 mya (95% CI 5 2.1-2.6 mya- Repenning et al. 1987;Fig. 3) We obtained the following genetic sequences of outgroup taxa representing these higher-order groups from GenBank (with accession numbers): Peromyscus maniculatus (GQ860675) as the representative outgroup taxon from the family Cricetidae; Eothenomys chinensis (FJ483847), Microtus mexicanus (AF251260), Myodes glareolus (Y07543), and Lemmus sibiricus (AF355450) as outgroup taxa belonging to the subfamily Arvicolinae; and Arborimus albipes (DQ198850) as the representative outgroup taxon belonging to the supergenus Phenacomyini. In addition, we sequenced 1 specimen of Microtus montanus (GU394082) to include with the Arvicolinae group. We used a Yule tree prior for this analysis because our data set included several species-level taxa. Therefore, we used a reduced data set for our ingroup sequences that included only a few individuals from each of the 3 major clades as determined by our maximum-likelihood and Bayesian inference phylogenetic analyses. This data set included 4 individuals from the Northern and Interior clade and 3 each from the Washington clade and the Oregon and California clade (see letters a-j in Appendix I). We also included 1 sequence from subclades identified in the phylogenetic analysis (Wallowa, Olympic, and the Sierra Nevada) as part of the reduced data set. The analysis was run 5 times with different random seeds, each for 10 million generations with samples logged every 10,000 generations, and also run for 50 million generations to test for convergence. The initial 10% of samples was discarded as burn-in after visual inspection of the trace revealed that these samples had not reached stationarity. We combined results from each run in LOGCOMBINER 1.5.2 from the BEAST package and analyzed them in TRACER 1.5 . Effective sample size values for all parameters exceeded 200.
Molecular diversity and demographic analyses.-We computed haplotype (h) and mean nucleotide (p) diversities (Nei 1987) for each clade using ARLEQUIN (Excoffier et al. 2005) for comparisons of genetic diversity. We also examined the demographic history of each major clade using Fu's F s test for neutrality (Fu 1997). Large and negative F s values suggest an excess of rare alleles, which can be an indication of population growth or selection. We also investigated patterns of historical demographic expansion by examining the distribution of pairwise differences (mismatch distributions) for each clade (Rogers and Harpending 1992) using ARLEQUIN. In general, populations undergoing recent and sudden expansion exhibit a Poisson-shaped mismatch distribution, and populations in equilibrium tend to have ragged distributions (Slatkin and Hudson 1991). We assessed the statistical significance of these distributions with sum of squared distances and Harpending's raggedness index (Harpending 1994).
We also explored demographic history of each major lineage using Bayesian skyline plots (Drummond et al. 2005) in BEAST. This coalescent-based inference method uses a Markov chain Monte Carlo sampling procedure with gene sequence data to estimate a posterior distribution of effective population size through time. We used divergence date estimates of time to the most recent common ancestor from our divergence analysis as root age priors for each of the major lineages so that we could model changes in effective population size against a timescale of years rather than expected substitutions per site. We also used the Hasegawa-Kishino-Yano substitution model with a gamma distribution and proportion of invariant sites, along with a relaxed clock model with an uncorrelated lognormal prior. We forced coalescent intervals into 5 groups for each lineage because test runs for higher group numbers did not converge as well. We used default settings for the remaining model parameter priors. We performed 5 runs with different random seeds, sampled every 10,000 generations, and discarded 10% of the samples as burnin for each major lineage. We varied the length of runs for each lineage with the Northern and Interior clade set at 100 million generations, the Washington clade set at 80 million generations, and the Oregon and California clade set at 40 million generations. Effective sample size values for all parameters in each run exceeded 200. We combined results for each lineage in LOGCOMBINER and produced skyline plots in TRACER.

RESULTS
Phylogenetic inferences and estimates of divergence dates.-Both the maximum-likelihood and Bayesian phylogenetic inferences reveal 3 major clades (Fig. 2), each associated with distinct geographic areas: a Northern and Interior clade represented by populations in the Northern Rockies, Wallowa Mountains, Cascades of northern Washington and southern British Columbia, and the northern interior mountains of British Columbia (Fig. 1,  Within each major clade we found 1 well-supported subclade associated with a geographically isolated mountain range, as follows. Within the Northern and Interior clade the isolated Wallowa Mountains of eastern Oregon (J) contain a lineage with high Bayesian nodal support (Fig. 2). Within the Washington clade the lineage in the isolated Olympic Mountains (N) shows strong nodal support. Within the Oregon and California clade the disjunct population in the Sierra Nevada (R) shows strong nodal support.
The estimation of time to the most recent common ancestor ( Fig. 3; Table 1) indicates that the tribe Phenacomyini diverged from its sister lineage in the subfamily Arvicolinae ,5.37 mya. The estimated mean time to the most recent common ancestor for the split of the genera Phenacomys and Arborimus was ,2.45 mya. This analysis estimated an early-Pleistocene origin of ,1.8 mya for the Oregon and California clade and an early-Pleistocene divergence of ,1.54 mya between the Washington clade and Northern and Interior clade. Sublineage diversification of isolated mountain populations (Wallowa, Olympic, and Sierra Nevada) occurred in the middle Pleistocene within each of the 3 major clades.
Molecular diversity and demographic analyses.-All 3 clades had high haplotype diversity and low nucleotide diversity (Table 2), and no haplotype was shared among the 3 regions occupied by the major clades. Only 4 haplotypes were found at .1 sampling locality. The significantly negative Fu's F s for the Northern and Interior clade suggests that this clade experienced recent rapid expansion. This is corroborated by the mismatch distributions, which were unimodal for this clade and did not possess a statistically significant Harpending's raggedness index or sum of squared deviation. The multimodal mismatch distributions and nonsignificant Fu's F s values for the other 2 major clades suggest population stability; however, this conflicts with the nonsignificant Harpending's raggedness index and sum of squared deviation of the mismatch distribution for both clades. This lack of correspondence could be due to the more conservative test based on pairwise sequence differences (mismatch distribution) as opposed to the Fu's F s (Ramos-Onsins and Rozas 2002).
The median estimates from the Bayesian skyline plots indicate that the populations of all 3 major lineages have rapidly declined following a long period of relative stability. This decline occurred after the last glacial maximum and at the beginning of the warming period in the early Holocene (Fig. 4).

DISCUSSION
The interaction between climate-mediated shifts in geographic distribution during the Pleistocene and variation in topography associated with regional occurrence of mountains has played a strong role in the phylogeographic structure of western heather voles in the Pacific Northwest. The initial southward range expansion of Phenacomys along the Pacific Coastal region into the conterminous United States (,1.8 mya) from its Arctic origins in Alaska and northeastern Siberia occurred during a warm interglacial period (Repenning 1990). Repeated intervals of glaciation through the remainder of the Pleistocene could have resulted in range expansions and contractions and subsequent secondary contact between diverging lineages in lowland refugia or in recolonized highelevation regions. However, the geographic distinctiveness of the 3 contemporary lineages suggests that barriers were persistent throughout these repeated cycles of climate change and did not allow secondary contact to occur.
We can attempt to identify geographic barriers that prevented secondary contact between the 3 major lineages. The Columbia River Gorge lies between the Oregon and California clade and the Washington clade and appears to be an obvious barrier. Less obvious is why populations belonging to the Northern and Interior clade did not mix with populations belonging to the other clades in lowland refugia or in recolonized alpine sites. Fossil records of Phenacomys indicate that populations existed in lowland regions of both western and eastern Washington during the late Pleistocene and early Holocene (Karrow et al. 1995;Lyman 2008;Rensberger and Barnosky 1993 ; Fig. 1). Palynological studies also indicate the presence of habitats (parkland tundra and steppe) during full glacial conditions that might have been able to support widespread populations of Phenacomys in Washington (Barnosky 1985;Heusser 1983). The contemporary presence of only Washington clade populations in the southern portion of the Washington Cascades suggests that upslope movement into these mountains occurred only from western refugial populations. The Columbia River runs along the eastern side of the Cascades in Washington and might have prevented recolonization of this portion of the Cascades by refugial populations from eastern Washington. However, the contemporary presence of only Northern and Interior clade populations in the North Cascades region of northern Washington and southern British Columbia suggests that dispersal into this region occurred from an eastern refugial source, possibly by populations from eastern Washington. The phylogeographic break between the Northern and Interior clade and the Washington clade in the Washington Cascades occurs over a distance of only 50 km and might have been maintained by increasing restriction of dispersal during recent warming in the Holocene combined with short neoglacial advances of montane glaciers in the Glacier Peaks region (Burke and Birkeland 1983;Thompson Davis et al. 2009). For the Oregon Cascades and Sierra Nevada populations no obvious barriers are apparent in the intervening Great Basin region that historically FIG. 3.-Bayesian tree with divergence dates from the relaxed molecular-clock analysis of haplotypes of Phenacomys intermedius. Circled numbers refer to estimated divergence dates, and gray horizontal bars indicate 95% credibility intervals (Table 1). Specimen numbers of Phenacomys indicated by letters a-j are provided in Appendix I, and outgroup sequence numbers and fossil calibration dates, indicated by F1 and F2, are provided in the ''Materials and Methods.'' All nodes possessed posterior probability support . 0.95 except for those indicated by gray stars.   (Diveley 1999;Grayson 1981;Thompson and Anderson 2000). The maintenance of geographic distinctiveness among major lineages of Phenacomys in the Pacific Northwest could be a product of the high degree of niche conservatism in this species that resulted from its Arctic origins  and restriction in modern and historical periods to cold environments. This could have played a role in the ability of Phenacomys to disperse across topographical barriers and persist through rapid environmental change. In contrast, Microtus longicaudus, another high-elevation arvicoline, possesses a wider niche breadth and was able to disperse across prominent geographic features during postglacial movements, which allowed its divergent lineages to come into secondary contact (Spaeth et al. 2009).
Postglacial routes for range expansion of western heather voles into northern areas previously occupied by continental ice sheets appear to have occurred along an interior axis rather than a coastal or Cascade Mountain axis. Examination of our data shows minimal northward expansion of Washington clade populations in comparison to Northern and Interior clade populations. In fact, only populations belonging to the Northern and Interior clade reoccupied areas previously overlaid by continental ice sheets, including the North Cascades. This spatial distribution suggests an interior dispersal route along or near the Rocky Mountain axis. The earlier retreat of the Laurentide Ice Sheet, as opposed to the Cordilleran Ice Sheet, during late-Pleistocene deglaciation (Siegert 2001) could have played a role in allowing northward expansion by interior refugial populations and the prevention of major northward expansion by the more-coastal Washington refugial populations. This pattern of asynchronous northward expansion by divergent lineages has been found in several low-elevation boreal mammals (Arbogast 1999;Demboski et al. 1999) and in another high-elevation arvicoline with similar habitat associations, the water vole (Microtus richardsoni- Carstens et al. 2005). The American pika (Ochotona princeps), however, displays a more symmetrical northward expansion between Cascade and Rocky Mountain lineages despite occupying similar high-elevation mountains in western North America (Galbreath et al. 2009). This demonstrates that climate-mediated movements did not necessarily occur in all assemblages of mammals but were likely more individualistic and dependent on each species' ecological breadth (Grayson 2006;Lawlor 1998).
Within the smaller geographic scale of each major clade we also found significant substructuring of populations of P. intermedius associated with shifts of elevation into isolated montane systems. The results from our divergence dating analyses indicate that the initial divergence of these sublineages began during the middle Pleistocene. As with the case for the 3 major clades, these sublineages persisted in distinct refugia despite repeated range expansions and contractions during glacial cycles associated with the remainder of the Pleistocene. Barriers that prevented mixing of lineages might have resulted from various prominent geographic features. The divergence between the Olympic Mountain and Washington Cascade populations possibly was reinforced by repeated advances of the Puget Lobe of the Cordilleran Ice Sheet during full glacial periods (Easterbrook 1986). The divergence between populations in the Wallowa Mountains and the rest of the populations in the Northern and Interior clade likely was due to the presence of the Snake River, whereas the divergence between populations in the Oregon Cascades and the Sierra Nevada could have been reinforced mostly by postglacial changes of vegetation in intervening regions from cool subalpine parkland to present-day mesic forests (Briles et al. 2005;Daniels et al. 2005;Hakala and Adam 2004;West 2004). Furthermore, lineage divergence also appears to be dependent on sufficient time for genetic drift to sort shared ancestral variation. None of the populations that recolonized northerly areas, which were covered by continental ice sheets, were shown to be divergent.
Results from our Bayesian skyline plots reveal population declines in all 3 major lineages of heather voles since the beginning of the Holocene. This pattern is also shared by pikas (Galbreath et al. 2009) and likely reflects the reduced amount of habitat in restricted mountain environments following postglacial upslope movements from more expansive lowland refugia. This pattern contrasts with those of low-elevation temperate mammals, which commonly show signatures of postglacial demographic expansion resulting from geographic expansion into large areas of recently available habitat (Lessa et al. 2003;Runck and Cook 2005). This pattern also might help explain why some of the more southerly interior populations, such as those that existed in the Great Basin Desert, are no longer extant. Nevertheless, some of our demographic analyses indicate that the Northern and Interior clade experienced rapid population growth. Although these results conflict with the outcome of the Bayesian Skyline Plots, they are consistent with the observation that this clade is the only major lineage to experience major geographic expansion since the end of the Pleistocene and the only case of substantial long-distance haplotype sharing (northern British Columbia and Montana; Fig. 2).
We have shown that the combination of major climatic fluctuations and topographical variation found within a regional montane system has had an important impact on the genetic variation of an alpine mammal, the western heather vole, in the Pacific Northwest. P. intermedius supports some of our predictions for high-elevation species, such as postglacial decline in population size and differentiation of populations restricted to isolated mountain ranges. However, in some ways P. intermedius showed similar patterns to those of low-elevation species; for example, the separation of refugial populations into distinct geographical regions and the asynchronous northward expansion by divergent lineages into areas previously glaciated by continental ice sheets. Studying the extreme environments and restricted ranges of alpine organisms is useful for understanding the tolerances of temperate mammals to strong environmental heterogeneity.