- Research article
- Open Access
Diversification of the Alpine Chipmunk, Tamias alpinus,an alpine endemic of the Sierra Nevada, California
© Rubidge et al.; licensee BioMed Central Ltd. 2014
- Received: 12 September 2013
- Accepted: 11 February 2014
- Published: 23 February 2014
The glaciation cycles that occurred throughout the Pleistocene in western North America caused frequent shifts in species’ ranges with important implications for models of species divergence. For example, long periods of allopatry during species’ range contractions allowed for the accumulation of differences between separated populations promoting lineage divergence. In contrast, range expansions during interglacial periods may have had homogenizing effects via increased gene flow following secondary contact. These range dynamics are particularly pronounced in the Sierra Nevada, California, given the complex topography and climatic history of the area, thus providing a natural laboratory to examine evolutionary processes that have led to the diversity patterns observed today.
Here we examined the role of late Pleistocene climate fluctuations on the divergence of the Sierra Nevada endemic Alpine Chipmunk (Tamias alpinus) from its sister taxon, western populations of the Least Chipmunk (T. minimus) from the Great Basin. We used one mitochondrial gene (cytochrome b) and 14 microsatellite loci to examine the evolutionary relationship between these species. Mitochondrial sequence data revealed that T. alpinus and T. minimus populations share mitochondrial haplotypes with no overall geneaological separation, and that diversity at this locus is better explained by geography than by species’ boundaries. In contrast, the microsatellite analysis showed that populations of the same species are more similar to each other than they are to members of the other species. Similarly, a morphological analysis of voucher specimens confirmed known differences in morphological characters between species providing no evidence of recent hybridization. Coalescent analysis of the divergence history indicated a late Pleistocene splitting time (~450 ka) and subsequent, though limited, gene flow between the two lineages.
Our results suggest that the two species are distinct and there is no contemporary introgression along their geographic boundary. The divergence of T. alpinus during this time period provides additional evidence that Pleistocene glacial cycles played an important role in diversification of species in Sierra Nevada and North America in general.
- Microsatellite Locus
- Markov Chain Monte Carlo
- Effective Population Size
- Great Basin
- Secondary Contact
Understanding processes that promote and maintain biodiversity is a key goal of evolutionary biology. Divergent natural selection resulting from resource heterogeneity and competitive interactions can drive population divergence and speciation [1, 2]. Nonadaptive divergence, operating via genetic drift due to isolation and founder effects, may also play a significant role in generating patterns of species diversity. Furthermore, hybridization (or reticulate evolution) during and subsequent to speciation can add novel genetic diversity to diverging lineages and affect the course of adaptive divergence [3–5]. The cyclical Pleistocene glacial and interglacial episodes have shaped the genetic architecture of taxa across the globe, as retraction to refugia facilitated the formation of distinct evolutionary lineages within species [6, 7] through both vicariant processes and shifting the spatial distribution and extent of ecological gradients. Recurrent fragmentation and expansion provides opportunities for initial divergence, repeated secondary contact, hybridization, and demographic fluctuation .
Range dynamics and the complex geological and climatic history of the Sierra Nevada, California has shaped diversity patterns of the central-western United States over the time span of most, if not all, extant species. The late Pleistocene in particular was a time of drastic climatic fluctuations in the Sierra Nevada  and the adjacent Great Basin and Mohave biomes  and thus a period linked to both intra- and interspecific diversification in many taxa (e.g., [11–14]).
The physical setting: Sierra Nevada and the Great Basin
The Sierra Nevada and the adjacent Great Basin is a biologically diverse region with a rich glacial history [10, 15–17]. The Sierra Nevada are a narrow, elongated, topographically complex, high, and relatively young mountain range that spans about 640 km from north to south and 110 km from west to east, and contains the highest peak in the continental United States (Mt. Whitney, 4,421 m). Given the steep elevational gradient on the east side of the Sierra Nevada, the mountain range casts a major rain shadow affecting the climate and ecology of the adjacent central Great Basin of the intermountain west. The Great Basin consists of over 10,500 km2 of valleys, basins, lakes and mountain ranges with extreme elevational relief throughout the region. It contains the lowest point in North America (-86 m, in Death Valley, California) as well as one of the highest (4342 m, White Mountain Peak, CA) . The complexity of small mammal distributions in relation to the high environmental diversity in the Great Basin, have made this region a focus for biogeographic research (e.g., [17–20]).
The interactions between Sierra Nevada with the Great Basin and Mojave Desert biomes during the cyclical glaciations of the Pleistocene created areas of recurrent isolation and secondary contact within and between species. The Sierra Nevada has experienced at least six cyclical episodes of glacial expansion and retreat, both regional and local in extent, beginning in the early Pleistocene (1.65 ± 0.7 mybp), with the most recent advance around 3.5 kbp , as well as three or more neoglacial advances ending with termination of the Little Ice Age, which spanned 1350 to 1850 AD . These episodes are intertwined with a similarly complex history of volcanism along the eastern flank of the range over the same time period . Few glaciers remain in the Sierra Nevada today, with most having shrunk by greater than half in mass and aerial extent within the latter half of the 20th century . This complex glacial history and associated range dynamics and steep environmental gradients in the Sierra Nevada and the adjacent Great Basin provide a natural laboratory to examine species histories including the combined effects of range fluctuations, oscillations between isolation and introgression, and adaptive divergence across environmental gradients (e.g., [25, 26]).
Chipmunk diversity in western North America
The biogeographic history, radiation, and evolutionary relationships of western North American chipmunks are complex and include several instances of historical introgression among species [27–32]. Chipmunk diversity is centered in western North America with 23 species [28, 33] and diversification has been linked to shifting ranges resulting from climatic cycles and subsequent shifts in habitat preference resulting from interspecific competition and niche partitioning over elevational gradients [34–36]. In this study, we examine the evolutionary history of the Alpine Chipmunk, Tamias alpinus, a narrow high-elevation endemic, relative to its widespread sister species, the Least Chipmunk, Tamias minimus.
The divergence history of T. alpinus is not well understood; however, recent work has shown that T. alpinus and T. minimus are paraphyletic . T. minimus represents a species complex and this study focuses on the western segment of the species, geographically adjacent to the range of T. alpinus. Previous work by Reid et al.  included geographically (and taxonomically) disparate sequences from T. minimus individuals from California, Nevada, Utah, Wyoming and Washington and clearly showed that T. alpinus is nested with western segments of T. minimus scrutator. Here, we focus on regions of allo/parapatry in the Sierra Nevada where these two species come into close proximity to better understand the evolutionary history of T. alpinus.
Our objective is to examine the relationship between T. alpinus and T. minimus, with the goal of gaining a better understanding of the evolutionary history of T. alpinus. More specifically, we focus on two questions 1) did T. alpinus diverge from T. minimus in association with late-Pleistocene glacial dynamics and 2) is there evidence for contemporary or historical introgression as reported in other species-pairs of chipmunks [30, 31]. To address these questions we examine morphological characters and genetic variation and population structure at one mitochondrial gene and 14 microsatellite loci. Contrasting patterns of nucleotide variation in the mitochondrial genome with patterns of genetic variation at microsatellite markers gives us a picture of both historical and contemporary processes respectively. And although differences in morphological characters between species are well defined (see below), testing for morphological intermediacy in individuals in adjacent versus distant populations will help distinguish between genetic similarity due to recent speciation or similarity as a result of secondary contact and hybridization.
T. alpinus is geographically restricted to the high elevations of the central to southern Sierra Nevada (Figure 1). Based on historical records and our own surveys, T. alpinus and T. minimus are allopatric throughout our study area. South of Yosemite, there are reported areas of sympatry in the central Sierra Nevada northwest of Bishop (D. Guiliani, pers. comm.), but none have been confirmed despite our targeted field surveys between 2009-2013. Collections made in 1911 in the southern Sierras recorded both species at the same locality (Little Brush Meadow, Tulare Co.), but apparently at different elevations (9750 ft for T. minimus and 10,000 ft for T. alpinus; field notes of the collectors, H. A. Carr and W. P. Taylor, MVZ archives). This difference in elevation signals distinctly different habitats, where lodgepole pine (Pinus contorta) with a sagebrush understory at the lower elevation is replaced by scattered foxtail pine (Pinus balfouriana) and arctic-alpine forbs near treeline. These habitats characterize the ecological distribution of the two species in other areas where their ranges come into close proximity.
The two species differ in morphology and habitat preferences. T. minimus, among other characters, is smaller in body mass, has a longer tail, shorter ears and darker coloration than T. alpinus and has a different bacular (penis bone) morphology [37, 40]; more specific details of these differences are summarized in Additional file 1: Table S1. In California, T. minimus is found in arid sagebrush habitat that ranges in elevation from 1500 m to above 3000 m in the Sierra Nevada and mountains to the immediate east (e.g., Sweetwater, White, and Inyo ranges). T. alpinus is restricted to the alpine zone of the Sierra Nevada at and above tree-line (2950 to 4100 m) where it occupies open granite habitat, meadow edges and talus slopes [39, 41].
Study site and samples
The study area is the central to southern Sierra Nevada, CA, USA, which includes the entire known range of T. alpinus and that of T. minimus to the immediate east and north (Figure 1). A total of 341 chipmunks were included in this study. The majority of samples were collected between 2003-2009, however 14 samples were taken from museum skins that were collected between 1911-1916 and housed in the Museum of Vertebrate Zoology (MVZ) at the University of California, Berkeley. For the samples collected between 2003-2009, we live-trapped animals using Sherman traps at 62 locations between 2003-2009 (Figure 1). We used non-lethal sampling (ear clips) and collected vouchered specimens, including liver samples, now catalogued in the MVZ (see Additional file 2: Table S4 for MVZ catalogue numbers and locations). Chipmunks collected in areas of potential sympatry were identified to species in the field based on distinct morphological differences including body size, pelage color, and ear and tail length. More rigorous morphological measurements on collected specimens were conducted in the lab to better document species-specific morphological differences. Sample collection was approved by the Animal Care and Use Committee (Protocol #R304-0509) at the University of California, Berkeley.
Morphological methodology and analyses
Despite well-characterized differences between species, these two taxa share a close genetic legacy, and in order to distinguish between the possibility that this similarity resulted from recent common ancestry versus reticulation subsequent to divergence (e.g., [30, 31, 42]), we examined the morphology of specimens taken from geographically adjacent versus distant localities to test for morphological intermediacy, using the same geographic groupings in the molecular comparison (Figure 1). We examined three data sets separately. The first containing external body metrics that have been used to separate the two species including: a) color and color pattern (b) tail length and bushiness; (c) ear length (N = 36 and 129 for T. alpinus and T. minimus respectively). The second dataset was based on comparing bacular dimensions (N = 15 and 33; shaft length, mid-shaft height, tip height, tip angle, and keel breadth following ). The third and final dataset examined craniodental features measured from preserved skulls (N = 161 and 159). We used separate conical variate analyses (CVA; JMP 5.1.1 statistical software) for each dataset, which generates a classification matrix of group membership, with individual posterior probabilities, based on multivariate discriminant functions.
DNA extraction, sequencing, and microsatellite genotyping
To extract DNA from the liver or ear tissue samples, we used the standard Qiagen DNAeasy kit following the manufacturer’s protocol (Qiagen). Tissue extractions were eluted in a total of 400 μl AE buffer. All extractions and PCR set-up on the 14 skin samples from museum specimens (“historical” samples) were conducted in a separate laboratory devoted to ancient DNA research. We followed the museum skin DNA extraction protocol described in Mullen and Hoekstra . After removing an approximately 3 mm × 3 mm square piece of skin from the lower lip, we placed the sample in 95% ethanol and refreshed the ethanol roughly every 3 hours over a 24-hour period to wash the sample of salts and PCR inhibitors. Following these washes, each skin sample was carefully removed of hair and shaved into smaller pieces with a scalpel and placed into a 1.5 ml locking Eppendorf tube. Between each sample, the forceps and scalpel were washed in 10% bleach, rinsed in 95% ethanol and flamed to avoid cross contamination. We extracted DNA using a Qiagen DNeasy Tissue Extraction kit with the following modifications. First we diluted the AE Buffer to 1:10 in RNAse-free H2O and warmed it to 70°C prior to elution. Second, we applied two elutions of 50 μl of warm 1:10 AE Buffer to the spin columns and allowed this elution step to incubate at room temperature for 5 min prior to the final spin. We conducted a negative extraction (sterilized forceps in extraction buffer) alongside all historical skin extractions. The negative extraction was run along with a negative PCR control in reactions to test for contamination between samples.
Mitochondrial DNA analysis
An 801 bp portion of the mitrochondrial gene, Cytochrome b (cyt b), was amplified using universal mammal primers, MVZ05 & MVZ16 . We sequenced 139 T. alpinus and 107 T. minimus samples for a total dataset of 246 sequences. For the 14 historical samples, because the DNA was more degraded we were only able to amplify a fragment less than 400 bp. Therefore, we developed three pairs of genus specific primers to amplify shorter fragments that could be pieced together to complete a 780 bp sequence for the historical DNA samples (Additional file 3: Table S2). To ensure the species-specific primers were not causing any irregularities in the sequences, we also used these primers on five modern DNA samples and compared the results with the sequences using the MVZ universal primer pair. The thermal cycler conditions for the mitochondrial PCRs were as follows: 94°C for 2mins, 35 cycles of 94°C for 30s, 47-50°C for 30s, 72°C for 60s and then a final extension at 72°C for 5mins. Historical samples were sequenced in both the forward and reverse direction and PCR’d and sequenced at least twice to ensure repeatability of resulting sequence. Amplicons were sequenced on an ABI 3730 Capillary Sequencer (Applied Biosystems, Inc.). Resulting sequences were edited and aligned using Sequencher 4.8 (Gene Codes Corp.).
Mitochondrial data analyses
We used two approaches to estimate genealogical relationships. First, to estimate the phylogenetic relationships of haplotypes, we used the Bayesian approach implemented in MRBAYES 3.1.2 . The best-fit model of nucleotide change was estimated using Akaike Information Criterion as implemented in jModeltest . The model of sequence evolution ranked highest by AIC for the dataset was the Tamura-Nei model (TrN + I + Γ) but because TrN is not an option in MrBayes, and this model is a special case of the general time reversible model (GTR) we used GTR + I + Γ. We ran four MCMC chains for 3,000,000 generations with trees sampled every 300 generations. We assessed convergence by examining the standard deviation of split frequencies, which were <0.01 after 3 × 106 generations. A burn-in period of 105 was discarded prior to calculating the consensus tree. Three individuals of the Panamint chipmunk (Tamias panamintinus) were used as an outgroup to T. alpinus and T. minimus. Traditional phylogeny reconstruction approaches such as described above, however, make several assumptions that make them inaccurate at the population level. For example, they assume ancestral haplotypes are no longer present in the population. Therefore, our second approach was to use haplotype networks to estimate the genealogical relationship using the statistical parsimony approach  as implemented in the program TCS 2.1 .
Genetic variation among sequences within species was quantified as haplotype diversity (hd), and nucleotide diversity (θπ. θS). To quantify mtDNA differentiation between species and/or populations we calculated the average and the net number of nucleotide substitutions per site (Dxy and Da, ). To visualize divergence patterns we clustered individuals by geography and species and used Dxy to produce a neighbour-joining tree in MEGA version 4 . All diversity and distance calculations were estimated in program DNAsp version 5 . We used an Analysis of Molecular Variance (AMOVA) implemented in the program ARLEQUIN  to examine the population structure of sequence diversity. F-statistic analogues (ϕ) were calculated to estimate the differentiation among groups (ϕCT) among populations within groups (ϕSC) and within populations (ϕST). Populations were grouped according to their species designation and sampling locality (Figure 1). We tested the statistical significance of the AMOVA with 10000 permutations and corrected the p-values associated with the ϕ values using Bonferroni correction for multiple tests. To test for historical population expansion or contraction in each species we calculated Tajima’s D  and Fu’s Fs statistic  and the 95% confidence interval around these statistics using the bootstrap method (with no recombination) offered in DNAsp  with 5000 replicates.
We amplified the DNA at 14 microsatellite loci in T. alpinus and T. minimus (Loci names: EuAmMS26, EuAmMS37, EuAmMS41, EuAmMS86, EuAmMS94, AC A2, AC A101, AC A108, AC B12, AC B111, AC C2, AC C122, AC D107 and AC D115). The first five were taken from the literature , and the remaining nine were developed in T. alpinus. Primer sequences of all 14 loci are available in Additional file 3: Table S2. Reverse primers were fluorescently labeled with one of the following dyes: PET, NED, FAM, or HEX, forward primers were unlabeled. PCR reactions with a volume of 8.0 μl contained reagents in the following concentrations: 0.5-1 μl DNA template, 0.25 μM each primer, 0.2 mM each dNTP, 0.8 μl 10X BSA, 0.8 μl 10X PCR buffer (Roche), 1.5 mM MgCl2 and 0.4U of Taq DNA polymerase (Roche). The thermal cycler consisted of 94°C for 2 min, followed by 30 cycles of 94°C for 40s, 51-60°C for 40s, and 72°C for 40s, and ending with a final extension at 72°C for 10mins. Locus-specific annealing temperatures are shown Additional file 1: Table S1). PCR products were sized by capillary electrophoresis on an ABI 3730 sequencer (Applied Biosystems, Inc.), and alleles were scored manually using program GENEMAPPER Ver. 4.0 software (Applied Biosystems, Inc.). Positive and negative controls as well as three replicate samples were run on each PCR plate for each locus. Repeat genotypes showed high repeatability.
We tested for linkage disequilibrium and deviations from Hardy-Weinberg equilibrium (HWE) in each locus, across populations and overall with an exact test using (10000 permutations; . Significant heterozygote deficiencies were used to identify the presence of null alleles as well as using the program FreeNA  to detect the frequencies of null alleles in our dataset. Bonferroni corrections for multiple tests were applied to p-values . To examine population structure, we applied the Bayesian approach implemented in the software Structure 2.3.3  to identify clusters of randomly mating individuals with minimum HW deviations and linkage disequilibrium. We ran the admixture model with correlated allele frequencies with five replicates of 106 Markov Chain Monte Carlo (MCMC) iterations after a burnin of 105 from K (number of parental populations) = 1 to K = 10. To provide the most accurate estimation of K, we used the statistic ∆K introduced by Evanno et al. . We averaged coefficients of membership across the five replicates using the software CLUMMP 1.1  and DISTRUCT 1.1  was used to plot the graphical representation of this membership. To further examine genetic structure we used the program Arelquin  to calculate pair-wise FST values. To visualize the genetic distance, we generated a neighbor-joining tree using the pairwise FST distances in the program MEGA version 4 .
To infer the divergence history between T. alpinus and T. minimus, we used the coalescent-based isolation-with-migration (IM) model  implemented in the program IMa2 . We estimated the following parameters: effective population size of T. alpinus (NeALP), T. minimus (NeMIN) and their common ancestor (NeA), the migration rate from T. alpinus into T.minimus (mALP->MIN) and from T. minimus to T. alpinus (mMIN->ALP), and finally time since divergence (t). IMa2 first uses a Bayesian Markov Chain Monte Carlo (MCMC) approach to integrate over the space of possible genealogies and divergence times then uses the genealogies to estimate the posterior distribution of effective population sizes and migration rates to calculate joint posterior probability of all model parameters [65–67]. We used 10 loci (cyt b sequences, and 9 microsatellite loci: EuAmMS26, EuAmMS41 EuAmMS86, EuAmMS94, EuAmMS37, ACA101, ACA108, ACC2, ACD115) partitioned by species in this analysis. Five microsatellite loci used in the population genetic analyses have complex repeat motifs and therefore may not follow a strict step-wise mutation model. We sub-sampled the entire dataset to improve computational efficiency. Thirty individuals from the microsatellite dataset and twenty individuals from the sequence dataset were randomly chosen from each species for the analysis, with assurance that each geographic area was represented. We used a two-population model where each species was considered a “population”. A series of preliminary runs were used to estimate upper bounds on priors and assess mixing. Our final run consisted of 60 chains (geometric heating scheme set ha = 0.980, hb = 0.50), a burnin of 3 × 105 steps followed by 30 × 106 steps sampling trees from each locus every 300 steps (ESS > 50). We did two replicates of the final run starting with a different random number seed. Each run took approximately 72 days to finish and both returned parameter estimates that were near identical. Two hundred thousand saved genealogies (100,000 from each run) were used to calculate the joint posterior probability of the parameters in L-mode of IMa2. We used a general mammalian nucleotide substitution rate weighted across sites of 8.2 × 10-9; see also ) to calculate the locus-wide mutation rate for the 801 bp segment of cyt b to be approximately 6.6 × 10–6, and the average mutation rate for the microsatellites as 1.0 × 10-4. These mutation rates were used to convert the parameter estimates into demographic units (i.e., time in years, population size in individuals and migration rates as individuals/generation). Finally, nested models were tested to determine if the full model fit the data significantly better than models when population sizes and/or migration rates were set to equal or zero.
Mitochondrial sequence data
Sample size, number of haplotypes detected, haplotype diversity (h d ) nucleotide diversity (θ π , θ S ) and tests of population expansion/contraction (Tajima’s D, Fu’ Fs statistics) in T. minimus and T. alpinus , and geographic groups of each species (Figure 1) at the mitochondrial gene, cyt b
No. of haplotypes
T. min- Wht/Iny
The statistical parsimony haplotype network for T. alpinus and T. minimus had a 95% parsimony limit of 12 steps (Figure 3a). The mtDNA phylogenetic tree estimated by the Bayesian analysis was weakly resolved but supports the network analysis and demonstrates a lack of clear genealogical separation between the two species (Additional file 5: Figure S4). There are four groups of haplotypes that are separated by at least 5 base pair changes and show some geographic structure (Haplogroups North1, North2, South1 and South2, Figure 3b) and three out of four of these contain individuals of both species (North1, North2 and South2). South1 is made up of only T. alpinus individuals from the southern portion of their range. A map of the haplogroups and shared haplotypes shows that genetic similarity is more defined by geography than by species identity (Figure 3b). Two southern T. alpinus haplotypes were more than 12 steps away from the others in the network and are not shown in Figure 3.
Pairwise comparisons of Tamias populations
Consistent with the above, the AMOVA attributed 28.63% of the genetic variation across haplotypes to be between species (ϕCT = 0.28, p = 0.34), 37.6% to be among populations within species (ϕSC = 0.52, p < 0.001) and 33.8% of the variation to differences within species (ϕST = 0.66, p < 0.0001). By comparison, the AMOVA for the analysis of nuclear loci (see next section) attributed 62% to variation among groups (ϕCT = 0.62, p = 0.14), 5.1% among populations within groups (ϕSC = 0.13, p < 0.001) and 32.9% to variation within populations (ϕST = 0.67, p < 0.001). The AMOVA analyses reveal that mtDNA variation is not explained by differences among species or geographic groups, but rather differences within species and populations.
Sample size (N), average allelic richness ( A ; corrected for differences in sample size), observed and expected heterozygosity (H o , H e & standard deviation (sd)) in T. alpinus and T. minimus at 14 microsatellite loci
The pairwise FST values between clusters showed significant differentiation across all clusters and ranged from 0.018 between T.min-N and T.min-C to 0.227 between T.alp-N and T.min-S (Table 2). In contrast to the results for mtDNA, the neighbor-joining tree based on FST shows that groups of the same species are genetically more similar to each other than to the other species (Additional file 8: Figure S7) although the southern populations of both species (T. alp-S and T.min-S) are both differentiated from their more northern conspecifics.
Coalescent analysis of divergence history
Parameter estimates from IMa2 runs for T. alpinus and T. minimus
2NmT.min to T. alp
2NmT.alp to T.min
We examined the evolutionary relationship of T. alpinus and T. minimus using cytochrome b and microsatellites to help elucidate the divergence history of T. alpinus in the Sierra Nevada. Microsatellite analyses were used to provide a contemporary view of this relationship and to examine details of population genetic structure within and across species. We found that T. alpinus and T. minimus populations share mitochondrial haplotypes with no overall geneaological separation, and that diversity at this locus is better explained by geography than by species’ boundaries. This pattern indicates either recent speciation of T. alpinus from T. minimus with retention of ancestral polymorphism, or extensive introgression subsequent to splitting. In contrast to the mtDNA sequence data, the analyses of nuclear microsatellite loci and morphology revealed that the two species are genetically distinct. Although there are highly differentiated populations within species, populations of the same species are more similar to each other than they are to members of the other species. This suggests that contemporary hybridization is not widespread along the geographic boundary between T. minimus and T. alpinus. Coalescent analysis of divergence history revealed mid to late Pleistocene divergence and low, but significant gene flow between the two lineages. Overall, our study suggests that: speciation between T. alpinus and T. minimus is relatively recent, secondary contact and historical introgression has occurred, and there is little evidence or opportunity for contemporary hybridization given the current distribution of these two species.
Differentiating between the genetic signature of incomplete lineage sorting and historical hybridization is difficult; however, distinguishing between the two is important in addressing non-concordance among characters in closely related species . The spatial pattern of genetic variation across species can help to provide an objective assessment of which process is more likely to have occurred because each should produce a specific spatial pattern [30, 72]. The results of our mitochondrial sequences show no strong spatial structure within species. Both the Bayesian tree and haplotype network indicate the two species are completely intermingled across the landscape. Recent hybridization should show a clustered pattern, where introgressed alleles are more common at or near the contact zone of the two species; in contrast, ancestral polymorphism should be diffuse and uniform across space . However, the spatial patterns described above assume a contact zone between species, and for this study, we were unable to obtain samples of both species from the same location. The one site where they were collected in close proximity in 1911 (Little Brush Meadow, Tulare Co., CA) showed distinct and divergent haplotypes between species. The single T. minimus taken from this locality had a unique haplotype that clustered with the southern minimus haplogroup in the network (South2 Figure 3, Additional file 5: Figure S4) whereas the T. alpinus from that location and three other T. alpinus from the general area collected at the same time, all clustered within southern T. alpinus haplogroup (South1). Overall, the distribution of shared haplotypes reflects geographic proximity, not current contact; northern T. alpinus and T. minmus shared three haplotypes and southern T. alpinus and central and southern T. minimus shared one. This geographic pattern, along with the results of the IMa2 analysis suggest that hybridization occurred at some point in the past but given the results of the microsatellite analysis and current distribution of the two species, hybridization is not currently ongoing.
There were potential problems with the parameter estimates using IMa2, in particular, that for divergence time. The first is that our microsatellite data exhibit genetic structure within lineages (recall that in this analysis, each species was considered a “population”), which may lead to an overestimation of divergence time . Second, the right tail of the posterior density plot reached low but non-zero values potentially making the 95% HPD high and low estimates unreliable (i.e. posterior density distribution did not converge within prior range). Given these caveats, the estimate of divergence time between T. alpinus and T. minimus from IMa2 is about 450 ka with some support of it being as recent as 110 ka (the lower 95% HPD of t). This timeframe occurs in the mid to late Pleistocene, a time of extreme climate fluctuations in the Sierra Nevada including several major glaciations with prolonged glacial rather than interglacial periods . The coalescent analysis also supports low but significant gene flow subsequent to splitting from T. minimus into T. alpinus. Given the relatively small population size estimated by the IMa2 analysis for T. alpinus compared to T. minimus and considering that T. alpinus is a range restricted endemic, it is conceivable that the T. alpinus population would receive more genes from T. minimus than the other way around.
The Pleistocene shaped the genetic structure and distributions of many species (reviewed in ) and its role in the speciation in several North American taxa, though sometimes debated, is clear (e.g., [74–76]). It is plausible that a founding population of T. alpinus became isolated in a refugium from an adjacent T. minimus population. A recent study that used sequence data from reproductive protein genes showed that T. alpinus, characterized by strongly divergent sequences, is a monophyletic group nested within T. minimus species complex . This nested pattern is consistent with expectations for a peripheral isolate of a widespread species providing further support to our results suggesting that speciation of T. alpinus is recent. The IMa2 analysis also supports this notion with a parameter estimate for the effective population size of T. alpinus (NeALP = 430,625), substantially smaller than the estimate of the effective population size of T. minimus (NeMIN = 1,448,317). Finally, although the estimates of migration between species were low (95% HPD high <1 from T. minimus to T. alpinus), there is significant evidence of historical genetic introgression between species, adding to the complexity of the divergence dynamics between these two species. The results of this study adds to the growing number of studies that have highlighted the importance of glacial-interglacial refugia in recently derived species across taxa in the central and southern Sierra Nevada (e.g., [12, 14]).
Hybridization can play an important role in the evolution of species [77–79]. It was previously accepted that the morphological differences in bacular morphology in western chipmunks mechanically prevented hybridization between species and was considered a strong pre-mating barrier to gene flow [43, 80]. However, several recent studies have documented both historical and ongoing hybridization in two non-sister species of Tamias and suggest that it may be more common in the genus than previously thought [30, 31, 42]. Our analyses provide little or no evidence for current introgression across the species boundaries; however, the potential for contemporary gene flow between these species where they come into close proximity in the southern portion of their range exists. One apparent hybrid individual, morphologically determined to be a T. alpinus with a divergent T. alpinus haplotype, was more similar to T. minimus than T. alpinus based on 14 microsatellite loci. The assignment of this individual to the T.min-N population is puzzling because the two localities are far apart geographically (at least 80 km). Despite extensive survey effort along the eastern flank of the Sierra Nevada, we have not found a locality of co-occurrence where there is direct potential for hybrid matings. Furthermore, our morphological analyses show no evidence of hybridization and species were easily discriminated based on the measured morphological characters.
Speciation is a prolonged process that likely has phases in different spatial contexts [81, 82]. Here, we provided evidence of recent speciation with limited yet significant post-divergence gene flow between the Alpine chipmunk and its closest relative. Our approach revealed an interesting and complex pattern of shared and intermingled haplotypes across species and highly differentiated populations within. The presence of geographic structure of shared haplotypes between these closely related lineages could be a result of allopatric speciation with secondary contact, but the same pattern could be a result of zones of primary parapatric speciation. As theory suggests, any mechanism that can cause divergence in allopatry can also occur in parapatry, as long as the selective gradient acting on differentiation is strong enough to counterbalance continuous gene flow . In order to further examine the genetic patterns that emerged in this study, future work comparing genome-level patterns of diversity between these two species may reveal regions of the genome under strong divergent selection; a pattern not expected in an allopatric speciation scenario . Additionally, environmental niche models (ENMs) used to predict historical ranges aid in the identification of areas of isolation and the potential for secondary contact in the past. The environmental modeling results can be combined with genetic data to test alternative hypotheses of allopatric speciation with secondary contact or primary parapatric speciation (e.g. ). A multifaceted approach including building historical ENMs and genomic data combined with increased sampling could further improve our understanding of the evolutionary history of Alpine chipmunk, a species that appears to be under threat due to recent climate change [70, 86]. A clear understanding of this species evolutionary history could help us understand its vulnerability in the face of environmental change.
We are grateful to Chris Conroy, Jessica Castillo, Karen Rowe, Kevin Rowe, Carol Patton, Les Chow and the YNP resurvey team for their role in collecting the field data and Marisa Lim for her assistance in the lab. Sonal Singhal provided valuable advice for implementing the IMa2 analyses. The manuscript was improved through feedback from Sonal Singhal, Sean Rovito, Juan Parra, Margarita Hadjistylli, Pauline Kamath, Kari Roesch-Goodman, Justin Brashares and three anonymous reviewers. E. M. Rubidge was supported by a National Science & Engineering Research Council (NSERC) PGS-D award, the Museum of Vertebrate Zoology, and the Environmental Science, Policy and Management Department at UC Berkeley, during this research. The project was funded by the Museum of Vertebrate Zoology at UC Berkeley, the Yosemite Fund, the National Geographic Society and the National Science Foundation.
- Dobzhansky T: Genetics and the Origin of Species. 1951, New York: Columbia University Press, 3Google Scholar
- Schluter D: Ecology and the origin of species. Trends Ecol Evol. 2001, 16: 372-380. 10.1016/S0169-5347(01)02198-X.PubMedView ArticleGoogle Scholar
- Mallet J: Hybrid speciation. Nature. 2007, 446: 279-283. 10.1038/nature05706.PubMedView ArticleGoogle Scholar
- Pinho C, Hey J: Divergence with gene flow: models and data. Annu Rev Ecol Evol Syst. 2010, 41: 215-230. 10.1146/annurev-ecolsys-102209-144644.View ArticleGoogle Scholar
- Abbott R, Albach D, Ansell S, Arntzen JW, Baird SJE, Bierne N, Boughman JW, Brelsford A, Buerkle CA, Buggs R, Butlin RK, Dieckmann U, Eroukhmanoff F, Grill A, Cahan SH, Hermansen JS, Hewitt G, Hudson AG, Jiggins C, Jones J, Keller B, Marczewski T, Mallet J, Martinez-Rodriguez P, Moest M, Mullen S, Nichols R, Nolte AW, Parisod C, Pfennig K, et al: Hybridization and speciation. J Evol Biol. 2013, 26: 229-246. 10.1111/j.1420-9101.2012.02599.x.PubMedView ArticleGoogle Scholar
- Avise JC, Walker D, Johns GC: Speciation durations and Pleistocene effects on vertebrate phylogeography. Proc R Soc B Biol Sci. 1998, 265: 1707-1712. 10.1098/rspb.1998.0492.View ArticleGoogle Scholar
- Hewitt GM: Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc B. 2004, 359: 183-195. 10.1098/rstb.2003.1388.View ArticleGoogle Scholar
- Avise JC: Phylogeography: The History and Formation of Species. 2000, Cambridge, Massachusetts: Harvard University PressGoogle Scholar
- Gillespie AR, Zehfuss PH: Glaciations of the Sierra Nevada, California, USA. In Quaternary Glaciations and Chonology. Volume 2b, Part II: North America, Developments in Quartenary Science. Edited by Ehlers J, Gibbard PL. Amsterdam: Elsevier; 2004.Google Scholar
- Grayson DK: The Great Basin: a Natural History. 2011, Berkeley, CA: University of California PressGoogle Scholar
- Feldman CR, Spicer GS: Comparative phylogeography of woodland reptiles in California: repeated patterns of cladogenesis and population expansion. Mol Ecol. 2006, 15: 2201-2222. 10.1111/j.1365-294X.2006.02930.x.PubMedView ArticleGoogle Scholar
- Schoville SD, Roderick GK: Alpine biogeography of Parnassian butterflies during Quaternary climate cycles in North America. Mol Ecol. 2009, 18: 3471-3485. 10.1111/j.1365-294X.2009.04287.x.PubMedView ArticleGoogle Scholar
- Hull JM, Keane JJ, Savage WK, Godwin SA, Shafer JA, Jepsen EP, Gerhardt R, Stermer C, Ernest HB: Range-wide genetic differentiation among North American great gray owls (Strix nebulosa) reveals a distinct lineage restricted to the Sierra Nevada, California. Mol Phylo Evol. 2010, 56: 212-221. 10.1016/j.ympev.2010.02.027.View ArticleGoogle Scholar
- Rovito SM: Lineage divergence and speciation in the Web-toed Salamanders (Plethodontidae: Hydromantes) of the Sierra Nevada, California. Mol Ecol. 2010, 19: 4554-4571. 10.1111/j.1365-294X.2010.04825.x.PubMedView ArticleGoogle Scholar
- Poage MA, Chamberlain CP: Stable isotopic evidence for a Pre-Middle Miocene rain shadow in the western Basin and Range: implications for the paleotopography of the Sierra Nevada. Tectonics. 2002, 21 (4):Google Scholar
- Mulch A, Sarna-Wojcicki AM, Perkins ME, Chamberlain CP: A Miocene to Pleistocene climate and elevation record of the Sierra Nevada (California). Proc Natl Acad Sci USA. 2008, 105: 6819-6824. 10.1073/pnas.0708811105.PubMedPubMed CentralView ArticleGoogle Scholar
- Grayson DK: The Late Quaternary biogeographic histories of some Great Basin mammals (western USA). Quat Sci Rev. 2006, 25: 2964-2991. 10.1016/j.quascirev.2006.03.004.View ArticleGoogle Scholar
- Brown JH: Mammals on mountaintops - nonequilibrium insular biogeography. Am Nat. 1971, 105: 467-10.1086/282738.View ArticleGoogle Scholar
- Grayson DK: Mammalian responses to Middle Holocene climatic change in the Great Basin of the western United States. J Biogeography. 2000, 27: 181-192. 10.1046/j.1365-2699.2000.00383.x.View ArticleGoogle Scholar
- Rickart EA: Elevational diversity gradients, biogeography and the structure of montane mammal communities in the intermountain region of North America. Global Ecol Biogeography. 2001, 10: 77-100. 10.1046/j.1466-822x.2001.00223.x.View ArticleGoogle Scholar
- Guyton B: Glaciers of California: Modern Glaciers, Ice Age Glaciers, Origin of Yosemite Valley, and a Glacial Tour in Sierra Nevada. 1998, Berkeley and Los Angeles: University of California PressGoogle Scholar
- Matthes F: Report of the committee on glaciers 1939-1940. Trans Am Geophys Union. 1940, 518: 396-406.View ArticleGoogle Scholar
- Glazner AF, Manley CR, Marron JS, Rojstaczer S: Fire or ice: Anticorrelation of volcanism and glaciation in California over the past 800,000 years. Geophys Res Lett. 1999, 26: 1759-1762. 10.1029/1999GL900333.View ArticleGoogle Scholar
- Basagic HJ, Fountain AG: Quantifying 20th century glacier change in the Sierra Nevada, California. Arct Antarct Alp Res. 2011, 43: 317-330. 10.1657/1938-4246-43.3.317.View ArticleGoogle Scholar
- Gompert Z, Fordyce JA, Forister ML, Shapiro AM, Nice CC: Homoploid hybrid speciation in an extreme habitat. Science. 2006, 314: 1923-1925. 10.1126/science.1135875.PubMedView ArticleGoogle Scholar
- Schoville SD, Stuckey M, Roderick GK: Pleistocene origin and population history of a neoendemic alpine butterfly. Mol Ecol. 2011, 20: 1233-1247. 10.1111/j.1365-294X.2011.05003.x.PubMedView ArticleGoogle Scholar
- Good JM, Sullivan J: Phylogeography of the red-tailed chipmunk (Tamias ruficaudus), a northern Rocky Mountain endemic. Mol Ecol. 2001, 10: 2683-2695. 10.1046/j.0962-1083.2001.01397.x.PubMedView ArticleGoogle Scholar
- Piaggio A, Spicer GS: Molecular phylogeny of the chipmunks inferred from mitochondrial cytochrome b and cytochrome oxidase II gene sequences. Mol Phylo Evol. 2001, 20: 335-350. 10.1006/mpev.2001.0975.View ArticleGoogle Scholar
- Demboski JR, Sullivan J: Extensive mtDNA variation within the yellow-pine chipmunk, Tamias amoenus (Rodentia: Sciuridae), and phylogeographic inferences for northwest North America. Mol Phylo Evol. 2003, 26: 389-408. 10.1016/S1055-7903(02)00363-9.View ArticleGoogle Scholar
- Good JM, Demboski JR, Nagorsen DW, Sullivan J: Phylogeography and introgressive hybridization: Chipmunks (genus Tamias) in the northern Rocky Mountains. Evolution. 2003, 57: 1900-1916.PubMedView ArticleGoogle Scholar
- Good JM, Hird S, Reid N, Demboski JR, Steppan SJ, Martin-Nims TR, Sullivan J: Ancient hybridization and mitochondrial capture between two species of chipmunks. Mol Ecol. 2008, 17: 1313-1327. 10.1111/j.1365-294X.2007.03640.x.PubMedView ArticleGoogle Scholar
- Reid N, Demboski JR, Sullivan J: Phylogeny estimation of the radiation of Western North American Chipmunks (tamias) in the face of introgression using reproductive protein genes. Syst Biol. 2012, 61: 44-62. 10.1093/sysbio/syr094.PubMedPubMed CentralView ArticleGoogle Scholar
- Levenson H, Hoffmann RS, Nadler CF, Deutsch L, Freeman SD: Systematics of the Holarctic chipmunks (Tamias). J Mamm. 1985, 66: 219-242. 10.2307/1381236.View ArticleGoogle Scholar
- Hoffmann RS: Different voles for different holes: environmental restrictions on refugial survival of mammals. Proceedings of the Second International Congress of Systematic and Evolutionary Biology. 1981, Vancouver, BC, 25-45.Google Scholar
- Hofmann RS, Anderson CG, Thorington RJ, Heaney L: Family Sciurdae. Mammal Species of the World: A Taxonomic and Geographic Reference. Edited by: Wilson DE, Reeder DM. 1993, Washington: Smithsonian Institute Press, 419-465.Google Scholar
- Bergstrom BJ, Hoffmann RS: Distribution and Diagnosis of 3 species of chipmunks (Tamias) on the front range of Colorado. Southwest Nat. 1991, 36: 14-28. 10.2307/3672112.View ArticleGoogle Scholar
- Verts BJ, Carraway LN: Tamias minimus. Mamm Species. 2001, 653: 1-10.View ArticleGoogle Scholar
- Sullivan RM: Phyletic, biogeographic, and ecologic relationships among montane populations of least Chipmunks (eutamias minimus) in the Southwest. Syst Zool. 1985, 34: 419-448. 10.2307/2413206.View ArticleGoogle Scholar
- Johnson DH: Systematic review of the chipmunks (genus Eutamias) of California. Univ Calif Publ Zool. 1943, 48: 63-143.Google Scholar
- Sutton DA: Tamias amoenus. Mamm Species. 1992, 390: 1-8.View ArticleGoogle Scholar
- Grinnell J, Storer TI: Animal Life in Yosemite. 1924, Berkeley, CA: University of California PressGoogle Scholar
- HIRD S, Sullivan J: Assessment of gene flow across a hybrid zone in red-tailed chipmunks (Tamias ruficaudus). Mol Ecol. 2009, 18: 3097-3109. 10.1111/j.1365-294X.2009.04196.x.PubMedView ArticleGoogle Scholar
- Patterson BD, Thaeler CSJ: The Mammalian Baculum: hypotheses on the nature of bacular variability. J Mamm. 1982, 63: 1-15. 10.2307/1380665.View ArticleGoogle Scholar
- Mullen LM, Hoekstra HE: Natural selection along an environmental gradient: classic cline in mouse pigmentation. Evolution. 2008, 62: 1555-1570. 10.1111/j.1558-5646.2008.00425.x.PubMedView ArticleGoogle Scholar
- Smith MF: Phylogenetic relationships and geographic structure in Pocket Gophers in the Genus Thomomys. Mol Phylo Evol. 1998, 9: 1-14. 10.1006/mpev.1997.0459.View ArticleGoogle Scholar
- Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.PubMedView ArticleGoogle Scholar
- Posada D: jModelTest: phylogenetic model averaging. Mol Biol Evol. 2008, 25: 1253-1256. 10.1093/molbev/msn083.PubMedView ArticleGoogle Scholar
- Templeton AR, Crandall KA, Sing CF: A cladistic-analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data 3. Cladogram estimation. Genetics. 1992, 132: 619-633.PubMedPubMed CentralGoogle Scholar
- Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000, 9: 1657-1659. 10.1046/j.1365-294x.2000.01020.x.PubMedView ArticleGoogle Scholar
- Nei M: Molecular Evolutionary Genetics. 1987, New York, NY: Columbia University PressGoogle Scholar
- Tamura K, Dudley J, Nei M, Kumar S: MEGA4: molecular evolutionary genetics analysis (MEGA) software version 4.0. Mol Biol Evol. 2007, 24: 1596-1599. 10.1093/molbev/msm092.PubMedView ArticleGoogle Scholar
- Librado P, Rozas J: DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009, 25: 1451-1452. 10.1093/bioinformatics/btp187.PubMedView ArticleGoogle Scholar
- Excoffier L, Laval G, Schneider S: Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinforma. 2005, 1: 47-50.Google Scholar
- Tajima F: Statistical-method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123: 585-595.PubMedPubMed CentralGoogle Scholar
- Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997, 147: 915-925.PubMedPubMed CentralGoogle Scholar
- Schulte-Hostedde AI, Gibbs HL, Millar JS: Microsatellite DNA loci suitable for parentage analysis in the yellow-pine chipmunk (Tamias amoenus). Mol Ecol. 2000, 9: 2180-2181. 10.1046/j.1365-294X.2000.105314.x.PubMedView ArticleGoogle Scholar
- Raymond M, Rousset F: An exact test for population differentiation. Evolution. 1995, 49: 1280-1283. 10.2307/2410454.View ArticleGoogle Scholar
- Chapuis MP, Estoup A: Microsatellite null alleles and estimation of population differentiation. Mol Biol Evol. 2006, 24: 621-631. 10.1093/molbev/msl191.PubMedView ArticleGoogle Scholar
- Rice WR: Analyzing tables of statistical tests. Evolution. 1989, 43: 223-225. 10.2307/2409177.View ArticleGoogle Scholar
- Pritchard JK, Stephens M, Donnelly P: Inference of population structure using multilocus genotype data. Genetics. 2000, 155: 945-959.PubMedPubMed CentralGoogle Scholar
- Evanno G, Regnaut S, Goudet J: Detecting the number of clusters of individuals using the software structure: a simulation study. Mol Ecol. 2005, 14: 2611-2620. 10.1111/j.1365-294X.2005.02553.x.PubMedView ArticleGoogle Scholar
- Jakobsson M, Rosenberg NA: CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007, 23: 1801-1806. 10.1093/bioinformatics/btm233.PubMedView ArticleGoogle Scholar
- Rosenberg NA: DISTRUCT: a program for the graphical display of population structure. Mol Ecol Resour. 2004, 4: 137-138.Google Scholar
- Nielsen R, Wakeley J: Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics. 2001, 158: 885-896.PubMedPubMed CentralGoogle Scholar
- Hey J: Isolation with migration models for more than Two populations. Mol Biol Evol. 2010, 27: 905-920. 10.1093/molbev/msp296.PubMedPubMed CentralView ArticleGoogle Scholar
- Hey J, Nielsen R: Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D-persimilis. Genetics. 2004, 167: 747-760. 10.1534/genetics.103.024182.PubMedPubMed CentralView ArticleGoogle Scholar
- Hey J, Nielsen R: Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proc Natl Acad Sci USA. 2007, 104: 2785-2790. 10.1073/pnas.0611164104.PubMedPubMed CentralView ArticleGoogle Scholar
- Pesole G, Gissi C, de Chirico A, Saccone C: Nucleotide substitution rate of mammalian mitochondrial genomes. J Mol Evol. 1999, 48: 427-434. 10.1007/PL00006487.PubMedView ArticleGoogle Scholar
- Schlötterer C: Evolutionary dynamics of microsatellite DNA. Chromosoma. 2000, 109: 365-371. 10.1007/s004120000089.PubMedView ArticleGoogle Scholar
- Rubidge EM, Patton JL, Lim M, Burton AC, Brashares JS, Moritz C: Climate-induced range contraction drives genetic erosion in an alpine mammal. Nat Clim Chang. 2012, 2: 285-288. 10.1038/nclimate1415.View ArticleGoogle Scholar
- Avise JC, Ball MR: Principles of genealogical concordance in species concepts and biological taxonomy. Oxford Surveys in Evolutionary Biology. Edited by: Futuyma D, Antonovics J. 1990, New York: Oxford University Press, 45-67.Google Scholar
- Goodman SJ, Barton NH, Swanson G, Abernethy K, Pemberton JM: Introgression through rare hybridization: a genetic study of a hybrid zone between red and sika deer (genus Cervus) in Argyll, Scotland. Genetics. 1999, 152: 355-371.PubMedPubMed CentralGoogle Scholar
- Wakeley J: The effects of subdivision on the genetic divergence of populations and species. Evolution. 2000, 54: 1092-1101.PubMedView ArticleGoogle Scholar
- Knowles LL, Otte D: Phylogenetic analysis of montane grasshoppers from western North America (Genus Melanoplus, Acrididae: melanoplinae). Ann Entomol Soc Am. 2000, 93: 421-431. 10.1603/0013-8746(2000)093[0421:PAOMGF]2.0.CO;2.View ArticleGoogle Scholar
- Knowles LL: Did the Pleistocene glaciations promote divergence? Tests of explicit refugial models in montane grasshopprers. Mol Ecol. 2001, 10: 691-701.PubMedView ArticleGoogle Scholar
- Johnson NK, Cicero C: New mitochondrial DNA data affirm the importance of Pleistocene speciation in North American birds. Evolution. 2004, 58: 1122-1130.PubMedView ArticleGoogle Scholar
- Arnold ML: Natural hybridization as an evolutionary process. Ann Rev Ecol Syst. 1992, 23: 237-261. 10.1146/annurev.es.23.110192.001321.View ArticleGoogle Scholar
- Seehausen O: Hybridization and adaptive radiation. Trends Ecol Evol. 2004, 19: 198-207. 10.1016/j.tree.2004.01.003.PubMedView ArticleGoogle Scholar
- Rieseberg LH: Hybrid origins of plant species. Ann Rev Ecol Syst. 1997, 28: 359-389. 10.1146/annurev.ecolsys.28.1.359.View ArticleGoogle Scholar
- Adams DR, Sutton DA: A description of bacullum and os clitoris of Eutamias townsendii ochrogenys. J Mamm. 1968, 49: 764-10.2307/1378743.View ArticleGoogle Scholar
- Butlin RK, Galindo J, Grahame JW: Sympatric, parapatric or allopatric: the most important way to classify speciation?. Philos Trans R Soc B. 2008, 363: 2997-3007. 10.1098/rstb.2008.0076.View ArticleGoogle Scholar
- Butlin R, Debelle A, Kerth C, Snook RR, Beukeboom LW, Cajas RFC, Diao W, Maan ME, Paolucci S, Weissing FJ, van de Zande L, Hoikkala A, Geuverink E, Jennings J, Kankare M, Knott KE, Tyukmaeva VI, Zoumadakis C, Ritchie MG, Barker D, Immonen E, Kirkpatrick M, Noor M, Macias Garcia C, Schmitt T, Schilthuizen M, Network MCS: What do we need to know about speciation?. Trends Ecol Evol. 2012, 27: 27-39.PubMedView ArticleGoogle Scholar
- Turelli M, Barton NH, Coyne JA: Theory and speciation. Trends Ecol Evol. 2001, 16: 330-343. 10.1016/S0169-5347(01)02177-2.PubMedView ArticleGoogle Scholar
- Slatkin M: Gene flow in natural populations. Ann Rev Ecol Syst. 1985, 16: 393-430. 10.1146/annurev.ecolsys.16.1.393.View ArticleGoogle Scholar
- Rossetto M, Allen CB, Thurlby KAG, Weston PH, Milner ML: Genetic structure and bio-climatic modeling support allopatric over parapatric speciation along a latitudinal gradient. BMC Evol Biol. 2012, 12: 1-16. 10.1186/1471-2148-12-1.View ArticleGoogle Scholar
- Rubidge EM, Monahan WB, Parra JL, Cameron SE, Brashares JS: The role of climate, habitat, and species co-occurrence as drivers of change in small mammal distributions over the past century. Global Change Biol. 2011, 17: 696-708. 10.1111/j.1365-2486.2010.02297.x.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.