Climate change during the late Quaternary has been implicated as the cause of both massive range shifts and extinction events. We combined molecular marker data and previously published fossil data to reconstruct the late Quaternary history of a grassland-dependent species, the black-footed ferret (Mustela nigripes), and to determine whether populations from Pleistocene refugia in the Columbia Basin, eastern Beringia, and Great Plains persisted into the Holocene and Recent eras. Using DNA extracted from 97 museum specimens of extirpated populations, we amplified 309 bp of the mtDNA control region, and 8 microsatellite markers from the nuclear genome. Overall haplotype diversity from 309 base pairs (bp) of mitochondrial DNA (mtDNA) control region was low (5 haplotypes, nucleotide diversity = 0.001 ± 0.001 SD) and was contained within a single phylogenetic clade. The star phylogeny and unimodal mismatch distribution indicated that a rapid range expansion from a single Pleistocene refugium occurred. Microsatellite data corroborated this genetic pattern: populations from the mixed grasslands of the Great Plains had significantly higher expected heterozygosity and allelic richness than populations to the west (HE = 0.66 versus 0.41, AR = 4.3 versus 2.7, respectively), and 𝜰, a measure of relative population size, was substantially greater in the east than west (2.4 versus 0.7). We infer from these data that black-footed ferrets rapidly colonized western ecoregions in a stepwise fashion from the Great Plains to the intermountain regions of the Rocky Mountains and the Colorado Plateau after the last ice age. It appears that glacial retreat and global warming caused both range expansion and localized extinction in this North American mustelid species.
The late Quaternary was a period marked by massive environmental and biotic change. At the end of the Pleistocene (14,500 years ago), mean annual temperatures increased by 5–7°C (Webb et al. 1993), ice sheets retracted from Holarctic regions of the Northern Hemisphere, and more than 70% of large mammal species became extinct (Martin 1984). Populations that had been confined to islands of suitable habitat expanded as ice sheets withdrew and ecosystems evolved (Graham et al. 1996). For many species, this period was important to the origination of genetic diversity and population structuring in the form of distinct genetic lineages generated from isolation in separate refugia (Hewitt 2000; Steele and Storfer 2006). These range shifts and the resulting genetic consequences molded the evolutionary history of species affected by global climate change (Hewitt 2004).
Fossil, palynological, and genetic data have all been used to document the range shift of plant and animal species from Pleistocene refugia to Holocene and eventually modern distributions. In particular, use of genetic data has documented the pattern of expansion of numerous boreal and montane mammalian species from low-latitude or low-elevation refugia following the warming trend, including long-tailed vole (Microtus longicaudus—Conroy and Cook 2000), dusky shrew (Sorex monticolus), northern flying squirrel (Glaucomys sabrinus), black bear (Ursus americanus—Lessa et al. 2003), European pine marten (Martes martes), and polecat (Mustela putorius—Davison et al. 2001). Ultimately, identifying patterns of refugia and expansion for multiple endemic or habitat-dependent species can lead to hypothesis testing on the formation and evolution of ecosystems (Carstens et al. 2005). To date, the vast majority of North American phylogenetic studies have focused on forest- or mesic-adapted species, which have greatly increased our understanding of how these ecosystems responded to late Quaternary events (Hewitt 2004).
Grassland ecosystems also responded to the movement patterns of glacial ice (Betancourt et al. 1990; Dort and Jones 1970); however, there is a dearth of corroborating genetic data. During the late Pleistocene glacial maxima, the distribution of grassland ecosystems in North America was drastically reduced (Bryson et al. 1970), and areas occupied by these ecosystems today held vastly different communities of plants and animals during the late Pleistocene (Pielou 1991). For example, the Great Plains was dominated by open-canopy park woodlands with areas of boreal forest in the central Plains (Corner and Voorhies 1994), and the lower elevations of the Colorado Plateau were composed of pine woodland and juniper–sagebrush desert scrub (Betancourt 1990). As the climate warmed during the Holocene, grasses and open-canopy habitat increased in abundance throughout the Great Plains (Graham and Mead 1987), Colorado Plateau (Betancourt 1990), and intermountain regions of the Rocky Mountains (Barnosky et al. 1987), bringing with them a variety of grassland-dependent fauna (Graham et al. 1996). We investigated the genetic consequences of a late Quaternary geographic range shift and expansion to a grassland-dependent species, the black-footed ferret (Mustela nigripes).
The black-footed ferret is a small-bodied mustelid whose modern distribution is tightly associated with prairie dogs (Cynomys); indeed, black-footed ferrets live in the burrow systems of prairie dogs, which are their primary prey (Clark 1994). Systematic eradication of prairie dogs throughout the 20th century led to the near extinction of this highly endangered mustelid. Phylogenetic evidence suggests that black-footed ferrets evolved from their sister taxon, Mustela eversmanii, and crossed the Beringian land bridge sometime between 2 million and 500,000 years ago (O'Brien et al. 1989). The earliest fossil evidence of black-footed ferrets in North America is mid-Pleistocene (∼800,000 years ago—Anderson 2004; Owen et al. 2000). Wisconsinian fossil records (30,000–14,500 years ago) indicate the persistence of black-footed ferrets during the last glacial maximum, and they are described from southern Columbia Basin, west of the Continental Divide, as well as from the Great Plains east of the Continental Divide (Anderson et al. 1986). For species with modern distributions east and west of the Rocky Mountains, it appears that populations on either side of the Continental Divide were isolated from one another to the extent that they formed distinct genetic lineages. Examples include red-tailed chipmunk (Tamias ruficaudus—Good and Sullivan 2001), the flightless montane grasshopper Melanoplus oregonensis (Knowles 2001), the swift and kit foxes (Vulpes velox and V. macrotis—Mercure et al. 1993), and the North American desert spider Agelenopsis aperta (Ayoub and Riechert 2004). We studied the available genetic evidence from modern black-footed ferret populations (known only from east of the Continental Divide and the Colorado Plateau) to identify extant genetic lineages formed during the Pleistocene. Using mitochondrial DNA (mtDNA), we inferred the number and location (east or west of the Continental Divide) of late Pleistocene refugia. Should the descendents of more than 1 refugium have persisted into modern times, we would expect distinct mitochondrial genetic lineages in modern samples. Using more rapidly evolving microsatellite markers from the nuclear genome, we further describe patterns of Holocene expansion in this grassland-dependent carnivore. As temperatures rose and grasslands expanded during the early Holocene, we hypothesized that the expansion of black-footed ferrets succeeded the expansion of grasslands and their graminivorous prey.
Materials and Methods
Sampling, extraction, and amplification
Black-footed ferrets are a highly endangered North American carnivore (IUCN 1988). Once extirpated from the wild, black-footed ferrets have since been reintroduced throughout their historical distribution. Because all reintroduced black-footed ferrets originated from 1 population, we used DNA from museum specimens (n = 86) or frozen skeletal muscle (n = 11) of modern individuals collected before their extirpation for our analyses. We successfully extracted DNA from 97 samples of black-footed ferrets from 8 North American museums (Appendix I). The average date of collection was 1926 ± 3.0 years (mean ± SE, range = 1883–1985) and sampling locations were representative of the modern range of the black-footed ferret (Fig. 1A). Black-footed ferrets were only known west of the Continental Divide in the Columbia Basin and Great Basin as recently as the Holocene.
We extracted DNA from 2 sources: frozen tissues and museum specimens. DNA was extracted from frozen skeletal muscle using a standard phenol–chloroform extraction procedure as previously described (Wisely et al. 2002). Extraction and amplification of DNA from these samples occurred in a molecular biology facility. We extracted DNA from epithelial and bone tissue of museum specimens following the protocol of Wisely et al. (2004b). Because DNA from these sources is of low yield, samples can easily be contaminated with modern DNA and polymerase chain reaction product; therefore, we used rigorous protocols and multiple negative controls (Miller et al. 2002). Based on criteria of Gilbert et al. (2005), we assessed the risk of erroneous data from contamination to be low. DNA was not from a common source (i.e., human or domestic animal), and samples were of relatively recent origin (<200 years old). Nonetheless, contamination from polymerase chain reaction product is always a concern. To lessen the risk of contamination, samples from museum specimens were extracted in a laboratory designed for the extraction and polymerase chain reaction preparation of DNA from low-quality sources. To minimize contamination, this laboratory was housed in a building where no other work on DNA was conducted. All personnel entering this facility showered before entry and wore protective apparel dedicated to the laboratory. All polymerase chain reaction amplification and downstream applications were conducted in another building. Negative controls were used in both extraction procedures (1 in 8) and polymerase chain reaction preparation (1 in 16) to ensure that contamination did not occur.
We amplified and analyzed a 309-base pair (bp) segment of the mtDNA control region. Because spurious results can be obtained when trying to amplify large fragments from low-quality DNA (Gilbert et al. 2005), we used 2 primer pairs that amplified shorter fragments; bffmt8 (5′-ATTCCCTGATTTTCTCACCA-3′) and bffmt9 (5′-AGAAGAGGTGCACGTTGAGG-3′), which yielded a product of 244–248 bp, and bffmt14 (5′-TCCACCTCACTTAGATCACGAG-3′) and bffmt15 (5′-ACCAAATGCATGACACCACA-3′), which yielded a product of 258 bp. Primers were designed based on published GenBank sequences for the black-footed ferret (accession nos. AF068571 [Davison et al. 1999] and AJ548475 [Michaux et al. 2004]). Polymerase chain reaction products were sequenced in both directions using an ABI 3100 automated sequencer (Applied Biosystems, Foster City, California) following the manufacturer's protocols. We excluded a 73- to 77-bp segment within the bffmt8–bffmt9 fragment that was difficult to resolve because of the presence of a mononucleotide, hypervariable repeat.
We amplified 8 microsatellite loci: Mvis075, Mer095, Mvis002 (Fleming et al. 1999), Mer049 (Wisely et al. 2002), Mvi57, Mvi87, Mvi232 (O'Connell et al. 1994), and G1A (Paetkau and Strobeck 1994) using polymerase chain reaction conditions previously published for DNA from museum specimens of this species (Wisely et al. 2004b). These conditions minimized the presence of stutter bands, which can be scored as false alleles (Taberlet et al. 1996). Still, genomic DNA extracted from museum specimens typically yields DNA of low quality and quantity and allelic dropout rates can be high (Taberlet et al. 1996; Wisely et al. 2004b). We quantified allelic dropout rate and determined whether allele specific bias in dropout rate occurred. To quantify allelic dropout, each genotype was amplified more than once (on average 2.4 times). If an allele was observed in one amplification, but not others, we assumed that allele to be real, because the probability of allelic dropout (no amplification of a true allele) is much greater than the probability of a false allele (artifact due to polymerase chain reaction conditions). Evidence suggests that false alleles occur only rarely (3–5% of total amplifications in low-quality DNA samples—Miller and Waits 2003; Taberlet et al. 1996), whereas allelic dropout is more common (17–50%—Miller et al. 2002; Miller and Waits 2003; Taberlet et al. 1996). Products were visualized using a ABI 373 automated sequencer (Applied Biosystems).
Mitochondrial DNA data analysis
We estimated gene diversity and nucleotide diversity for the 309-bp segment of the control region. The relationship among haplotypes was mapped using a minimum spanning network generated by ARLEQUIN 2.0 (Schneider et al. 2000). The proportion of total genetic variation attributable to among-population versus within-population variation (φST) was not estimated because several populations were monotypic.
To examine population dynamics, we tested the hypothesis of a recent population expansion using a mismatch analysis (Rogers and Harpending 1992). We estimated the number of nucleotide differences between every pair of individuals. Under an equilibrium expectation, the curve generated by the relative frequency of nucleotide differences should be multimodal. The distribution of values should be unimodal when a recent population expansion has occurred. We compared the observed distribution of values to the expected equilibrium distribution using 2 indices of fit: sum of squares differences (SSD— Schneider and Excoffier 1999) and an index of raggedness (r—Harpending 1994); we estimated the probability that the distribution was unimodal by bootstrapping with 10,000 resamples. We further examined population dynamics using Tajima's test of neutrality (Tajima 1989); assuming that sequence variation in the control region is neutral, a significant negative value of Tajima's D indicates a population expansion, whereas a value significantly greater than zero represents a population bottleneck.
We calculated the percent variation contained within and among populations using an analysis of molecular variance (AMOVA—Excoffier et al. 1992). We used population definitions based on ecoregion; we then grouped data into 3 northern (ecoregions Nmix, Cmix, and WYshrub, as defined in Fig. 1A) and 3 southern populations (ecoregions Smix, Short, and COshrub). All analyses were conducted using ARLEQUIN 2.0.
Microsatellite data analysis
Examination of data indicated that allelic dropout did occur, but without bias to allele identity (see “Results”). When dropout events are random, it should not affect estimates of allele frequency within or among populations (Miller et al. 2002; Miller and Waits 2003). Because the focus of our research was not on individual genetic identity, but rather on population-level processes, a modest level of allelic dropout was tolerated. To reduce analytical imprecision, bias, or inaccuracy, we minimized analyses that assumed genetic equilibrium or required individual genotypic data. Instead, we focused on analyses that used allele frequency or gene genealogy data, which are more robust to violations of assumptions of Hardy–Weinberg equilibrium (Pearse and Crandall 2004).
To compare regional diversity, differentiation, and gene flow among northern and southern ecoregions, and among eastern (Nmix, Cmix, and Smix) and western (Short, COshrub, and WYshrub) ecoregions (Fig. 1A), we a priori assigned individuals to ecoregions using the Bayesian assignment method of Rannala and Mountain (1997) with the program GENECLASS 1.0.02 (Cornuet et al. 1999). Grouping samples into ecoregions (hereafter called populations) allowed us to calculate population genetic parameters among regions with different late Quaternary histories. GENECLASS calculates the probability densities of allele frequencies for each population, and then uses them to generate the probability of a particular genotype in a particular population. GENECLASS assigns an individual to the population that has the highest log likelihood of assignment. This assignment test did not assume genetic equilibrium in the population; rather, probabilities are based on allele frequencies.
We also used a Bayesian clustering method that assumes no a priori knowledge of sample location or number of populations using the program STRUCTURE (version 2.0—Pritchard et al. 2000). The program estimates the probability of a set of allelic frequencies given K populations (Pr(X/K)), and allows for an admixture model, that is, a proportion, qk, of an individual's genome comes from population k, such that Σ qk = 1 per individual. Individuals were assigned to the population with the highest qk value. We estimated log Pr(X/K) for K = 1–7 using a Markov chain Monte Carlo simulation approach. For each value of K, we ran 107 simulations; the first 106 simulations were used as burn-in for the model. We ran each set of simulations 3 times per K to ensure the stability of the parameter estimation model. We chose K based on the value of the posterior probability, that is, the value at which K reached an asymptote.
We tested Hardy–Weinberg equilibrium for populations created by both assignment procedures, GENECLASS and STRUCTURE, using the Markov chain method and for linkage equilibrium using the program GENEPOP (Raymond and Rousset 1995). We calculated gene diversity and allelic richness for ecoregions with the program FSTAT (Goudet 2001). We tested for differences in the mean and variance of gene diversity and allelic richness between individuals pooled from western populations versus eastern populations with an independent sample t-test and Levene's test for equality of variances (Zar 1984) using the statistical program SPSS (version 1.01; SPSS Inc., Chicago, Illinois).
We compared 2 distance measures with different underlying assumptions: chord distance (Dc—Cavalli-Sforza and Edwards 1967) is a geometric estimate of distance with no biological assumptions, and FST (Weir and Cockerham 1984) is a measure of population differentiation that assumes that loci are in Hardy–Weinberg equilibrium. Because FST remains the standard for interspecific comparisons of genetic structure (Pearse and Crandall 2004), we calculated FST values. We tested how allelic dropout rate could differentially affect genetic distance values by testing for similarity of pairwise patterns of differentiation between FST and chord distance using a Mantel test.
We used a coalescent model of microsatellite evolution to calculate 𝜰, an estimate of relative population size using the program MIGRATE (Beerli and Felsenstein 2001). The microsatellite model of evolution was approximated with Brownian motion. In an initial run of the program, we used estimates of FST to calculate starting values of 𝜰. Resulting values from the initial run were used as starting values for 3 subsequent runs. We used default search parameters for every run; all runs (except for the initial one) gave similar results. For maximum likelihood estimation, we selected a full migration model to estimate 𝜰.
Results
Mitochondrial DNA control region
We obtained forward and reverse consensus sequence data on 69 of 97 sequences and found 5 haplotypes (GenBank AM746622–AM7466226; Table 1). All but 1 sequence varied by 1 base transition from the most common haplotype to produce a star phylogeny (Fig. 1B). A star phylogeny describes a shallow haplotype tree where most haplotypes coalesce to a common, abundant haplotype (Avise 2000; Slatkin and Hudson 1991). The abundant haplotype comprised 91% of the individuals sampled and was found in all ecoregions that we sampled. Each of the 4 other haplotypes was unique to an ecoregion (Fig. 1A); thus, AMOVA results indicated that 100% of the molecular variance resided within populations (Table 2). All 5 haplotypes were found in eastern ecoregions (Nmix, Cmix, and Smix); however, only the most common haplotype was found in western ecoregions (Short, COshrub, and WYshrub; Fig. 1A).
Gene diversity (Nei 1987) was 0.19 ± 0.06; nucleotide diversity was also low (0.001 ± 0.001 SD). Two lines of evidence suggested that a sudden range expansion occurred. We found a significant departure from the expected distribution of values generated by Tajima's test of selective neutrality indicative of a range expansion (D = −1.547, P = 0.05). In addition, the observed mismatch distribution fit the predicted unimodal distribution of a sudden range expansion (Fig. 2) using both indices of fit (sum of squares differences, SSD < 0.001, P = 0.48; and index of raggedness, r = 0.44, P = 0.70, where P was the probability that the mismatch distribution was unimodal). The starlike relationship of haplotypes to one another (Fig. 1B) further suggests that there was 1 ancestral haplotype (a) from which the other haplotypes were recently derived to form 1 haplogroup. Haplotypes b, c, and d differed from haplotype a by 1 nucleotide substitution, whereas e differed by 2 substitutions.
Microsatellite markers
Across the 8 loci, we successfully amplified unambiguous DNA product on average 2.4 times per genotype out of an average of 3.8 attempts per genotype. We found evidence of allelic dropout in our data set. Of the 220 genotypes that we scored as heterozygous at least once in 572 replicate amplifications of those genotypes, 24% of the 572 amplifications (135 amplifications) amplified only 1 allele. In 50% of the 34 cases where dropout was observed more than once in a replicate amplification, both alleles dropped out in separate replicate amplifications. Dropout, therefore, appeared to be random and not biased with respect to allele designation.
Results from the program GENECLASS indicated that 60% of individuals correctly assigned to the ecoregion from which they were sampled (Table 3). The northern and central mixed grasslands (Nmix and Cmix) had high rates of reciprocal misclassification (50% and 13%, respectively) and the 3 southern ecoregions (Smix, Short, and COshrub) had high rates of misassignment within those 3 regions (28 of 51 individuals assigned to an ecoregion that was different from its a priori assignment), suggesting that the ecoregion subdivision had weak genetic support. Individuals from the Wyoming Basin shrub ecoregion were correctly assigned to their ecoregions, which suggests that this population was isolated from other ecoregions.
Using the program STRUCTURE, we found that log Pr(X/K) exhibited very low variation among simulation sets, thus it appeared that our number of simulations was sufficient to accurately estimate the parameters. Values of −log (X/K) asymptoted at K = 6; therefore, we assumed 6 populations in further analyses using STRUCTURE. Forty-four of 97 individuals exhibited admixture of population ancestry; we considered an individual's genome to be admixed if the largest qk of an individual was less than 0.80. The program STRUCTURE confirmed population substructure between individuals from the Wyoming Basin shrub ecoregion and the rest of the ecoregions; all individuals from the Wyoming Basin were assigned to population 2 (Table 4). The other 5 populations, k = 1 and k = 3–6, showed no clear geographic patterning; population designations of individuals, as assigned by STRUCTURE, were distributed throughout ecoregions.
We found no evidence of linkage disequilibrium. When individuals were assigned to ecoregions, 4 of 6 populations had an excess of homozygotes (Table 5); however, when no a priori assignment of individuals to populations was made, only 1 population (population 3) was out of Hardy–Weinberg equilibrium at 3 loci that had a heterozygote deficit (Mvis002, Mer095, and Mvi57; data not shown).
Eastern ecoregions (Nmix, Cmix, and Smix) had the highest number of private alleles, whereas the western ecoregions (Short, COshrub, and WYshrub) had 16 times fewer unique alleles than ecoregions to the east (Table 5). Western ecoregions also had significantly lower measures of diversity than eastern ecoregions; in the 3 western ecoregions, expected heterozygosity was 62% of the 3 eastern mixed grassland ecosystems (Student's t = 4.6, d.f. = 36, P < 0.001, d.f. adjusted for unequal variance; Levene's test F = 7.4, P = 0.009) and western ecoregions had significantly less allelic richness (corrected for variation in sample size) than eastern ecoregions (Student's t = 5.9, d.f. = 46, P < 0.001, d.f. not adjusted; Levene's test was not significant; Table 5).
The FST and chord distance were significantly correlated (P < 0.001; Table 6). The effective number of migrants per generation, as calculated and averaged from pairwise FST values was 1.4 migrants per generation. The highest gene flow occurred among the 3 southern ecoregions, Smix, Short, and COshrub. Estimates of gene flow confirmed assignment tests; WYshrub was the most isolated population, having the largest estimates of genetic differentiation and the lowest estimate of 𝜰. 𝜰 was substantially lower in the western ecoregions (Short, COshrub, and WYshrub); the average 𝜰 of these regions was 0.7, which was 29% of the average value for the 3 eastern ecoregions (Table 6).
Discussion
Pleistocene refugium
Because the noncoding mtDNA control region evolves 5–10 times faster than the coding regions of the mitochondrial genome (Gemmell et al. 1996), it is often used to infer late Pleistocene and Holocene events. We found few haplotypes and little difference among haplotypes in 309 bases of the mtDNA control region. The short branches of the minimum-spanning network are indicative of a Holocene-aged expansion from a single population in which haplotype a predominated, and further suggest that modern black-footed ferrets are descended from a single Pleistocene refugium. The lack of diversity and shortness of the branches stand in contrast to the 800,000-year fossil record of black-footed ferrets in North America. Should more than 1 lineage have persisted, we would have expected to find deep divergence among haplogroups. Given that the museum collections from which we sampled were representative of the modern distribution of black-footed ferrets, this result is most parsimoniously explained by the persistence of 1 haplogroup into modern times.
Low nucleotide diversity and the absence of ancient lineages also have been reported for other species of holarctic mustelids both in North America and Europe, including European polecat–domestic ferret (M. putorius–M. furo—Davison et al. 2001), wolverine (Gulo gulo—Walker et al. 2001), European otter (Lutra lutra—Cassens et al. 2000), and fisher (Martes pennanti—Drew et al. 2003). For instance, nucleotide diversity that we found (0.001 ± 0.001) in black-footed ferrets was similar to that found in European polecats throughout Europe (0.009 ± 0.001) and similarly limited in the number of base-pair substitutions (Davison et al. 1999, 2001). Davison et al. (2001) argued that the genetic pattern of low haplotype diversity in present-day European polecat and pine marten (Martes martes) to be evidence of colonization from a single Pleistocene refugium. Our pattern of low diversity and a star-patterned phylogeny also is consistent with a single late Pleistocene refugium, descendants of which persisted into modern times (Avise 2000). Because all haplotype diversity originated east of the Continental Divide, we hypothesize that this refugium resided in the Great Plains of North America.
Despite molecular evidence for a single refugium, fossil evidence of black-footed ferrets whose remains date between 14,500 and 30,000 years old (the span of the Wisconsinian glacial maximum) indicated that multiple refugia persisted during the glacial maximum. In addition to fossil remains from the Great Plains, fossils were found in the southern Columbia Basin (ColShrub; Fig. 1), and eastern Beringia, which indicated that this species persisted well into the Holocene in these regions (Anderson et al. 1986). These refugia would have been separated from the Great Plains by the Cordilleran ice sheet and a series of mountain top glaciers throughout the Rocky Mountains south of the ice sheet (Hollin and Schilling 1981; Pielou 1991). Had multiple refugial lineages persisted into modern populations of black-footed ferrets, more deeply divergent mitochondrial haplotypes that formed into east–west haplogroups should have been detected. We found no evidence of divergent haplotypes indicative of multiple refugia or secondary contact in modern populations, but a nearly monotypic haplotype that extended throughout the modern distribution. We postulate that the greater Columbia Basin refugial population of black-footed ferrets (ColShrub and GBshrub) was extirpated before modern times based on examination of our data and the emerging cross-species phylogeographic pattern of the region. For black-footed ferrets, it appears that modern intermountain populations and populations in the Colorado Plateau to the south and west of the Rocky Mountains (areas that were under ice or boreal forest during the last glacial maximum) were colonized from a refugial population east of the Continental Divide. Multiple records of black-footed ferrets from the late Pleistocene from Nebraska, eastern Wyoming, and Alberta, Canada, confirm the presence of black-footed ferrets in this eastern refugium (Anderson et al. 1986).
Holocene expansion
Evidence of postglacial expansion from a single Great Plains refugium is tripartite. First, the star phylogeny pattern produced by control region haplotypes is indicative of a recent range expansion from a single source (Avise 2000). Second, the unimodal pattern of the frequency of pairwise mismatch distribution was consistent with rapid range expansion (Fig. 2). A significant negative value for Tajima's D provided further evidence that the mismatch distribution had been affected by a rapid range expansion. Finally, 100% of the genetic variation in mitochondrial haplotypes was found within, not among, populations, suggesting little differentiation among populations, which is consistent with a rapid range expansion.
The pattern of genetic diversity and gene flow found in the microsatellite data further supports the idea of a Holocene range expansion from a single source and provides evidence for the location of this refugium. Comparisons of genetic diversity between eastern and western ecoregions (average allelic richness: 4.3 versus 2.7; average expected heterozygosity: 0.66 versus 0.41; average number of private alleles: 11 versus <1) indicate higher genetic diversity in the eastern mixed prairie grassland ecosystems (Nmix, Cmix, and Smix), as would be expected for the source of a population expansion. A clinal decrease in microsatellite diversity has previously been used to infer rapid postglacial expansion in another mustelid, the fisher, from western coastal mountain ranges (Wisely et al. 2004a). The east-to-west pattern of decreasing genetic diversity found in black-footed ferrets also is indicative of successive founder events as populations expanded west from the mixed grassland prairie of the Great Plains. The high relative population size, 𝜰, of eastern populations compared to western populations (average 𝜰 = 2.4 versus 0.7; Table 6) further supports the hypothesis of multiple founder events as populations spread westward.
Despite the rapid expansion of black-footed ferrets into the emerging grassland and shrub ecosystems, population substructure was evident among black-footed ferret samples. Pairwise FST values ranged from 0.03 between COshrub and Smix to 0.69 between WYshrub and Short, and the program STRUCTURE identified 6 clusters into which individuals were assigned. Nonetheless, these clusters had no clear geographic pattern and 45% of the 97 samples had admixed genomes from different population clusters. This pattern suggests that population differentiation occurred across the ecoregions, yet limited but ongoing gene flow contributed to the admixture of the differentiating populations. This pattern of divergence with gene flow is expected for patchily distributed populations, which are heavily influenced by dispersal events to produce populations with incomplete isolation (Mayr 1963; Wright 1931). The patchy distribution of prairie dog complexes (Koford 1958) produced a similar pattern of genetic structure within the species (Antolin et al. 2006) and may explain the same pattern of genetic structure in black-footed ferrets, the obligate predators of prairie dogs.
The Wyoming Basin shrub population, found at the periphery of the modern distribution of black-footed ferrets and the source of all extant ferrets, merits extra consideration. Both FST and chord distance values were highest between the Wyoming Basin shrub and all other ecoregions and 𝜰 was smallest for this population (Table 6). Both assignment tests indicated that this population was clearly defined and differentiated from other populations, and this population also had the lowest genetic diversity (Table 5). The increased isolation and genetic drift and the low genetic diversity we found in this population are consistent with the pattern of genetic structure found for populations at the periphery of a species' distribution (Lesica and Allendorf 1995). Examination of the microsatellite data suggests that once colonized, this population became isolated during the Holocene. Yet, despite differentiation suggested by microsatellite data, all WYshrub animals had the most common mtDNA haplotype, suggesting that this population was part of the eastern refugium and that isolation occurred after the last glacial maximum.
Our findings have conservation implications. All extant black-footed ferrets are descended from animals from the Wyoming Basin shrub. The high level of isolation, low level of diversity, and small population size of the source population suggest that genetic factors may diminish the long-term evolutionary potential and viability of this species (Frankel and Soulé 1980). Efforts to conserve the remaining genetic diversity, as outlined in the species recovery plan, should be continued in earnest, including maximizing the retention of diversity in the captive population, genetic augmentation of reintroduced populations with captive individuals, and a plan for incorporating undiscovered remnant populations should the situation arise.
Conclusions
Analysis of DNA from museum specimens is a powerful tool for understanding phylogeographic and population-level processes of populations for which there are no modern representatives. Extensive holdings of black-footed ferret specimens at museums across North America allowed for detailed analysis of past population dynamics and processes. Despite evidence of allelic dropout and nonequilibrium conditions, we were able to elucidate phylogeographic patterns using analyses that relied on allele frequency data rather than equilibrium conditions.
Analysis of the location and origin of forest refugia has dominated phylogeographic literature of North America. By contrast, relatively little is known about the time-transgressive process that replaced forest and parkland with grassland in the Great Plains (Pielou 1991). Our study is one of a handful to discuss late Quaternary dynamics of a grassland species. Examination of the data suggests that, as grassland habitat expanded westward during the ensuing warming trend, a grassland-dependent species, the black-footed ferret, expanded westward as well. Furthermore, the eastern Great Plains were an important Pleistocene refuge for grassland species that now occupy the western Great Plains, Colorado Plateau, and intermountain basins of the Rocky Mountains.
The volatile climatic changes that occurred at the end of the Pleistocene heralded massive changes in ecosystem composition and function. These changes contributed to the extinction of large numbers of mammalian taxa and a redistribution of the taxa that remain today (Guthrie 1984). Two ecoregions in which black-footed ferrets were distributed during the late Quaternary experienced substantial decreases in graminoid species during the Holocene. Eastern Beringia experienced a decrease in graminoid plant community diversity (Guthrie 1984) and large portions of the Columbia Basin's grassland community were replaced by forest ecosystems during the Holocene (Blinnikov et al. 2002; Whitlock et al. 2000). No black-footed ferrets exist there now, and our analysis suggests that these populations left no modern descendents. These northern and western refugia for black-footed ferrets in the Columbia Basin and eastern Beringia appear to have followed the trend of other grassland mammal populations from the late Quaternary and were extirpated.
Acknowledgments
We thank J. Maldonado, S. Lance, and J. Ortega for statistical and laboratory advice and C. Dorsey for laboratory assistance. Thanks are also extended to C. McIntosh for assistance with the program MIGRATE, and to 2 anonymous reviewers for their valuable input. We thank the curators and collection managers at the institutions that allowed us to sample museum material, especially R. Humphrey at University of Colorado Museum, C. Jones at Denver Museum of Science and Nature, C. Conroy at University of California Museum of Vertebrate Zoology, S. McLaren at Carnegie Museum of Natural History, J. Yearous at North Dakota State Historical Society, C. Ramotnik at Museum of Southwestern Biology, M. Gosselin at Canadian Museum of Nature, and B. Fisher at National Museum of Natural History. We thank P. Else for providing line art. Funding for the work was provided by the Smithsonian Institution Fellowship Program, the Smithsonian Institution's Abbott Fund, Friends of the National Zoo, and Kansas State University Division of Biology.
Literature Cited
Appendices
Appendix I
Accession numbers of museum specimens from which we acquired DNA for this study. Carnegie Museum of Natural History: CMNH19392, CMNH20628; Denver Museum of Science and Nature: DMNH257, DMNH1208, DMNH1558, DMNH1559, DMNH1684, DMNH1883, DMNH1987, DMNH2024, DMNH2247, DMNH2248, DMNH3206, DMNH3644, DMNH3703, DMNH5199, DMNH5792, DMNH6033, DMNH8274; North Dakota State Historical Society: NDSHS5159, NDSHS11689; Canadian Museum of Nature: NMC11700, NMC11703, NMC11752, NMC12682, NMC14078, NMC14079, NMC14095, NMC24235; National Museum of Natural History: NMNH15470, NMNH34977, NMNH35011, NMNH65061, NMNH83992, NMNH83994, NMNH168741, NMNH180719, NMNH188450, NMNH188451, NMNH188452, NMNH188453, NMNH188454, NMNH188455, NMNH188456, NMNH188457, NMNH188458, NMNH201945, NMNH209150, NMNH211513, NMNH228233, NMNH228789, NMNH232400, NMNH234138, NMNH234970, NMNH234971, NMNH234972, NMNH234973, NMNH241014, NMNH243818, NMNH243819, NMNH243820, NMNH243990, NMNH245641, NMNH247073, NMNH248973, NMNH251453, NMNH285877, NMNH289498, NMNH013996 A21066, NMNH015471 A22427, NMNH019263 A35376, NMNH019294 A35017, NMNH019295 A35018, NMNH022537 A30064, NMNH025358 A32771, NMNH19263-35016; Museum of Vertebrate Zoology: UCB55212, UCB95039; University of Colorado Museum: UCM59, UCM287, UCM458, UCM493, UCM553, UCM963; Museum of Southwestern Biology: MSB 39933, MSB 39932, MSB 39515, MSB 39386, MSB 39385, MSB 39384, MSB 39383, MSB 39382, MSB 39381, MSB 39378, MSB 39379, MSB 39822, MSB 39373, MSB 39377.
Fig. 1.—A) Distribution of grassland and shrubland ecoregions in North America that contained black-footed ferrets (Mustela nigripes) either in Quaternary or modern times. Nmix: northern mixed grassland, Cmix: central mixed grassland, Smix: southern mixed grassland, Short: western short grassland, COshrub: Colorado plateau shrubland, WYshrub: Wyoming basin shrubland, ColShrub: southern Columbia Basin shrubland, GBshrub: Great Basin shrubland. Dark shaded areas are ecoregions that are associated with forests. Ecoregions are based on Ricketts et al. (1999). Nmix, Cmix, Smix, Short, COscrub, and WYscrub ecoregions represent the modern distribution of black-footed ferrets. No black-footed ferrets are known from ColShrub or GBshrub after the Holocene. Circle centroids represent locations at which black-footed ferrets were collected; circle size represents sample size: small (n = 1 or 2), medium (n = 3–5), large (n = 6–10), largest (n = 12). The bold black line indicates the Continental Divide. B) Minimum spanning network based on 309 base pairs of the mitochondrial control region. Numbers in brackets are the sample sizes for each haplotype, which is represented by a circle. Circles with shading represent haplotype a, vertical striped areas represent haplotype b, horizontal striped areas represent haplotype c, stippled areas represent haplotype d, and white areas represent the relative frequency of haplotype e. Circle size is indicative of the relative frequency of each haplotype.
Fig. 2.—Distribution of observed (filled circle, dashed line) and expected (open diamond, solid line) pairwise differences in substitutions among haplotype sequences of Mustela nigripes. Expected values were generated in ARLEQUIN 3.01 using a sudden stepwise expansion model using the estimated population sizes before (θ0 = 0.19) and after (θ1 = 0.24) the expansion, estimated to have occurred in τ = 3.0 generations. Parameters are expressed in units of mutational time.
Table 1.—Variable sites of 5 D-loop haplotypes (a–e) for black-footed ferret (Mustela nigripes). M. eversmanii (steppe polecat, ME in table) is included for reference (Genbank AJ548477—Michaux et al. 2004). The 309-base pair sequences of M. eversmanii compared here correspond to base numbers 55–363.
Table 2.—Components of an analysis of molecular variance (AMOVA—Excoffier et al. 1992) of DNA control region sequences of mitochondrial DNA from 6 populations of Mustela nigripes (n = 69). We used population definitions based on ecoregion; we then grouped data into 3 northern (ecoregions Nmix, Cmix, and WYshrub, as defined in Fig. 1A) and 3 southern populations (ecoregions Smix, Short, and COshrub). All analyses were conducted using ARLEQUIN 2.0.
Table 3.—Proportion of individual black-footed ferrets (Mustela nigripes) that were sampled in an ecoregion (see Fig. 1 for ecoregion definitions) and the subsequent ecoregion assignment using the Bayesian assignment procedure in GENECLASS.
Table 4.—Average proportion of membership (Qk, the population average of individual qk; calculated in STRUCTURE) of individual genotypes of black-footed ferrets (Mustela nigripes) from each ecoregion given K = 6 populations.
Table 5.—Sample size (n), expected heterozygosity (HE), observed heterozygosity (HO), number of unique alleles per population (U), average number of alleles per locus corrected for sample size (allelic richness [AR]), and the percentage of individuals correctly classified into populations by the assignment test (%C) for 6 historical populations of the black-footed ferret (Mustela nigripes).
Table 6.—Pairwise population calculations of chord distance (Dc, below diagonal) of Cavalli-Sforza and Edwards (1967) and an estimate of population differentiation (FST, above diagonal—Weir and Cockerham 1984) for 6 populations of Mustela nigripes. Estimates of relative population size, 𝜰, calculated from MIGRATE (Beerli and Felsenstein 2001) are given in bold on the diagonal. All parameters were estimated from the microsatellite marker data set.
