http://orcid.org/0000-0001-9580-6987
Figures
Abstract
Genetic variation is the basis upon which natural selection acts to yield evolutionary change. In a rapidly changing environment, increasing genetic variation should increase evolutionary potential, particularly for small, isolated populations. However, the introduction of new alleles, either through natural or human-mediated processes, may have unpredictable consequences such as outbreeding depression. In this study, we identified a contact zone and limited gene flow between historically separated genetic lineages of American pikas (Ochotona princeps), representing the northern and southern Rocky Mountain subspecies, within Rocky Mountain National Park. The limited spatial extent of gene flow observed may be the result of geographic barriers to dispersal, selection against hybrid individuals, or both. Our fine-scale population genetic analysis suggests gene flow is limited but not completely obstructed by extreme topography such as glacial valleys, as well as streams including the Colorado River. The discovery of two subspecies within this single protected area has implications for monitoring and management, particularly in the light of recent analyses suggesting that the pikas in this park are vulnerable to fragmentation and local extinction under future projected climates. Future research should focus on the fitness consequences of introgression among distinct genetic lineages in this location and elsewhere, as well as within the context of genetic rescue as a conservation and management strategy for a climate sensitive species.
Citation: Castillo Vardaro JA, Epps CW, Frable BW, Ray C (2018) Identification of a contact zone and hybridization for two subspecies of the American pika (Ochotona princeps) within a single protected area. PLoS ONE 13(7): e0199032. https://doi.org/10.1371/journal.pone.0199032
Editor: Tzen-Yuh Chiang, National Cheng Kung University, TAIWAN
Received: December 20, 2017; Accepted: May 30, 2018; Published: July 11, 2018
Copyright: © 2018 Castillo Vardaro et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All raw sequence data is available publicly from NCBI (https://www.ncbi.nlm.nih.gov/genbank/), GenBank accession numbers MG700594-MG700613. All microsatellite genotypes are within S3 File.
Funding: This project was funded by the National Parks Service (https://www.NPS.gov) Climate Change Response Program (PMIS #163377), awarded to CWE and CR. Additional support was provided by the National Science Foundation Graduate Research Fellowship Program (to JACV). Publication of this paper was supported, in part, by the Henry Mastin Graduate Student Fund, awarded to CWE. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Intraspecific genetic diversity is the most fundamental element of biodiversity and provides the basis for natural selection to yield evolutionary change [1–3]. When confronted with rapid environmental change, populations must either adapt in situ to new conditions, shift their distribution to more favorable environmental conditions, or face extinction [4]. Therefore, understanding patterns of genetic diversity and population structure is essential for developing and implementing effective conservation and management strategies for at risk populations [5, 6]. Local, regional, and historical processes all contribute to contemporary patterns of genetic diversity and population differentiation. Populations at the extremes of a species’ range, either geographic or climatic, may be most vulnerable to rapid, contemporary climate change. However, they may also represent reservoirs of adaptive potential if there is local adaptation to extreme environmental conditions [5–7].
Anthropogenic manipulation of genetic structure, whether intentional or accidental, also may influence adaptive potential. Translocations and augmentations, whereby individuals from one locality are moved to another location to found a new population or to supplement an existing population [8, 9] have been used to combat declining populations for decades. There are many well-known examples of mammalian reintroductions stemming from translocations, including wolves [10] and bighorn sheep [8, 11], as well as augmentations including, panthers [12] and bighorn sheep [8]. In addition to traditional attempts at demographic rescue by increasing population numbers, some more recent interventions have focused on genetic rescue which aims to increase population resilience by increasing overall genetic diversity and or targeting specific traits such as resistance to disease [12]. Best practices dictate that translocated individuals should come from closely-related genetic stock to avoid admixture among distinct evolutionary units that may have negative consequences such as maladaptation and outbreeding depression [11, 13]. Yet, some have raised the question of whether, in the face of extinction, more dramatic interventions such as interbreeding distantly related populations with different biogeographic history could prove beneficial [9]. For example, interbreeding with distantly related genetic lineages could introduce novel alleles or reintroduce ancestral alleles that were lost through drift, possibly having positive consequences for fitness [14]. Evolutionary rescue [15] and the potential genetic consequences of translocations [11, 16] are active and extensive areas of research. Nevertheless, these are difficult processes to study empirically outside a laboratory setting.
Naturally-occurring hybrid zones among distinct and formerly isolated genetic lineages provide opportunities to learn more about introgression, such as how genes spread through admixed populations and the subsequent effects on physiology and behavior [17, 18]. Among the well-studied natural hybrid zones, many are the result of secondary contact that occurred following the last ice age and thus occur in geographic clusters sometimes referred to as “suture zones” [17, 19–23]. Well-studied zones include mountain ranges in Europe and western North America [17, 19], and such taxa as grasshoppers in the French Alps [21] and gophers in the Rocky Mountains [24]. These studies shaped our early understanding of genetic introgression in hybrid zones. Barton and Hewitt [25, 26] argued that many hybrid zones are clines maintained by a balance between dispersal and selection against hybrids, which they refer to as “tension zones”. The characteristics of the cline (e.g., steepness and width) are therefore determined by the rate of gene flow, dispersal distance, and relative fitness of alleles between the two genetic sources [25, 26]. Recent developments in sequencing technologies and the inclusion of population genomic data into studies of hybrid zones has greatly increased our understanding of the evolutionary consequences of hybridization, but the increase in studies and associated data has also illuminated the complexity and variability of possible outcomes of hybridization [18].
In this study, we describe a previously undocumented hybrid zone between two distinct genetic lineages of American pikas, Ochotona princeps, within Rocky Mountain National Park (ROMO) and discuss the implications for conservation of the species within this management unit and elsewhere. American pikas are often considered a sentinel species with respect to climate change due to numerous observations of local extinctions, particularly in lower, drier, and hotter habitats [27–30]. Recent extirpations at higher elevations among habitat thought to be more favorable to pikas have further worried biologists and managers [31]. Pikas are small lagomorphs (121–176 g, Fig 1) found throughout much of the intermountain western United States [32]. They are restricted to fractured rock habitats, such as talus slopes and lava flows, which provide refuge from predators and thermal buffering [32–36]. They cannot tolerate prolonged exposure to high temperatures and are therefore typically found at high elevations, but may persist at lower elevations and in hotter climates if there are suitable microclimatic refugia [37–40]. Predictive modeling has suggested widespread losses in the species’ distribution particularly in, but not limited to, low elevations [41]. More recently, models accounting for shifts in functional connectivity as well as climatic variables suggest a highly variable and idiosyncratic response to climate change [42]. Some populations are likely to persist, while others such as those in ROMO may be at high risk of extirpation [42]. American pikas are therefore considered a climate indicator species [43] and were petitioned to be listed under the Endangered Species Act in 2007. In 2010, the United States Fish and Wildlife Service (USFWS) concluded that listing of the species, or any portion of the species, was not currently warranted due to lack of scientific information [44]. Thus, understanding the factors that shape the distribution, genetic diversity, and adaptive potential for American pikas is of immediate conservation concern and underscores the need to explore innovative management strategies.
Photo credit: Dick Orleans.
Previous work based on morphology, dialect, and both mtDNA and nuclear coding sequence data identified a potential historic contact zone between the northern and southern Rocky Mountain lineages in the vicinity of ROMO [45–47]. Establishing the existence of such a contact zone within a single protected area would inform management of the species within ROMO and, potentially, encourage investigation of such contact zones for other species. Furthermore, as conservationists are increasingly exploring genetic rescue as a management option [9, 48], natural hybrid zones present an opportunity to evaluate concerns around anthropogenic manipulation of genetic structure. Previous phylogeographic studies of American pikas utilized a combination of molecular markers that reflect historical gene flow [47, 49], whereas the current study investigates contemporary gene flow among historically separated populations using markers with mutation rates that are relatively moderate (mtDNA) and high (microsatellite loci) as compared to nuclear coding sequences. Here we present a fine-scale population genetic study evaluating this potential contemporary contact zone and evidence for gene flow between the Northern Rocky Mountain (O. p. princeps) and Southern Rocky Mountain (O. p. saxatilis) lineages [50]. Given the spatial distribution of behavioral evidence for a hybrid zone [45], the American pika’s philopatric behavior and low dispersal ability [32, 51, 52], as well as the potential for extreme topography and streams to limit dispersal [53], we expect any gene flow among the two lineages to be limited to locations separated by no more than a few kilometers and shaped by landscape features.
Methods
Study sites and genetic sampling
In this study, we focus on ROMO as the context for contact between the two Rocky Mountain subspecies, in comparison with data from two additional study sites, Grand Teton National Park (GRTE) and Great Sand Dunes National Park (GRSA). GRTE and GRSA fall within the geographic range of the northern and southern Rocky Mountain lineages, respectively [47] (Fig 2). ROMO, however, does not fall within either subspecies range as defined in Galbreath et al. [47] (Fig 2). Detailed study design and microsatellite genotype data for these and other sites were previously reported for related studies [53]. Briefly, we collected fecal samples for genetic analyses between June 2010 and August 2014 through a combination of random, targeted, and opportunistic sampling. We collected random samples during standardized occupancy surveys conducted for other related studies [42, 54]. These survey locations were determined according to a generalized random stratified tessellation design [55] within potential pika habitat identified via remote sensing [54]. In addition to random sampling, we collected fecal samples opportunistically while traveling between survey locations, as well as through targeted searches of areas found to have pikas. We avoided collecting old fecal pellets by preferentially collecting pellets with green plant material inside to avoid degraded DNA. The color of the plant material fades from green to yellow within a few weeks to months [56]. We only collected distinct clusters of fecal pellets that were not contacting other previously deposited pellets to avoid contamination with DNA from other individuals. We collected fecal samples in paper coin envelopes and stored them dried until extraction. All field work was conducted as part of the Pikas in Peril? Project (PMIS #163377) and collections made under Scientific Research and Collection Permits (ROMO-2011-SCI-0032, GRSA-2010-SCI-0004, GRTE-2010-SCI-0079).
Study sites (numbered) and major mitochondrial lineages (black outlines) redrawn from Galbreath et al. 2010. Predicted distribution of American pikas (gray shading) derived from Kuchler potential natural vegetation [57, 58] according to Hafner and Sullivan [46]. Numbers refer to localities from Galbreath et al. 2010 and correspond to Table B in S1 File. Stars labeled with letters refer to sites from this study and correspond to Table B in S1 File. They are as follows: C) Grand Teton National Park, D) Rocky Mountain National Park north, F) Rocky Mountain National Park south, and G) Great Sand Dunes National Preserve.
Laboratory methods
We extracted genomic DNA from fecal samples using a modified AquaGenomic DNA extraction protocol (MultiTarget Pharmaceuticals LLC, Salt Lake City, UT, USA). We genotyped individuals at 24 microsatellite loci in four multiplex polymerase chain reactions (PCR) using a Qiagen Multiplex PCR kit (Qiagen, Valencia, CA, USA). Detailed PCR protocol, primer sequences, and methods for calling and screening microsatellite genotypes are provided in Castillo et al. [53]. In order to compare with previous phylogenetic analyses, we used the two primer pairs described in Galbreath et al. [49] to amplify the cytochrome-b oxidase (Cyt-b) and D-loop mtDNA region for a subset of our samples from ROMO, plus one from GRSA and GRTE. Due to the low quality and quantity of fecal DNA template as compared to fresh tissue, we designed an additional five primers to amplify smaller regions ranging from 483–566 bp (Table A in S1 File). All fragments were amplified in 10.5 μl reactions with final reagent concentrations of 2.25 mM MgCl2, 0.1 nM primers, 0.16 mM each dNTPs, 0.7 U Taq polymerase, and 0.5 μl template DNA. All reactions included a 15 minute (95°C) initial denaturation; 39 cycles of 30 sec. denaturation (95°C), 45 sec annealing (60°C), and 30 sec (72°C) extension; and a final 5 minute (72°C) extension. We sequenced all DNA fragments in both directions and visually inspected sequence alignments using Geneious 6.1.2 [59].
Genetic structure within Rocky Mountain National Park
We performed a Bayesian clustering analysis in program Structure [60] to infer population structure within ROMO. We ran 10 replicates for each inferred number of populations (K = 1 to 10) totaling 100 replicates, with 100,000 MCMC steps each of burnin and run steps. We used eight regions within the national park, identified by geographic features, as sampling localities for location prior information in the admixture model, as well as correlated allele frequencies among populations as run parameters. We analyzed the model output using Structure Harvester [61] according to the ΔK method proposed by Evanno et al. [62]. We compared the optimal K from both ΔK and mean Ln Pr(X|K) methods to determine the best K value. We then used program Clumpp [63] to determine the optimal assignment of individuals across the ten runs for the best supported value of K. Finally, we visualized the output from Clumpp spatially using ArcMap 10.0 (ESRI, Redlands, California).
Given the limitations of Structure and methods to interpret such analyses, particularly in the case of two genetic clusters [64], we also performed a principal components analysis (PCA) and discriminant analysis of principal components (DAPC) [65] using the package “adegenet” [66] in R [67]. DAPC maximizes variation between groups while minimizing variation within groups and has the benefit of not relying on assumptions of Hardy-Weinberg proportions [65]. We performed the DAPC on all individual genotypes with no prior population assignment information by using the successive K-means approach, implemented by the find.clusters function, to identify the optimal number of groups based on Bayesian Information Criterion (BIC).
Phylogenetic analysis
We obtained previously published sequence data covering the Cyt-b and D-loop region of the mitochondrial genome (c. 1650 bp) for 112 Ochotona princeps individuals, plus one O. collaris, from GenBank (Table B in S1 File). For Bayesian analysis, we selected the best data partitioning scheme and substitution models for each partition using a greedy algorithm in PartitionFinder 2.1.1 [68] with four a priori partitions for each codon of Cyt-b and one for D-Loop. We implemented the BIC corrected for small sample sizes to identify the best substitution model scheme (Table 1). We performed a Bayesian phylogenetic analysis of the partitioned matrix with MrBayes 3.1.2 [69] using the partitioning and models as detailed in Table 1. We performed two runs of four independent Markov chain Monte Carlo (MCMC) chains with 10M replicates each, sampling every 1000 generations or until the standard deviation of split frequencies between the two runs was less than 0.01. We discarded the first 25% of generations as burn-in for each run and then concatenated tree files. We found a maximum credibility tree (MCC), the tree with the highest product of posteriors for all nodes, using TreeAnnotator v1.8.2 [70]. We then edited the MCC tree in FigTree v.1.4.0 (http://tree.bio.ed.ac.uk/software/figtree/).
Population differentiation and genetic diversity
Once we identified two genetic clusters within ROMO (see results), we quantified genetic diversity and differentiation among the two groups, as well as GRTE and GRSA, from microsatellite genotypes. To measure deviation from panmixia assuming Hardy-Weingerg proportions, we calculated pairwise FST [71] using the “hierfstat” package [72] in R. Additionally, to measure population differentiation we calculated D [73, 74] using the “mmod” [75] package in R. We calculated geographic distance among study sites from the centroid of genotyped samples within each study site. We calculated expected heterozygosity (He) and allelic richness (Ar) corrected for sample size using the “hierfstat” package [72] in R. Finally, we calculated haplotype diversity (h) and nucleotide diversity (π) for the ROMO mtDNA dataset using DnaSP [76].
Geographic cline analysis and hybrid detection
We fit the admixture proportion values from the Structure analysis for K = 2 to equilibrium cline models using the Metropolis-Hastings Markov chain Monte Carlo algorithm implemented in the R package HZAR [77]. We fit 15 candidate models that varied in the number of cline shape parameters estimated and selected the best model according to AIC corrected for small sample size. We were thus able to estimate the geographic center and width of the cline. Finally, we performed a hybrid detection analysis to estimate the probability that individual genotypes reflect genotype frequency categories corresponding to pure individuals, F1 or F2 hybrids, or backcrosses. We categorized individuals with Q > 0.99 in either northern or southern cluster as “pure” individuals and categorized all others as of unknown origin. We performed the hybrid detection analysis using the program newhybrids [78], implemented in the R package parallelnewhybrid [78], with 100,000 burnin reps and 500,000 sweeps. We performed ten replicate runs and averaged the posterior probabilities for each individual across all ten runs.
Results
Genetic data
After removing individual samples that either failed to amplify, were contaminated (contained more than 2 microsatellite peaks at any locus), or gave inconsistent genotypes, our final dataset included 230 genotyped individuals from ROMO, 194 from GRTE, and 54 from GRSA. After screening loci for significant deviations from expected Hardy Weinberg proportions and removing loci that failed to amplify consistently across sites, we included 22 microsatellite loci. Number of alleles per locus across all four sites ranged from 5 to 28 (mean = 14).
Genetic structure
The Structure analysis supported two genetic clusters within ROMO (K = 2, Fig 3, Table A in S2 File and Figure A in S2 File). The DAPC likewise identified two clusters (Fig 4). All individuals were assigned to the same population clusters based on the DAPC (posterior membership probabilities > 0.99) and Structure analysis (Q ≥ 0.6, where Q is the proportion of the genome that originates from population K, also known as the admixture proportion [60]). Geographically, the two genetic clusters were roughly segregated into a northern and southern cluster. The northern cluster included individuals found north of Mt. Chapin or west of the Colorado River (Fig 3). The Structure analysis suggested some admixture occurred between the clusters (Fig 3). This was also supported by the PCA, where individuals identified as admixed in the Structure analysis (0.2 < Q < 0.8) had more intermediate values along the first principal component axis compared to individuals with higher Q values (Fig 5).
Individuals are shown as bar plots representing probability of assignment (q values) from the Structure analysis for K = 2. Individuals cluster geographically as a northern and southern population, separated by Mt. Chapin and the Colorado River. Concentric, black squares and circles indicate the placement of sequenced individuals into the northern and southern mtDNA lineages, respectively (Fig 6). Red shading indicates potential pika habitat [54]. Inset Figure shows ΔK values, indicating support for K = 2. Hillshade background was derived from the USGS National Elevation Dataset, streams and lakes were from the National Hydrography Dataset (https://nationalmap.gov).
Darker circles indicate multiple samples from that locality. Hillshade background was derived from the USGS National Elevation Dataset (https://nationalmap.gov).
Points are color-coded according to geographic location and correspond to sampling locations in inset map of ROMO. Triangles represent individuals identified from the STRUCTURE analysis as having notably mixed ancestry (0.2 ≤ Q < 0.8). Points in the dashed oval correspond to sampling localities within the dashed oval in the inset map, which in turn correspond to red points in Fig 4. These points separate along PC1 and are geographically isolated by topography. The sampling localities west of the Colorado River (dark green squares west of blue line in inset) partially separate along PC2.
We subsequently performed two separate Structure analyses for the individuals assigned to the northern and southern clusters with Q ≥ 0.6. Within the northern group, there was support for three clusters, with one primarily in the northeast and another mostly restricted to the west of the Colorado River (Table B in S2 File and Figure B in S2 File). This was supported by the PCA, where individuals west of the Colorado River segregated along the 2nd principal component axis (Fig 5). There was less support for population substructure within the southern group where there was some support for K = 2 and K = 6 (Table C in S2 File and Figure E in S2 File). Both results suggest isolation by distance with some restricted gene flow across streams and steep topography, but no major barriers to gene flow (Figure F in S2 File). Again, this was supported by a PCA of individuals within the southern cluster showing a latitudinal gradient along the first PCA (Figure G in S2 File).
Phylogenetic analysis
We included 19 individuals from within ROMO and one from GRSA that resulted in good quality sequence data covering the 1650 bp target region, along with the 113 individuals from Galbreath et al. (2009) (Table B in S1 File). PartitionFinder supported the a priori scheme with substitution models summarized in Table 1. The MrBayes runs achieved convergences within the first 10M generations. We recovered a MCC phylogeny with strong support (posterior probability > 0.9) for each of the five previously identified clades of O. princeps (Fig 6). Among the samples from ROMO, 9 individuals grouped within the northern Rocky Mountain lineage, O. p. princeps, and 10 individuals were grouped within the southern Rocky Mountain lineage, O. p. saxatilis (Fig 6). The results were consistent with both the Structure and DAPC analyses based on microsatellite genotypes (Fig 3).
Black circles on branches indicate Bayesian posterior probabilities >95%. Samples from ROMO appear in both the northern (blue) and southern (red) Rocky Mountain lineages. Numbers and letters in parentheses refer to Fig 2 and are listed in Table B in S1 File.
Population differentiation and genetic diversity
Pairwise FST calculated from microsatellite genotypes was greater between lineages than within lineages (Table 2) and increased with geographic distance (Fig 7). FST was greater between ROMO N and ROMO S, at a distance of <20 km, than within-lineages estimates hundreds of km apart, but less than between-lineage comparisons involving GRTE or GRSA. Pairwise D followed the same pattern (Table 2). Microsatellite allelic richness was greater within ROMO than either GRTE or GRSA, but these differences were not significant (-1.23 < t 1.38, p > 0.1). Likewise, microsatellite heterozygosity was similar among the study sites (Table 3) and relatively high compared to 9 other sites from a previous study [79]. Haplotype (h) and nucleotide diversity (π) calculated from mtDNA sequences was similar between ROMO N and ROMO S. We identified 7 haplotypes in each site. Genetic diversity was slightly higher in ROMO N (h = 0.00199, π = 0.92) than ROMO S (h = 0.00163, π = 0.91).
Pairs of sites within the same genetic lineage (filled circles) and between genetic lineages (open circles).
Geographic cline analysis and hybrid detection
The cline width for the best supported model was approximately 8.5 km with a two log-likelihood support range of 6.5–11.25 km (Table E in S2 File and Figure H in S2 File). The cline width estimates from the top 13 models were also within this range (Table E in S2 File). The estimated cline center was approximately 27.75 km north of the southernmost sampled individual, just south of Mt. Chapin (Table E in S2 File and Figure I in S2 File). Admixed individuals from the geographic north were characterized as either F2 hybrids or north backcrosses with posterior probability ≥ 0.6 (Figure J in S2 File). The one admixed individual from the geographic south was characterized as a south backcross with posterior probability 0.93. There was some support for backcrosses within the individuals with admixture proportions Q ≥ 0.8 (Figure J in S2 File).
Discussion
We identified a previously undescribed contact zone between the northern and southern Rocky Mountain lineages within Rocky Mountain National Park. Our results were consistent across all three different types of analyses (STRUCTURE, DAPC, and phylogenetic) and therefore robust to the K = 2 conundrum common in STRUCTURE analyses [64]. Moreover, we determined there was contemporary gene flow between the two lineages. Galbreath et al. [47] identified shared alleles among northern and southern Rocky Mountain populations in sequences of two nuclear introns (protein kinase C iota and mast cell growth factor) and determined that this was the result of gene flow since the last glacial maximum when receding montane glaciers were no longer barriers to dispersal. However, based on more-rapidly mutating mtDNA loci, they determined that the Colorado River represents a relatively recent barrier preventing contemporary gene flow between the two lineages [47]. In their taxonomic revision, Hafner and Smith [50] described the northern and southern subspecies’ ranges as occurring on either side of the Colorado River. The microsatellite markers used in this study have a considerably higher mutation rate than mtDNA, reflecting evolutionary processes within a few tens of generations [80–82], as compared to many hundreds to thousands of generations for mtDNA [82]. We determined that the Colorado River is not an adequate delineation of the geographic boundary between the two subspecies (Figs 2–4 and S2 File). However, the river does appear to be at least a partial barrier to gene flow (Fig 5 and S2 File).
Our results indicate contemporary gene flow between the two subspecies within ROMO, evidenced by intermediate levels of population differentiation as compared to estimates between populations from different or the same genetic lineage (Table 2), as well as the PCA and Structure analyses. However, admixture appears geographically limited to within a <10 km zone (Fig 3, Table E in S2 File, Figures H and I in S2 File). The geographic cline analysis should be interpreted with caution because 1) the area within and immediately surrounding the inferred hybrid zone was not exhaustively sampled, and 2) the analysis assumes a linear cline with minimal variation perpendicular to the cline [77]. We observed genetic structure on either side of the Colorado River (Figure C in S2 File), therefore this analysis should be repeated in the future with more extensive sampling and possibly omitting those individuals west of the Colorado River. Nevertheless, the results of the analysis were consistent with low dispersal ability further limited by geographic features observed in ROMO and other study areas [53]. Gaps in pika habitat, the Colorado River, and glacial valleys appear to contribute to genetic structure within (Figures C and E in S2 File) and among genetic lineages (Figs 2 and 3) despite their close geographic proximity. Alternately, the observed restricted gene flow among genetic lineages may be the result of only relatively recent contact, or possibly selection against hybrid individuals. This last hypothesis in particular warrants further investigation.
We did not find individuals with notably mixed ancestry in the southern cluster, with the exception of a single, likely backcrossed individual, collected just east of the continental divide (Fig 5). This pattern may be the result of more frequent dispersal of individuals from south to north. Previous research identified hybrid vocalizations at the headwaters of the Colorado River along the continental divide (Somers 1973), <5 kilometers from where we identified admixture among individuals. That study suggested that vocalizations were indicative of ancestry rather than learned behavior; and a subsequent study revealed that hybrid individuals reared in captivity do in fact produce hybrid vocalizations (Somers, personal communication). Additional observations of hybrid vocalizations outside ROMO suggest that further investigation should be made into the geographic extent and degree of admixture between these two genetic lineages. Dispersal is not sex-biased in pikas [83, 84], therefore differences in observed patterns between mtDNA and microsatellite analyses most likely do not reflect differences in dispersal between males and females. However, the sex of both parents and offspring may affect hybrid survival and fitness if for example traits are sex-linked or play a role in mate choice, such as in chemosensory behavior or vocalizations [14, 85–87]. Future research should investigate whether there is directional gene flow and if so, is it the result of dispersal patterns or selection for or against particular hybrid combinations.
Recent research has shown variation in population persistence across the species range over the past century [36, 88], as well variation in predicted future trends for American pikas that do not necessarily conform to a “colder is always better” scenario [36, 42, 54]. This may be in part driven by fine-scale microhabitat characteristics or other factors that are not captured by most models. Local adaptation and past biogeographic history likely also play important roles in population persistence such that populations that have experienced hotter climates in the past may be better prepared to deal with a warming climate in the future. We therefore might expect different responses to climate change between populations of the northern and southern Rocky Mountain lineages currently living under similar environmental conditions. What does this mean for hybrid lineages? Our research suggests, based on higher allelic richness, that gene flow between ROMO N and ROMO S has potentially increased genetic diversity within ROMO. In general, higher genetic diversity is thought to increase evolutionary potential in the face of rapid, environmental change [89]. However, gene flow may counteract local adaptation [90, 91]. Possible future scenarios for ROMO under a changing climate include: 1) resilience as a result of greater genetic diversity increasing adaptive potential through novel gene combinations, 2) decreased resilience as a result of the spread of maladapted genes from one lineage to the other, 3) increased resilience as a result of the spread of adaptive genes from one lineage to the other, or 4) little or no effect of admixture on resilience. Evidence from translocated bighorn sheep suggest adaptation to local environmental conditions was a strong determinant of translocation success in some cases [92], while other examples suggest high phenotypic plasticity and translocation success even within populations including subspecies hybrids [11]. Our study relied on neutral genetic markers to describe underlying genetic processes, therefore to address these hypotheses future research should seek to identify potential adaptive variation among these and other pika populations. Such studies could inform the feasibility of management actions such as translocating pikas within ROMO or among populations separated by greater geographic and environmental conditions, as has been proposed [93, 94].
The presence of two distinct genetic lineages as well a hybrid zone within a single park boundary presents some interesting, and potentially complicated, implications for management of this species within ROMO. Previous research demonstrated that there are five major American pika genetic lineages, each with independent evolutionary trajectories [49]. Galbreath et al. further suggest that these lineages should be considered distinct evolutionarily significant units [95] for management purposes. The 2010 decision by the USFWS to not list the American pika, or any subpopulations, under the ESA was based largely on subspecies revisions [50] according to these five independent lineages [44]. Should certain subspecies be listed in the future, pikas in ROMO could, in theory, be subject to different federal regulations, despite coexisting within a single management unit. One recent study suggested that pika in ROMO may be at particularly high risk of extirpation [42], but did not consider those subspecies separately. To further complicate this scenario, there are not clear guidelines under the ESA for treatment of hybrid individuals [96–98]. One example of successful translocations among subspecies and protection of hybrid offspring under the ESA is the introduction of Texas panthers (Puma concolor stanleyana) to populations of Florida panthers (P. c. coryi) [12]. In contrast, Allendorf et al. [99] recommended only pure Westslope cutthroat trout (Oncorhynchus clarki lewisi) be protected under the ESA and not hybrids with Yellowstone cutthroat trout (O. c. bouvieri) or rainbow trout (O. mykiss). While the ESA is a powerful management tool, it is often challenging to reconcile complex ecological and evolutionary processes with structured legal decisions.
Supporting information
S1 File. Mitochondrial DNA primer and sequence details.
https://doi.org/10.1371/journal.pone.0199032.s001
(DOCX)
S2 File. Results from population genetic Structure, hzar, and newhybrids analysis of American pikas within Rocky Mountain National Park, based on 22 microsatellite loci.
https://doi.org/10.1371/journal.pone.0199032.s002
(DOCX)
S3 File. Microsatellite genotypes and sampling localities.
https://doi.org/10.1371/journal.pone.0199032.s003
(CSV)
Acknowledgments
We thank all the National Parks Service staff, field technicians, and volunteers who contributed. We particularly thank R. Anthony, C.S. Baker, M. Baldner, M. Jeffress, D. Madison, C. Millar, A. Moffett, T. Rodhouse, D. Schwalm, W. Schweiger, and S. Wolff.
References
- 1.
Fisher RA. The genetical theory of natural selection. Oxford: Clarendon Press; 1930.
- 2. Hughes AR, Inouye BD, Johnson MTJ, Underwood N, Vellend M. Ecological consequences of genetic diversity. Ecology Letters. 2008;11(6):609–23. pmid:18400018
- 3. Pauls SU, Nowak C, Bálint M, Pfenninger M. The impact of global climate change on genetic diversity within populations and species. Molecular Ecology. 2013;22(4):925–46. pmid:23279006
- 4. Merilä J, Hendry AP. Climate change, adaptation, and phenotypic plasticity: the problem and the evidence. Evolutionary Applications. 2014;7(1):1–14. pmid:24454544
- 5. Eckert CG, Samis KE, Lougheed SC. Genetic variation across species’ geographical ranges: the central–marginal hypothesis and beyond. Molecular Ecology. 2008;17(5):1170–88. pmid:18302683
- 6. Hoffmann AA, Sgro CM. Climate change and evolutionary adaptation. Nature. 2011;470(7335):479–85. pmid:21350480
- 7. Hampe A, Petit RJ. Conserving biodiversity under climate change: the rear edge matters. Ecology Letters. 2005;8(5):461–7. pmid:21352449
- 8. Singer FJ, Papouchis CM, Symonds KK. Translocations as a tool for restoring populations of bighorn sheep. Restor Ecol. 2000;8(4):6–13. PubMed PMID: WOS:000166077600002.
- 9. Weeks AR, Sgro CM, Young AG, Frankham R, Mitchell NJ, Miller KA, et al. Assessing the benefits and risks of translocations in changing environments: a genetic perspective. Evolutionary Applications. 2011;4(6):709–25. pmid:22287981
- 10. Beschta RL, Ripple WJ. Riparian vegetation recovery in Yellowstone: The first two decades after wolf reintroduction. Biological Conservation. 2016;198:93–103. PubMed PMID: WOS:000377735400012.
- 11. Malaney JL, Feldman CR, Cox M, Wolff P, Wehausen JD, Matocq MD. Translocated to the fringe: genetic and niche variation in bighorn sheep of the Great Basin and northern Mojave deserts. Diversity and Distributions. 2015;21(9):1063–74. PubMed PMID: WOS:000360099700007.
- 12. Johnson WE, Onorato DP, Roelke ME, Land ED, Cunningham M, Belden RC, et al. Genetic Restoration of the Florida Panther. Science. 2010;329(5999):1641–5. pmid:20929847
- 13.
Allendorf F, Luikart G, Atiken S. Conservation and the Genetics of Populations. Blackwell Publishing. 2013.
- 14. Barton NH. The role of hybridization in evolution. Molecular Ecology. 2001;10(3):551–68. pmid:11298968
- 15. Bell G, Gonzalez A. Evolutionary rescue can prevent extinction following environmental change. Ecology Letters. 2009;12(9):942–8. pmid:19659574
- 16. Biebach I, Leigh DM, Sluzek K, Keller LF. Genetic Issues in Reintroduction. Reintroduction of Fish and Wildlife Populations. 2016:149–83. PubMed PMID: WOS:000388303700009.
- 17. Hewitt GM. Hybrid zones-natural laboratories for evolutionary studies. Trends in Ecology & Evolution. 1988;3(7):158–67. https://doi.org/10.1016/0169-5347(88)90033-X.
- 18. Gompert Z, Mandeville EG, Buerkle CA. Analysis of Population Genomic Data from Hybrid Zones. Annual Review of Ecology, Evolution, and Systematics. 2017;48(1):207–29.
- 19. Hewitt GM. Speciation, hybrid zones and phylogeography—or seeing genes in space and time. Molecular Ecology. 2001;10(3):537–49. pmid:11298967
- 20. Swenson Nathan G, Howard Daniel J. Clustering of Contact Zones, Hybrid Zones, and Phylogeographic Breaks in North America. The American Naturalist. 2005;166(5):581–91. pmid:16224723.
- 21. Hewitt GM. A sex-chromosome hybrid zone in the grasshopper Podisma pedestris (Orthoptera: Acrididae). Heredity. 1975;35:375. pmid:1061710
- 22. Anderson E. INTROGRESSIVE HYBRIDIZATION. Biological Reviews. 1953;28(3):280–307.
- 23.
Remington CL. Suture-Zones of Hybrid Interaction Between Recently Joined Biotas. In: Dobzhansky T, Hecht MK, Steere WC, editors. Evolutionary Biology: Volume 2. Boston, MA: Springer US; 1968. p. 321–428.
- 24. Thaeler JCS. Four Contacts Between Ranges of Different Chromosome Forms of the Thomomys Talpoides Complex (Rodentia: Geomyidae). Systematic Biology. 1974;23(3):343–54.
- 25. Barton NH, Hewitt GM. Analysis of hybrid zones. Annual Review of Ecology and Systematics. 1985;16:113–48. PubMed PMID: WOS:A1985AUL3900006.
- 26. Barton NH, Hewitt GM. Adaptation, speciation, and hybrid zones. Nature. 1989;341(6242):497–503. PubMed PMID: WOS:A1989AU72100044. pmid:2677747
- 27. Stewart JAE, Perrine JD, Nichols LB, Thorne JH, Millar CI, Goehring KE, et al. Revisiting the past to foretell the future: summer temperature and habitat area predict pika extirpations in California. Journal of Biogeography. 2015;42(5):880–90.
- 28. Jeffress MR, Van Gunst KJ, Millar CI. A surprising discovery of American pika sites in the northwestern Great Basin. Western North American Naturalist. 2017;77(2):252–68. PubMed PMID: WOS:000408018500012.
- 29. Beever EA, Brussard PE, Berger J. Patterns of apparent extirpation among isolated populations of pikas (Ochotona princeps) in the Great Basin. Journal of Mammalogy. 2003;84(1):37–54. PubMed PMID: ISI:000181414700004.
- 30. Beever EA, Dobrowski SZ, Long J, Mynsberge AR, Piekielek NB. Understanding relationships among abundance, extirpation, and climate at ecoregional scales. Ecology. 2013;94(7):1563–71. pmid:23951716
- 31. Stewart JAE, Wright DH, Heckman KA. Apparent climate-mediated loss and fragmentation of core habitat of the American pika in the Northern Sierra Nevada, California, USA. Plos One. 2017;12(8):17. PubMed PMID: WOS:000408693600008. pmid:28854268
- 32. Smith AT, Weston ML. Ochotona Princeps. Mammalian Species. 1990;352:1–8.
- 33. Smith AT. Distribution and dispersal of pikas: influences of behavior and climate. Ecology. 1974;55(6):1368–76. PubMed PMID: WOS:A1974V242000017.
- 34. Ivins BL, Smith AT. Responses of pikas (Ochotona princeps, Lagomorpha) to naturally-occurring terrestrial predators. Behavioral Ecology and Sociobiology. 1983;13(4):277–85. PubMed PMID: WOS:A1983RR47500006.
- 35. Holmes WG. Predator risk affects foraging behavior of pikas- observational and experimental evidence. Animal Behaviour. 1991;42:111–9. PubMed PMID: WOS:A1991FY40500012.
- 36. Millar CI, Heckman K, Swanston C, Schmidt K, Westfall RD, Delany DL. Radiocarbon dating of American pika fecal pellets provides insights into population extirpations and climate refugia. Ecological Applications. 2014.
- 37. Millar CI, Westfall RD. Distribution and climatic relationships of the American pika (Ochotona princeps) in the Sierra Nevada and western Great Basin, USA; periglacial landforms as refugia in warming climates. Arctic Antarctic and Alpine Research. 2010;42(1):76–88. PubMed PMID: ISI:000275560500008.
- 38. Collins GH, Bauman BT. Distribution of Low-Elevation American Pika Populations in the Northern Great Basin. Journal of Fish and Wildlife Management. 2012;3(2):311–8.
- 39. Millar CI, Westfall RD, Delany DL. New records of marginal locations for American pika (Ochotona princeps) in the western Great Basin. Western North American Naturalist. 2013;73(4):457–76.
- 40. Varner J, Dearing MD. The Importance of Biologically Relevant Microclimates in Habitat Suitability Assessments. PLoS ONE. 2014;9(8):e104648. pmid:25115894
- 41. Calkins MT, Beever EA, Boykin KG, Frey JK, Andersen MC. Not-so-splendid isolation: modeling climate-mediated range collapse of a montane mammal Ochotona princeps across numerous ecoregions. Ecography. 2012;35(9):780–91. PubMed PMID: WOS:000305941800002.
- 42. Schwalm D, Epps CW, Rodhouse TJ, Monahan WB, Castillo JA, Ray C, et al. Habitat availability and gene flow influence diverging local population trajectories under scenarios of climate change: a place-based approach. Global Change Biology. 2016;22(4):1572–84. pmid:26667878
- 43. Hafner DJ. North-American pika (Ochotona princeps) as a late Quaternary biogeographic indicator species. Quaternary Research. 1993;39(3):373–80. PubMed PMID: ISI:A1993LA95200012.
- 44.
USFWS. 12-month Finding on a Petition to List the American Pika as Threatened or Endangered. Federal Register2010. p. 6437–71.
- 45. Somers P. Dialects in southern Rocky Mountain pikas, Ochotona princeps (Lagomorpha). Animal Behaviour. 1973;21(1):124–37. http://dx.doi.org/10.1016/S0003-3472(73)80050-8.
- 46. Hafner DJ, Sullivan RM. Historical and ecological biogeography of Nearctic pikas (Lagomorpha, Ochotonidae). Journal of Mammalogy. 1995;76(2):302–21. PubMed PMID: ISI:A1995QZ50200003.
- 47. Galbreath KE, Hafner DJ, Zamudio KR, Agnew K. Isolation and introgression in the Intermountain West: contrasting gene genealogies reveal the complex biogeographic history of the American pika (Ochotona princeps). Journal of Biogeography. 2010;37(2):344–62. PubMed PMID: ISI:000273771100012.
- 48. Love Stowell SM, Pinzone CA, Martin AP. Overcoming barriers to active interventions for genetic diversity. Biodiversity and Conservation. 2017;26(8):1753–65.
- 49. Galbreath KE, Hafner DJ, Zamudio KR. When cold is better: Climate-driven elevation shifts yield complex patterns of diversification and demography in an alpine specialist (American pika, Ochotona princeps). Evolution. 2009;63(11):2848–63. PubMed PMID: ISI:000271031900007. pmid:19663994
- 50. Hafner DJ, Smith AT. Revision of the subspecies of the American pika, Ochotona princeps (Lagomorpha: Ochotonidae). Journal of Mammalogy. 2010;91(2):401–17. PubMed PMID: ISI:000277120600011.
- 51. Peacock MM. Determining natal dispersal patterns in a population of North American pikas (Ochotona princeps) using direct mark-resight and indirect genetic methods. Behavioral Ecology. 1997;8(3):340–50. PubMed PMID: ISI:A1997XB63800013.
- 52. Castillo JA, Epps CW, Davis AR, Cushman SA. Landscape effects on gene flow for a climate-sensitive montane species, the American pika. Molecular Ecology. 2014;23(4):843–56. PubMed PMID: WOS:000330264000009. pmid:24383818
- 53. Castillo JA, Epps CW, Jeffress MR, Ray C, Rodhouse TJ, Schwalm D. Replicated landscape genetic and network analyses reveal wide variation in functional connectivity for American pikas. Ecological Applications. 2016;26(6):1660–76. pmid:27755691
- 54. Jeffress MR, Rodhouse TJ, Ray C, Wolff S, Epps CW. The idiosyncrasies of place: geographic variation in the climate-distribution relationships of the American pika. Ecological Applications. 2013;23(4):864–78. pmid:23865236
- 55. Stevens DL, Olsen AR. Spatially balanced sampling of natural resources. Journal of the American Statistical Association. 2004;99(465):262–78. PubMed PMID: WOS:000220638200025.
- 56. Nichols LB. Fecal pellets of American pikas (Ochotona princeps) provide a crude chronometer for dating patch occupancy. Western North American Naturalist. 2010;70(4):500–7. PubMed PMID: ISI:000286649600010.
- 57. Kuchler AM. Potential natural vegetation of the conterminous United States. Am Geogr Soc, Spec Publ. 1964;36.
- 58.
Conservation Biology Institute. U.S. potential natural vegetation, original Kuchler types, v2.0 (spatially adjusted to correct geometric distortions). 2012. Available: http://databasin.org/datasets/1c7a301c8e6843f2b4fe63fdb3a9fe39
- 59. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28(12):1647–9. pmid:22543367
- 60. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–59. PubMed PMID: ISI:000087475100039. pmid:10835412
- 61. Earl D, vonHoldt B. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genet Resour. 2012;4(2):359–61.
- 62. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology. 2005;14(8):2611–20. PubMed PMID: WOS:000229961500029. pmid:15969739
- 63. 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(14):1801–6. pmid:17485429
- 64. Janes JK, Miller JM, Dupuis JR, Malenfant RM, Gorrell JC, Cullingham CI, et al. The K = 2 conundrum. Molecular Ecology. 2017;26(14):3594–602. pmid:28544181
- 65. Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC genetics. 2010;11(1):94.
- 66. Jombart T. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24(11):1403–5. pmid:18397895
- 67.
Team RC. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2015.
- 68. Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. PartitionFinder 2: New Methods for Selecting Partitioned Models of Evolution for Molecular and Morphological Phylogenetic Analyses. Molecular Biology and Evolution. 2017;34(3):772–3. pmid:28013191
- 69. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: Efficient Bayesian Phylogenetic Inference and Model Choice Across a Large Model Space. Systematic Biology. 2012;61(3):539–42. pmid:22357727
- 70.
Rambaut A, Drummond A. TreeAnnotator v1. 7.0. 2013. Available as part of the BEAST package at http://beast.bio.ed.ac.uk
- 71. Weir BS, Cockerham CC. Estimating F-Statistics for the Analysis of Population Structure. Evolution. 1984;38(6):1358–70. pmid:28563791
- 72. Goudet J. Hierfstat, a package for r to compute and test hierarchical F-statistics. Molecular Ecology Notes. 2005;5(1):184–6.
- 73. Jost L. GST and its relatives do not measure differentiation. Molecular Ecology. 2008;17(18):4015–26. pmid:19238703
- 74. Whitlock MC. and D do not replace FST. Molecular Ecology. 2011;20(6):1083–91. pmid:21375616
- 75. Winter DJ. mmod: an R library for the calculation of population differentiation statistics. Molecular Ecology Resources. 2012;12(6):1158–60. pmid:22883857
- 76. Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, et al. DnaSP 6: DNA Sequence Polymorphism Analysis of Large Datasets. Molecular Biology and Evolution. 2017:msx248-msx. pmid:29029172
- 77. Derryberry EP, Derryberry GE, Maley JM, Brumfield RT. hzar: hybrid zone analysis using an R software package. Molecular Ecology Resources. 2014;14(3):652–63. pmid:24373504
- 78. Wringe BF, Stanley RRE, Jeffery NW, Anderson EC, Bradbury IR. parallelnewhybrid: an R package for the parallelization of hybrid detection using newhybrids. Molecular Ecology Resources. 2017;17(1):91–5. pmid:27617417
- 79.
Castillo JA. Population genetics of American pika (Ochotona princeps): Investigating gene flow and genetic diversity across multiple, complex landscapes. Corvallis, OR: Oregon State University; 2015.
- 80. Spear SF, Storfer A. Landscape genetic structure of coastal tailed frogs (Ascaphus truei) in protected vs. managed forests. Molecular Ecology. 2008;17(21):4642–56. PubMed PMID: ISI:000260345200007. pmid:19140987
- 81. Landguth EL, Cushman SA, Schwartz MK, McKelvey KS, Murphy M, Luikart G. Quantifying the lag time to detect barriers in landscape genetics. Molecular Ecology. 2010;19(19):4179–91. pmid:20819159
- 82. Wang IJ. Choosing appropriate genetic markers and analytical methods for testing landscape genetic hypotheses. Molecular Ecology. 2011;20(12):2480–2.
- 83. Smith AT, Ivins BL. Colonization in a pika population—dispersal vs philopatry. Behavioral Ecology and Sociobiology. 1983;13(1):37–47. PubMed PMID: WOS:A1983RA45800006.
- 84. Zgurski JM, Hik DS. Polygynandry and even-sexed dispersal in a population of collared pikas, Ochotona collaris. Animal Behaviour. 2012;83(4):1075–82.
- 85. Smadja C, Butlin RK. On the scent of speciation: the chemosensory system and its role in premating isolation. Heredity. 2008;102:77. pmid:18685572
- 86. Gerhardt HC. Sound Pattern Recognition in Some North American Treefrogs (Anura: Hylidae): Implications for Mate Choice1. American Zoologist. 1982;22(3):581–95.
- 87. den Hartog PM, de Kort SR, ten Cate C. Hybrid vocalizations are effective within, but not outside, an avian hybrid zone. Behavioral Ecology. 2007;18(3):608–14.
- 88. Erb LP, Ray C, Guralnick R. On the generality of a climate-mediated shift in the distribution of the American pika (Ochotona princeps). Ecology. 2011;92(9):1730–5. PubMed PMID: WOS:000294152800004. pmid:21939069
- 89. Frankham R. Genetics and extinction. Biological Conservation. 2005;126(2):131–40. http://dx.doi.org/10.1016/j.biocon.2005.05.002.
- 90. Funk WC, McKay JK, Hohenlohe PA, Allendorf FW. Harnessing genomics for delineating conservation units. Trends in Ecology & Evolution. 2012;27(9):489–96. http://dx.doi.org/10.1016/j.tree.2012.05.012.
- 91. Slatkin M. Gene Flow in Natural Populations. Annual Review of Ecology and Systematics. 1985;16(1):393–430.
- 92. Wiedmann BP, Sargeant GA. Ecotypic Variation in Recruitment of Reintroduced Bighorn Sheep: Implications for Translocation. Journal of Wildlife Management. 2014;78(3):394–401. PubMed PMID: WOS:000336027400003.
- 93. Mathewson PD, Moyer-Horner L, Beever EA, Briscoe NJ, Kearney M, Yahn JM, et al. Mechanistic variables can enhance predictive models of endotherm distributions: the American pika under current, past, and future climates. Global Change Biology. 2017;23(3):1048–64. pmid:27500587
- 94. Wilkening JL, Ray C, Ramsay N, Klingler K. Alpine biodiversity and assisted migration: the case of the American pika (Ochotona princeps). Biodiversity. 2015:1–13.
- 95. Moritz C. Defining ‘Evolutionarily Significant Units’ for conservation. Trends in Ecology & Evolution. 1994;9(10):373–5. http://dx.doi.org/10.1016/0169-5347(94)90057-4.
- 96.
Haig SM, Allendorf FW. Hybrids in Policy. In: Scott JM, Goble DD, Davis F, editors. The Endangered Species Act at Thirty: Conserving Biodiversity in Human-dominated Landscapes. 2. Washington, DC: Island Press; 2006. p. 150–63.
- 97. vonHoldt BM, Brzeski KE, Wilcove DS, Rutledge LY. Redefining the Role of Admixture and Genomics in Species Conservation. Conservation Letters. 2017.
- 98. Ellstrand NC, Biggs D, Kaus A, Lubinsky P, McDade LA, Preston K, et al. Got Hybridization? A Multidisciplinary Approach for Informing Science Policy. BioScience. 2010;60(5):384–8.
- 99. Allendorf FW, Leary RF, Hitt NP, Knudsen KL, Lundquist LL, Spruell P. Intercrosses and the U.S. Endangered Species Act: Should Hybridized Populations be Included as Westslope Cutthroat Trout? Conservation Biology. 2004;18(5):1203–13.