Life on the edge - a changing genetic landscape within an iconic American pika metapopulation over the last half century

doi:10.21203/rs.3.rs-1735950/v1

Download PDF

Research Article

Life on the edge - a changing genetic landscape within an iconic American pika metapopulation over the last half century

https://doi.org/10.21203/rs.3.rs-1735950/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 27 Sep, 2023

Read the published version in PeerJ →

Version 1

posted

You are reading this latest preprint version

Declines and extirpations of American pika (Ochotona princeps) populations at historically occupied sites started being documented in the literature during the early 2000s. Commensurate with global climate change, many of these losses at peripheral and lower elevation sites have been associated with changes in ambient air temperature and precipitation regimes. Here, we report on a decline in available genetic resources for an iconic American pika metapopulation, located at the southwestern edge of the species distribution in the Bodie Hills of eastern California, USA. Composed of highly fragmented habitat created by hard rock mining, the ore dumps at this site were likely colonized by pikas at the end of the 19th century from nearby natural talus outcrops. Genetic data extracted from both contemporary samples and archived natural history collections allowed us to track population and patch-level genetic diversity for Bodie pikas across three distinct sampling points during the last half- century (1948, 1988–1991, 2013–2015). In addition to declines in within-population allelic diversity and expected heterozygosity, we observed an increase in population structure and a reduction in effective population size from more extensive sampling of extant patches during 1988–1991 and 2013–2015, respectively. Furthermore, census records from the last 50 years as well as archived museum samples collected in 1947 from a nearby pika population in the Wassuk range (NV, USA) provide further support of the increasing isolation and genetic coalescence occurring in this region. This study highlights the importance of museum samples and long-term monitoring in contextualizing our understanding of population viability.

American pika

climate change

genetic monitoring

evolutionary potential

conservation science

museum specimens

The spatial distribution of an organism’s habitat can vary at local, regional, and global scales (Shen et al. 2013: Anderson et al. 2016; Klingler, Jahner et al. 2021). How an organism navigates the spatial structure of its habitat will depend upon its body size, mode of locomotion and life history (Peacock and Smith 1997; Pinto and MacDougall 2010; Wauters et al. 2010; Neam and Lacher 2018). Climate change is likely to alter not only the extent of suitable habitat, but also the spatial organization of occupiable habitat patches for many organisms. As some habitat patches become unsuitable, the distance among inhabitable patches may increase, potentially impeding an organism’s ability to move among otherwise suitable habitat areas. Such a dynamic may spiral into an extinction vortex (Benson et al. 2019) as populations or subpopulations become increasing isolated and subject to demographic and genetic stochasticity (Lande 1993), thereby reducing both persistence probability and evolutionary potential (Blomqvist et al. 2010).

Here we examine the temporal dynamics of patch occupancy and the maintenance of genetic diversity within a population of American pika (Ochotona princeps), a small, thermally sensitive alpine lagomorph considered vulnerable to global warming (Beever et al. 2010; Calkins et al. 2012; Stewart et al. 2015; Schwalm et al. 2016). O. princeps is broadly distributed across mountain ranges of the western United States and Canada (Smith and Weston 1990; Hafner and Smith 2010). Pikas are known to occupy a diversity of rocky habitats including talus patches, lava flows and even hard-rock mining ore dumps (Severaid 1955; Peacock 1997; Peacock and Smith 1997; Rodhouse et al. 2010). Such habitat is naturally fragmented and varies widely in patch size, proximity to other patches and therefore levels of connectivity across the range of this species. Because pikas already sustain a high resting body temperature (40.1°C) very close to their upper lethal temperature (43.1°C), the thermal consistency of the interstitial spaces of rocky talus habitat allows them to adjust their internal temperature through behavioral thermoregulation (MacArther and Wang 1973; Wilkening et al. 2015). This physiological sensitivity combined with a limited dispersal capacity (Smith and Ivins 1983; Peacock 1997) increases pika vulnerability to changing environmental conditions under climate change (Stewart et al. 2017).

Over the past several decades, pika populations have been disappearing, particularly in the mountain ranges of the Great Basin physiographic region, southern Utah, and lower elevation sites on the eastern slopes of the Sierra Nevada (Beever et al. 2003; 2010; 2011; 2016; Stewart and Wright 2012; Nichols et al. 2016; Stewart et al. 2015; 2017). These declines disproportionately affect the subspecies O. princeps schisticeps, whose distribution includes the Sierra Nevada and Great Basin mountain ranges of western and central Nevada (Galbreath et al. 2009; 2010). In addition to evidence of shifts in the elevational distribution of pikas as a result of prehistoric climate change (late Pleistocene (Wisconsinan) to early Holocene; Hafner 1994; Grayson 2005), metrics of temperature measured during 20th and 21st century human mediated climate change have repeatedly been identified as predictors of extant pika distributional change and persistence (Beever et al. 2003; Beever et al. 2010; Calkins et al. 2012; Stewart et al. 2015; Schwalm et al. 2016). However, at smaller spatial scales, non-climatic factors such as habitat configuration, connectivity, and suitable dispersal habitat are also likely to contribute to site-specific persistence probabilities that vary across the species range (Peacock and Smith 1997; Jeffress et al. 2013; Castillo et al. 2014; Schwalm et al. 2016). Indeed, pika populations continue to persist in some low-elevation sites (Beever et al. 2008; Rodhouse et al. 2010; Manning and Hagar 2011; Jeffress et al. 2017), which may reflect the nuances of local habitat suitability and potential for behavioral plasticity in this species (Erb et al. 2011; Varner and Dearing 2014; Moyer-Horner et al. 2015; Mathewson et al. 2017).

This study examines changes in patch occupancy, levels of genetic diversity, and dispersion of genetic variation across the landscape for an iconic pika metapopulation in the Bodie Hills (Mono County, California) (Severaid 1955; Smith 1974; Peacock and Smith 1997), a low elevational spur of the main Sierra Nevada mountain range on the border with Nevada. Stewart et al. (2015) recently investigated climate change and pika persistence in the Sierra Nevada, which included the Bodie population, and showed that average summer temperature and talus extent within a one km radius were predictors of both current occupancy and extirpation for historically occupied sites in the Sierra Nevada. Now protected as a California State Historic Park (BSHP) (elevation, ~ 2500 m), the BSHP population is one of the only extant pika populations remaining within the Bodie Hills as almost all the natural rock outcroppings found here are currently unoccupied by pikas despite evidence of recent 20th and 21st century occupation (Nichols 2011; Stewart et al. 2015; Nichols et al. 2016).

Embedded within Great Basin sagebrush-scrub plant community, pika habitat at BSHP is composed of ore dumps created by hard-rock mining during the late 19th and early 20th centuries, and pikas are believed to have colonized these ore dumps from neighboring natural talus in the early 1900s (Fig. 1; Severaid 1955). The pika population at BSHP was first reported by J.H. Severaid (Severaid 1955) and has been under continuous study since 1969 (Smith 1974a, 1974b, 1978; Peacock and Smith 1997; Moilanen et al. 1998; Nichols 2011; Stewart et al. 2015; Nichols et al. 2016; Klingler, Jahner et al. 2021). This population has been the subject of studies on metapopulation dynamics (Peacock and Smith 1997; Moilanen et al. 1998), population genetic structure (Peacock and Smith 1997), species distribution modeling (Stewart et al. 2015) and most recently a range wide comparison of genomic variation among multiple lineages and subspecies (Klingler, Jahner et al. 2021).

We use the results of two early census periods in 1972 and 1977, together with 30 years of annual BSHP population census survey data (Moilanen et al. 1998; Nichols, this study) as well as contemporary population genetic data (1988–1991 and 2013–2015), to test the hypothesis that genetic diversity has eroded at this site due in part to changes in patch occupancy likely attributable to 20th and 21st century climate change. Furthermore, we provide historical context of the population genetic diversity once present at this site by comparing these two sets of more contemporary genetic samples to museum samples collected at BSHP in 1948 and 1949 by J.H. Severaid and curated at the Museum of Vertebrate Zoology (MVZ), University of California at Berkeley. In addition, we also provide a regional context by comparing BSHP samples to historic samples collected at Big Indian Mountain (BIH), Wassuk Range, Mineral County Nevada, ~ 35 km northeast of BSHP, by S.B. Benson and O.P. Pearson in 1947 (curated at MVZ).

Recent census data estimates the extant population to be < 100 individuals, representing about 100 pika generations since colonization from natural outcrops circa 1900. Ore dump patches that had been previously occupied and were therefore considered suitable habitat in the southern and middle sections of BSHP experienced a population collapse sometime in the 1980s which was later documented by the 1989 census (Fig. 1; Moilanen et al. 1998). The southern half of the study area has never fully recovered with only a few of the extirpated talus patches being intermittently recolonized by the extant population in the northern half of the mining district (Moilanen et al. 1998; Nichols, this study). The increasing isolation of the BSHP pika population within the Bodie Hills during recent years (Nichols et al. 2016) occurred simultaneously with severe drought (2012–2014) in California that led to increasingly warm air temperatures and low levels of precipitation (Griffin and Anchukaitis 2014; Diaz and Wahl 2015). The frequency of drought years during the last two decades has increased compared to the previous 100 years and the likelihood of extremely warm, low precipitation years is projected for the near future (AghaKouchak et al. 2014; Diffenbaugh et al. 2015; Williams et al. 2015).

We hypothesize that the warming air temperatures and decreased precipitation documented for this region have negatively affected survivorship, recruitment and patch occupancy of pikas at BSHP over the latter half of the 20th and early 21st centuries. We also hypothesize that changes in patch occupancy, in terms of both the duration of occupation and turnover dynamics, have fundamentally altered the metapopulation dynamic of this site and led to genetic coalescence and concomitant losses in genetic diversity (Gilpin 1991). Therefore, we expect to observe evidence of a disrupted metapopulation dynamic (Peacock and Smith 1997) that includes: (1) a steady reduction in the proportion of patches occupied (as measured by proportion of patches occupied/total patches surveyed), (2) evidence of genetic coalescence and genetic erosion via a decline in allelic diversity and heterozygosity, (3) a decrease in the effective population size (Ne) over the 20th and 21st centuries and (4) increasing regional isolation.

Study species

Pikas are small alpine lagomorphs (~ 150 g), which inhabit rocky substrates (talus) in the mountain ranges of the American west. Adult pikas are individually territorial and defend their territories from conspecific intrusions (Smith and Ivins 1983; Peacock 1997). As lagomorphs, they do not hibernate, and store vegetation collected during the summer months in one or two centrally placed “haypiles” on their territories that they subsequently feed on during the winter months (Smith and Weston 1990). Haypile sites, which define the effective center of an adult's home range, are maintained in the same location even when territory ownership changes (Peacock 1997). Territories are found at the talus-vegetation interface where there is easy access to forage and escape cover from predators. Suitable locations for territories on a talus slope are limited and therefore the number of territories found in any one location is relatively static. Juveniles are born in late spring/early summer and remain on their natal territory until they are weaned at which time they begin exploratory movements to locate open territories (Smith and Ivins 1983). Territory ownership is essential for overwinter survival and juveniles typically colonize the first available territory they encounter during their natal summer (Peacock 1997; Peacock and Smith 1997).

BSHP Population census surveys

The majority (~ 80 patches) of the total suitable habitat patches at BSHP were surveyed in 1972, 1977, 1989, 1991–2001, 2003–2006, and 2008–2018 (Smith 1974a; Moilanen et al. 1998; Nichols, this study). In order to characterize patch occupancy over time as well as provide a crude estimate of annual population size, we calculated: 1) the number of active territories per habitat patch per year determined by presence of fresh fecal pellets in latrine sites and/or active haypiles indicated by fresh vegetation (Dearing 1997); 2) the proportion of patches occupied for each survey period; and 3) a census size for the entire patch network for each survey period based upon individual patch census data. We then conducted linear regression analyses (SYSTAT 1998) to assess the relationship between the proportion of patches occupied per survey year over all years as well as changes in census size for the entire patch network as well as the southern (all patches south of Silver Hill), middle (patches on Silver Hill south of Mono Shaft) and northern (all patches north of and including Mono Shaft) patch networks separately over time (see Fig. 1).

Genetic sample collection

Mid-20th century museum samples (1947–1949). Approximately 200 museum specimens collected during the mid-20th century at BSHP by J.H. Severaid (1955) are curated in the mammal collections in the Museum of Vertebrate Zoology, University of California at Berkeley. Additional museum preserved samples are available from a number of historically occupied pika sites across Nevada. K. Klingler was given access to 16 specimens from BSHP collected by Severaid 1948–1951 (Supplementary material, Table S1). E. Hekkala was given access to an additional 19 samples collected from BHSP (1948–1949) by Severaid and Big Indian Mountain, Wassuk Range, Nevada (BIH, 1947) collected by S.B. Benson and O.P. Pearson (Table S1). A 2–5 mm square piece of dried museum pelt was cut from the ventral lip and placed in a 1.5 ml locking tube. We included the BIH historical samples in order to test for a pattern of increasing regional isolation of the BSHP population over the timeframe from 1948–2015.

Late 20th century samples (1988–1990). Pika ear tissue was collected from live-trapped individuals in BSHP as part of a metapopulation dynamics study of 16 occupied patches in the northern portion of the study area (1988–1990; Peacock and Smith 1997). After completion of the study, the remaining DNA samples were stored in a -20°C freezer until present (n = 114). All adult samples (n = 75) from the 1988–1991 sampling period were used to estimate population level genetic metrics for this study.

Early 21st century samples (2013–2015). Fecal DNA samples were collected from 26 occupied habitat patches in the northern half of the BSHP site during June-September 2013–2015 to be used for comparison to the earlier time periods. Adult pikas defecate in well-defined latrine sites within their territories (Nichols 2011). Differences in size between adult and juvenile animals make identification of adult versus juvenile fecal pellets unambiguous. Multiple fresh fecal pellets (~ 6–8) from single latrine sites identified using methods of Nichols (2011), were collected per individual using stainless steel serrated tipped forceps. Any cellular contamination among samples was minimized by wiping down the forceps after use, then submerging them in 70–100% ethanol and allowing the ethanol to volatilize by air drying between sample collections. Two to five fecal pellets per individual animal were immediately transferred to a 2.0 ml vial with 500 µl of ATL tissue lysis buffer (QIAGEN DNeasy96 Blood and Tissue kits) after collection in the field. At the end of each collection period (~ 2–4 hours), 50 µl of ice-cold proteinase K was added to 400 µl of supernatant from each sample for storage until DNA extraction. We did not homogenize the pellets, but inverted them multiple times once placed in the lysis buffer in order to wash the epithelial cells from the pellet surface prior to pipetting the supernatant and adding the proteinase K.

DNA isolation and PCR amplification

Total genomic DNA was isolated according to the manufacturer’s protocol (DNeasy96 Blood and Tissue kits, QIAGEN) with modifications specific to each sample type. Details regarding modifications to DNA extraction protocol for fecal and museum genetic samples can be found in Peacock et al. (2017) and Hekkala et al. (2011). All DNA extractions and polymerase chain reactions (PCR) for early 20th century museum samples (1947–1949) were conducted in a separate laboratory in order to avoid DNA contamination from contemporary pika samples that have been isolated in our laboratory.

We used eight di, tri, and tetranucleotide repeat, polymorphic microsatellite loci of the 28 that have been developed for the American pika (OCP4, OCP6, OCP7, OCP8, OCP9, Peacock et al. 2002; OCP10, OCP17, OCP21; Peacock and Kirchoff 2009, direct submission Genbank). We chose microsatellite loci that were smaller (in terms of fragment length) and that produced a consistent and scorable PCR product from both the fecal and museum DNA samples. PCRs were carried out in 12 µl reaction volumes on a MBS Satellite 0.2G thermal cycler (Thermo Electron Corporation), and included a 15 minute hot start at 94°C, followed by 41 cycles at 94°C, 30s each for denaturing, a touch down annealing temperature for 90 seconds (seven cycles of 65°C, followed by seven cycles of 61°C, and seven cycles of 58°C each were completed with a final 20 cycles at 55°C) and an elongation step at 72°C for 30 seconds (modified from Peacock et al. 2002). The first 21 cycles amplified the locus specific primer and the final 20 cycles added a fluorescently labeled M13 tail to the PCR product. All PCR products were diluted to an appropriate concentration for use with an automated sequencer and 1 µl of PCR product was added to 19 µl of HiDye Formamide/LIZ500 size standard for fragment analysis (Applied Biosystems). Fragment size analysis was carried out on an Applied Biosystems 3730 DNA Genetic Analyzer at the UNR Nevada Genomics Center (http://www.ag.unr.edu/genomics/). All alleles were scored, binned, and genotyped using ABI Prism GeneMapper software (version 5.0).

Allelic dropout and genotyping error rates for fecal and museum samples

DNA from low quality fecal and preserved museum samples can lead to genotyping errors, which may bias estimates of population genetic metrics (Creel et al. 2003; McKelvey and Schwartz 2005). Common genotyping errors include allelic dropout, preferential allele amplification, and false alleles. For contemporary fecal samples, each individual was genotyped multiple times (2–4) at all loci using separate PCR reactions until a consensus genotype was reached. Individuals were only included in genetic analyses if a consensus genotype was achieved through repeated, successful amplification. For museum samples, DNA was amplified in five separate PCR reactions for all loci and individuals were included in genetic analyses only if a consensus genotype was reached at each locus for three of the five PCR reactions.

For both the fecal and museum DNA samples, we did not score any questionable low amplitude alleles and eliminated individuals that had three or more scorable alleles at any locus from the analysis. We tested for both heterozygote excess and deficiency per locus per patch subpopulation for the fecal DNA samples using the Fisher Exact test in GENEPOP (version 4.2; Rousset 2008) in order to rule out multiple individuals in DNA isolated from multiple fecal pellets and to also control for possible contamination during pellet collection. Once consensus genotypes were obtained, we identified and removed duplicate samples from the 2013–2015 fecal DNA dataset using the Multilocus Matches DS module in GenAlEx (version 6.501; Peakall and Smouse 2006; 2012). In addition, we used the probability-of-identify function in GenAlEx for the 2013–2015 data to assess our power to detect individuals using the 8 microsatellite loci dataset. As the 1988–1991 data were tissue samples obtained from individuals who were individually marked with ear tags, individual identity was not an issue.

Locus-specific error rates for allelic dropout (ε₁), false alleles (ε₂), and genotyping success rate were estimated using replicate genotypes for fecal DNA and museum specimens through a maximum-likelihood based approach (program Pedant, Delphi version 7.0) (Johnson and Haydon, 2007). In addition, the programs MicroChecker (version 2.2.3; Van Oosterhout and Hutchinson 2004) and FreeNA (Chapuis and Estoup 2007) were used to test for the presence of null alleles, and genotyping errors. If a locus exhibited systematic patterns of deviation from HWE, showed evidence of null alleles or allelic dropout it was removed from the analysis.

Population genetic analyses

We used FSTAT (version 2.9.3.2; Goudet et al. 2002) to quantify the number of alleles (Na), to test for linkage disequilibrium and deviations from Hardy-Weinberg equilibrium (F_IS) per locus per time period and to calculate a global θ for each time period (1948–1949, 1988–1991, 2013–2015). We also used FSTAT to calculate pairwise F_ST estimates for individual patches sampled in both 1988–1991 and 2013–2015 that had sufficient data for comparisons (patch occupancy n ≥ 3 individuals). We also calculated a pairwise F_ST at the BSHP population scale for comparison between the two time periods (1988–1991 and 2013–2015) that had data from extensive sampling.

In order to correct for differences in sample size between the three sampling periods, a hierarchical rarefaction approach was used in the program HP-RARE (version 1.0; Kalinowski 2005) to calculate estimates of allelic richness (Rs) and private alleles (Pa). The average observed (Ho) and expected (He) heterozygosity per locus per sampling period were calculated using Microsatellite Toolkit in Excel (Park 2001). Differences in Rs, Pa, Ho and He among time periods were tested with the time period as a categorical variable using ANOVA in SYSTAT together with post hoc Tukey tests (version 8.0, 1998). Due to the small sample size for 1948–1949, effective population size (Ne) was estimated for the 1988–1991 and 2013–2015 time periods only using the linkage disequilibrium module based on a single point sample in NeEstimator (version 1.3; Peel et al. 2004; Waples and Do 2010). We conducted an AMOVA in GenAlEx for the 1988–1991 and 2013–2105 population wide datasets to assess the partitioning of genetic variation within and among habitat patches. We tested for evidence of genetic bottlenecks for patches with sufficient data (n ≥ 3 individuals) for both time periods using the two-phase (TPM), and single step (SMM) mutation models in BOTTLENECK (Cornuet and Luikart 1996).

A Principal Coordinates Analysis (PCoA) in GenAlEX, which uses a pairwise genetic distance matrix calculated from multi-allelic, multi-locus microsatellite data, was used to characterize any changes in the spatial partitioning of genetic diversity within BSHP across the three time periods and between BSHP and BIH. We then used STRUCTURE (version 2.3.4; Pritchard et al. 2000) to assess differences between genotype frequencies and therefore assignment probability to distinct genotype clusters (k) across the three time periods using a 50,000-iteration burn-in followed by ten 500,000 Markov Chain Monte Carlo (MCMC) replications per k, for k = 1–10 and the Δk method to determine the optimal number of genotype clusters (Evanno et al. 2005).

Results from an earlier genetic study, which used variable number of tandem repeat minisatellite nuclear DNA data (Peacock and Smith 1997), revealed frequent dispersal among patches together with high within patch mortality that resulted in a fluid movement dynamic among patches and rapidly changing population genetic structure across years. Therefore, we did not use an assignment test approach, such as Bayesian genotype clustering analysis, which is not well suited to capture such shifting dynamics, to assess population genetic structure across habitat patches within time periods. Instead, to further characterize genetic changes within and among patches between the two time periods, 1988–1990 and 2013–2015 respectively, we calculated population level relatedness (r) for each time period using the Queller and Goodnight method (1989) in GenAlEx, as well as r per patch for all patches with n > 3 individuals. We also calculated a genetic similarity metric among all individuals across the entire population, among individuals per patch and pairwise comparisons of genetic similarity among patches for each time period for which we had sufficient data (patch occupancy n ≥ 3 individuals) using the allele sharing matrix module in Microsatellite toolkit in Excel (Park 2001). We then conducted paired t-tests among patch genetic similarities for those patches with sufficient data from both 1988–1990 and 2013–2015 to test for statistically significant changes in genetic similarity within patches over time.

Patch occupancy

Long-term census survey data indicate a steady decline in total patch occupancy (Fig. 2). However, the slope is only significant when the early survey data from 1972 and 1977 are included (R² = 0.345, F_1,27 = 14.326, p = 0.001; Fig. 2a). Annual surveys since 1989 (excepting 1990, 2002 and 2007) reveal a highly fluctuating population size and although a decline in occupancy is also observed, it is non-significant (R² = 0.069, F_1,25 = 1.86, p = 0.185; Fig. 2b). This non-significant decline is likely due to the fact that by 1989, the bulk of the remaining population was found in the generally larger patches of the northern patch network. The decline in patch occupancy in the southern and middle patches was observed even for patches with sizes similar to those found in the northern portion of the patch network (i.e., same number of territories, Fig. 2c, d). Subpopulations in the southern patch network (all patches south of Silver Hill, see Fig. 1) were extirpated sometime between the 1977 and 1989 survey periods and although this part of the study area experienced a subsequent increase in occupancy in 2003 and 2004 (16% and 14.8%, respectively) occupancy declined to 3.7% in 2005 and was completely extirpated by 2008 (Fig. 3a). In 2015 and 2016, one and five individuals respectively were observed on the Red Cloud, a larger patch in the southern patch network with approximately 12 territories (Fig. 1). In addition, one individual was observed on the Booker Consolidated patch in 2016, also in the southern patch network (Fig. 1). The talus patches in the middle patch network (patches on Silver Hill south of Mono Shaft, see Fig. 1) were all mostly unoccupied by the early 2000s and have remained so (Fig. 3b). In contrast, the northern patch network (all patches north of and including Mono Shaft, see Fig. 1) has maintained consistently occupied patches over the same time period and although population size has fluctuated on these patches, the northern network has not experienced a decline in total occupancy either between 1972–2018 or 1989–2018, respectively (Fig. 3c).

Allelic dropout, null alleles and genotyping error rates

Null alleles were probable at multiple loci in 1988–1991 (OCP10 and 21) and 2013–2015 (OCP10 and 17), but there were no consistent patterns across loci, patch network, or time. The fecal DNA samples had an amplification success rate of 95.5% with average rates of allelic dropout and false alleles at 2.9% and 1.6%, respectively (Table 1). The museum specimens had a lower amplification success rate (90.2%), higher allelic dropout (3.9%) and more false alleles (4.4%) compared to fecal DNA samples from extant animals, which was expected considering age and degradation of samples. The ear tissue samples collected during 1988–1990 were not repeatedly amplified due to the high-quality nature of these samples.

Table 1

Locus-specific rates of allelic dropout (ε₁) and false alleles (ε ₂) of for museum and fecal DNA samples. Genotyping success rate per locus averaged across all samples (%).
	1949 Museum Sample DNA			2015 Fecal DNA
Locus	ε₁ (allelic dropout rate)	ε ₂ (false allele rate)	% success	ε₁ (allelic dropout rate)	ε ₂ (false allele rate)	% success
OCP 4b	0.031	0.138	0.923	0.000	0.010	0.930
OCP 6b	0.091	0.036	0.846	0.000	0.032	0.987
OCP 7	-	-	-	0.095	0.000	0.917
OCP 8	0.050	0.050	0.969	0.006	0.000	0.974
OCP 9b	0.000	0.015	0.892	0.006	0.000	0.939
OCP 10	0.031	0.046	0.954	0.023	0.000	0.996
OCP 17	0.015	0.000	0.877	0.052	0.070	0.956
OCP 21	0.033	0.050	0.877	0.014	0.033	0.934
Average	0.035	0.039	0.909	0.034	0.015	0.956

Genetic diversity

A total of 298 unique adult individuals were successfully genotyped from the contemporary and historic BSHP mining district as well as the historical samples from the Big Indian Mountain (BIH) site (tissue samples 1947–1949, n = 22; 1988–1990, n = 75; fecal samples 2013–2015, n = 191). The 1947–1949 museum samples failed to amplify at OCP7, which was subsequently removed from all analyses that included museum samples. The 1948–1949 dataset was too small to calculate heterozygous excess and deficiency on a per patch basis and we did not observe heterozygous excess or deficiency for any locus on any per patch for the 1988–1991 dataset. However, we observed both heterozygous excess and deficiency in the 2013–2015 fecal DNA dataset, but primarily heterozygous excess, for subpopulations on five of the 26 patches sampled, which is consistent with a genetic bottleneck signature (see genetic bottleneck results below) (Bull Run Χ² = 56.08, df = 14; High Peak Χ² = 52.86, df = 16; Hobart Χ² = 36.55, df = 16; Tioga Shaft Χ² = 36.83, df = 14; Vindicator Hoover Χ² = 49.77, df = 14; p < 0.000). Because fecal samples were collected fresh and individual collections were made from distinct latrine sites on individual territories, we conclude that we sampled single individuals per fecal sample collected. The probability of identity calculated for the 2013–2015 fecal DNA dataset after duplicate samples were removed using the Multilocus Matches DS module in GenAlEx was 1.9E-05.

We observed linkage disequilibrium at one locus pair (OCP10 x OCP21, p = 0.0001, 7280 permutations) in the 1988–1991 dataset for the High Peak subpopulation only. For the 2013–2015 dataset, three locus pairs were in linkage disequilibrium, but only for two patches, High Peak (OCP8 x OCP9) and Hobart (OCP4 x OCP8, OCP8 x OCP10) (p = 0.00009, 10,800 permutations). On a BSHP population wide basis, a significantly larger F_IS was observed for OCP10 in the 1948–1949 and 1988–1991 datasets (F_IS = 0.683, 0.501, p = 0.0018, 560 randomizations) and significantly smaller F_IS than random for OCP4 (F_IS=-0.286), OCP9 (F_IS=-0.236) and OCP 21 (F_IS = -0.232) (p = 0.0018, 560 randomizations) in the 2013–2015 dataset (Table S2). All genetic diversity metrics, number of alleles (Na), allelic richness (Rs), observed (Ho) and expected (He) heterozygosity and inbreeding coefficient (F_IS) are summarized per locus per time period in Table S2.

We observed significant declines in Rs (F = 10.228, p = 0.001), Pa (F = 20.89, p = 0.000), and He (F = 8.774, p = 0.006), but not Ho (F = 0.087, p = 0.917) over the three time periods at the BSHP site (1948–2015; Fig. 4a-d). Post hoc Tukey tests show that the 1948–1949 values for Rs and Pa were significantly higher than those for 1988–1991 and 2013–2015 (p ≤ 0.006), whereas the 1988–1991 and 2013–2015 values did not differ (p = 0.911, 0.732 respectively). The 1948–1949 He estimate was significantly higher than 2013 − 1015, but not 1988–1991 (p = 0.005, 0.107 respectively; Fig. 4c), whereas the 1988–1991 and 2013–2015 He values did not differ (p = 0.296).

Population genetic analyses

AMOVA results show 91% of the genetic variance was within individuals, 7% among individuals and 2% among patches in 1988–1991, with 91% within individuals and 9% among patches in 2013–2015 (999 permutations). The small sample size from 1948–1949 makes it difficult to accurately assess relatedness on a population wide basis for this early time period, but there was a clear trend in the increase in relatedness over time with low average relatedness in 1948–1949 (r = -0.068, 95% CI -0.113-0.128) and 1988–1991 (r = 0.024, 95% CI -0.060-0.060) and an r value consistent with an average relatedness among first cousins (r = 0.158, 9% CI -0.028-0.021) in 2013–2015 (Fig. 4f). Consistent with the relatedness results, pairwise genetic similarity among all BSHP individuals per time period show a significant increase in population level genetic similarity across time (Kruskal-Wallis Test Statistic = 1780.085, df = 2, p = 0.000; Fig. 5a). For the subset of patches for which we had sufficient data from both 1988–1991 and 2013–2015, comparisons among individuals also showed a significant increase in genetic similarity within patches (Mann-Whitney U, 10.5–97595.5, df = 1, p = 0.000-0.009; Fig. 5b), as well as increases in average relatedness for all paired patch comparisons with the exception of High Peak, which had an average r of 0.01 and − 0.003 for the 1988–1991 and 2013–2015 respectively (Fig. S1). We note that High Peak is centrally located within the northern patch network and also one of the largest patches present in the study area.

We did not observe significant pairwise F_ST estimates among patch subpopulations for the 1988–1991 dataset, but did observe multiple significant F_ST estimates for patch pairs during the 2013–2015 time period (n = 12 significant F_ST estimates, p = 0.0006, 1560 permutations; Table S3). The majority of the significant pairwise F_ST estimates were between core and peripheral populations in the northern patch network (Table S3). Global estimates of θ for the latter two time periods, with extensive sampling across all occupied patches, show an increase in population genetic structure (1988–1991, θ = 0.021, 95% CI 0.012–0.03; 2013–2015, θ = 0.087, 95% CI 0.067–0.107). Pairwise F_ST between BSHP and BIH were significant across all time periods, with increasing pairwise F_ST over time (BIH 1947 and BSHP 1948–1949 F_ST = 0.152, BIH 1947 and BSHP 1988–1991 F_ST = 0. 246, BIH 1947 and BSHP 2013–2015 F_ST = 0.361; p = 0.008)

STRUCTURE results for samples from all three time periods including the historical BIH samples indicate that k = 3 was the best fit of the data (Fig. 6a, b). Assignment to one of the genotype clusters was composed of 1947–1949 individuals from BSHP and BIH only. No individuals from BIH and BSHP historic samples (1947–1949) assigned to the predominant genotype cluster identified for individuals sampled in 2013–2015 and only limited membership was observed in this cluster for individuals from 1988–1991, indicating allelic loss and changes in allele frequency, commensurate with observed declines in allelic richness and private alleles over this same time period (Fig. 5c). These results were also supported by the PCoA analysis, which separated the 2013–2015 dataset on coordinate axis 1 from the 1947–1949 and 1988–1991 datasets (Fig. 5d), whereas the BSHP 1948–1949 and 1988–1991 samples largely overlapped in PCoA space and the BIH 1947 samples overlapped completely with BSHP 1948–1949 samples.

Effective population size declined between the two later time periods 1988–1991 (Ne = 76, 95% CI 47.4-149.8) and 2013–2015 (Ne = 47.1, 95% CI 38.4–58.3; Fig. 4e). Genetic bottlenecks were observed for multiple patch subpopulations in 1988–1991 and all but one patch showed evidence of genetic bottlenecks in the 2013–2015 dataset (Table 2).

Table 2

Bottleneck results per patch for 1988–1991 and 2013–2015, calculated for individual patch populations with n ≥ 3 individuals. Significant p values are bolded and italicized.
	1988–1991	p values		2013–2015	p values
Patch Identification	N	TPM	SMM	N	TPM	SMM
BECHTEL (BEC)	-	-	-	9	0.008	0.009
BLACK HAWK (BHS)	-	-	-	10	0.011	0.015
BLUE POINT (BPS)	-	-	-	8	0.011	0.011
BODIE55 (BD55)	-	-	-	5	0.010	0.011
BODIE56 (BD56)	-	-	-	7	0.016	0.018
BODIE68 (BD68)	6	0.075	0.272	4	0.031	0.043
BODIE OLD SHAFT (BOS)	-	-	-	6	0.012	0.016
BULL RUN (BULL)	12	0.072	0.220	23	0.008	0.008
HIGH PEAK70-72 (HP70-72)	-	-	-	9	0.102	0.103
HIGH PEAK (HP)	23	0.010	0.012	41	0.014	0.104
HOBART (HOB)	-	-	-	12	0.007	0.008
LUCKY JACK (LJ)	3	0.013	0.015	6	0.009	0.010
SCANDINAVIAN (SCAN)	3	0.079	0.082	-	-	-
SECURITY TUNNEL (SEC)	-	-	-	8	0.025	0.012
SUMMIT (SUM)	13	0.008	0.010	4	0.034	0.043
TIOGA SHAFT (TGA)	3	0.013	0.013	15	0.011	0.011
VINDICATOR HOOVER (VHO)	-	-	-	12	0.015	0.015

The results of this study show that the BSHP population has experienced a steady decline in patch occupancy since the mid-20th century along with a corresponding loss of genetic diversity (as measured by Rs, Pa and He), a decrease in effective population size and an increase in population genetic structure. Both population wide and individual patch levels of relatedness and genetic similarity have also increased. We observed not only significant genetic differentiation among patch subpopulations as of the last time period sampled (2013–2015), but also significant genetic differentiation among the entire patch network population across all three time periods. Although the northern half of the study area has remained consistently occupied since at least the mid-20th century (Severaid 1955), the decline in overall effective population size suggests the pika metapopulation at BSHP may be undergoing genetic coalescence as described by Gilpin (1991) despite the continuously occupied patches in the northern half of the study area.

The range-wide distribution and structure of genetic variation enables species to tolerate or adapt to changing conditions and is often highly dependent upon movement patterns of individuals within and among local populations (Booy et al. 2000; Bálint et al. 2011; Pauls et al. 2013). The pika is a small alpine mammal with limited dispersal capabilities (Smith and Ivins 1983; Peacock 1997) and a very narrow range of thermal tolerance (MacArther and Wang 1973). Patterns observed in this study suggest that the dispersal potential of pikas may be changing at BSHP. The movement dynamics among patches in the earlier Peacock and Smith (1997) study, which used both mark-recapture and genetic data to document movement patterns and population genetic structure, suggested a fluid movement dynamic across the northern patch network. Although juvenile pikas tended to settle on the first available territory in the Peacock and Smith (1997) study, territory availability varied over space and time such that levels of genetic similarity within and among patch subpopulations was transitory and shifted among years. Both higher and lower levels of genetic similarity among individuals than would be expected by chance were observed on both single patches and for neighboring groups of patches, which tended to be the more peripheral subpopulations (Peacock and Smith 1997). In 2013–2015, we observe statistically significant differentiation among not only groups of neighboring peripheral subpopulations but also between large core patches such as High Peak and Tioga Shaft with more peripherally located subpopulations, a pattern that was not observed in the late 20th century (1988–1991).

Despite past evidence of an active metapopulation dynamic (Peacock and Smith 1997; Moilanen et al. 1998) and significant inter-annual flux in both patch occupancy and dispersion of genetic variation among patches, patch population dynamics at BSHP appear to be changing. Genetic variation within the habitat patch network has a larger among-patch partition now than it did approximately 30 years ago (Peacock and Smith 1997). These results suggest that occupied patches at BSHP are becoming more isolated, perhaps due to lack of either initiated or successful dispersal among patches. Under this scenario of reduction in dispersal and connectivity, losses of genetic variation will accelerate over time via random genetic drift – a result of demographic and genetic stochasticity inherent to small populations (Lande 1993; Reed and Frankham 2003; Whiteley et al. 2015). However, despite significant declines in allelic richness and effective population size, the BSHP pika metapopulation has retained similar levels of observed heterozygosity (Ho) since the mid-20th century, although our limited sample size from 1949 likely masks the true variation that once existed. This result, in addition to broad distributions for genetic similarity on each patch, suggests ongoing juvenile dispersal and movement of alleles among extant patches via an active albeit altered metapopulation dynamic, which may in fact act to somewhat decrease the extinction probability of the entire patch network. It is possible that territorial turnover on natal patches may have become more frequent, leading to increased probability of recruitment to the natal patch for juvenile animals and the likelihood of increased genetic similarity and relatedness within patches as we observed. Although individuals avoid an increased mortality risk associated with between-patch dispersal by recruiting to their natal patch, decreased rates of among-patch dispersal limits the opportunity for the spread of novel genetic variants and increases loss of genetic variation through random genetic drift. Increased rates of patch turnover, especially on smaller patches, may turn these patches into sink habitats and lead to genetic coalescence of the genetic variation found in the larger, more stable patches. Metapopulation theory predicts accelerated population decline and coalescence of genetic variation as the probability of re-colonization of extirpated patches and dispersal among extant patches decreases (Gilpin 1991). Indeed, the among-patch comparisons in relatedness and genetic similarity show a shift to higher mean values in 2013–2015 and increases in overall population wide r and genetic similarity when compared with patch dynamics from 1988–1991.

The BSHP population is declining with little to no recolonization potential from natural talus patches within the Bodie Hills, as these patches are largely extirpated. Even if still occupied by pikas, these natural talus patches are themselves isolated in expanses of sage scrub habitat (Nichols 2011). It is unlikely that the BSHP pika metapopulation has experienced any gene flow with other pika populations in the Bodie Hills since the mid-20th century increasing both the regional isolation of the BSHP population and the impact of random genetic drift. This hypothesis of isolation is further supported by the changing pattern of genetic overlap between the mid-century BIH samples and the three time periods sampled at BSHP. Although we observed complete overlap between BIH with mid-century BSHP samples, by 2015 this gene pool has drifted and appears distinct from earlier time periods.

Decreases in mean annual precipitation and concomitant increases in minimum January temperature over the last 65 years may be indicative of a long-term warming trend within this region and the importance of mean annual precipitation given the severe drought that recently occurred in California during 2012–2014 (Griffin and Anchukaitis 2014; Diffenbaugh et al. 2015). Although not directly tested in this study, several lines of evidence suggest climate change may have played a role in the decline of the BSHP population. The middle and southern portions of the population are dominated by small habitat patches with very few medium or larger sized patches available. The loss of subpopulations on patches in the middle and southern portions of the study site suggests those patches no longer provide either suitable microclimate and/or sufficient forage habitat for pikas, thereby increasing the probability of overwinter mortality, resulting in the collapse of the metapopulation dynamic in this half of the study area. In the near term, with projections suggesting increasing ambient temperatures and changes in precipitation regimes for this region (Cayan et al. 2008; Jepsen et al. 2016), the smaller ore dumps in the northern half of BSHP are also likely to lose their thermal buffering capacity, thus no longer providing adequate refugia for pikas and becoming unsuitable sink habitat. This, in turn, will alter the spatial structure and proximity of the remaining talus patches that can be occupied. Indeed, the observed increase in genetic similarity and relatedness within extant patches in the northern half of the study area suggests that some underlying mechanisms related to habitat suitability and a changing climate may already be at work.

Population resilience and the role of neutral genetic variation in species’ persistence

Population level responses to local environmental conditions will ultimately determine the ability of species to persist across their ranges. Populations on the edge of their distribution are often considered, at least theoretically, to be less capable of mounting an evolutionary response largely due to limited genetic variation, stressful environments, and smaller population sizes (Rogers and Peacock 2012; Bouchard et al. 2017). Indeed, previous research investigating population genetic structure in low elevation sites in the Columbia River Gorge found that pika populations at the edge of their bioclimatic envelope were characterized by relatively low allelic richness and heterozygosity with very little evidence of active gene flow (Robson et al. 2016). However, despite being isolated, the BSHP population has maintained unique allelic variants into the 21st century, suggesting potentially adaptive differences from other populations in the region (Klingler, Jahner et al. 2021). Some differences have already been noted including vocalizations (Conner 1982; Trefry and Hik 2010) and litter sizes, which tend to be larger on average in BSHP than those in the Sierra Nevada (3.68 and 3.12 respectively; Smith 1978). Levels of intraspecific and even subspecific lineage diversity may differ significantly between pika populations in the range core compared to sites that occur on the trailing edge or range peripheries (Lesica and Allendorf 1995; Robson et al. 2016; Waterhouse et al. 2017; Klingler, Jahner et al. 2021).

Despite the fact that multiple studies have shown that persistence on the edge can be fairly common, especially when extirpation across more central parts of the range is occurring, (Gibson et al. 2009), the question becomes as climate change accelerates and genetic variation erodes, do isolated populations like BSHP have enough genetic resources to adapt in situ? The rapid loss of genetic variation due to random genetic drift associated with increasing isolation and climatic correlates does not bode well for an adaptation scenario and very likely could result in extinction of unique lineages (Robson et al. 2016; David and Wright 2017).

Although pika populations across their range have experienced repeated fragmentation and isolation during previous periods of climate change, the current rates of warming may preclude adaptive responses (Grayson 2005; Galbreath et al. 2009; 2010). While regional cohesion has been maintained due to admixture during interglacial periods, Holocene era fragmentation has generated significant genetic differentiation among populations and regions with a potential increase in extirpation risk (Galbreath et al. 2009; 2010; Robson et al. 2016; Klingler, Jahner et al. 2021). Indeed, over the past 25–30 years, coincident with warming temperatures, the incidence of extirpations across the American pika’s range, primarily in the Great Basin and within the Bodie Hills regions, has been particularly pronounced (Beever et al. 2003; Stewart et al. 2015; Nichols et al. 2016).

Projections suggest that in the Sierra Nevada and Great Basin, temperatures will continue to warm significantly especially during the summer season even under low-emission scenarios (+ 2°C by 2100; IPCC, 2014) with significant declines in snowpack (Cayan et al. 2008; Jepsen et al. 2016). The sudden extirpation of the New York Hill pika population (Nichols et al. 2016), located approximately 20 km from BSHP in sagebrush-scrub and in similar man-made subdivided ore dump habitat, suggests a loss of the metapopulation dynamic and an ensuing extinction vortex. The potential vulnerability of these populations due to unpredictable or stochastic factors such as greater exposure to the expanding distributions of novel disease vectors (e.g. plague and lagomorph hemorrhagic fever), potential predators (e.g. western rattlesnake, Crotalus viridis), extreme weather events, and competitors for forage and habitat (e.g. Bushy-tailed woodrat, Neotoma cinerea or Yellow-bellied marmot (Marmota flavivientris) (Lowell et al. 2015; Benson et al. 2019; Ray et al. 2016; Foley et al. 2017) should be considered in discussions concerning pika persistence. The continued loss of genetic diversity through population extirpation or long-term genetic erosion in small, isolated populations may represent a significant evolutionary cost at the species level by reducing its capacity to respond to ongoing environmental change. This study is the first to characterize temporally stratified patterns of population genetic variation and structure in an iconic American pika metapopulation. Although our study indicates a loss of neutral genetic variation across time for one pika metapopulation, future research would benefit from analyses of adaptive trait loci for multiple populations across the full range of the American pika (Waterhouse et al. 2018; Schmidt et al. 2021).

Online Resources The online version contains supplementary material available at http://XXX

Acknowledgements We thank Chris Conroy (UC Berkeley-MVZ) for approval of specimen loans, access to specimens and assistance with sample preparation and shipping. We thank UNR Nevada Genomics Center (http://www.ag.unr.edu/genomics/) for DNA sequencing and fragment analysis. Special thanks to BSHP staff and volunteers for continued logistical support especially California State Park Supervising Rangers, J. Heitzmann, T. Peters, and R. Peek. Field work would not have been possible without K. King and P. Nichols. Thank you to V. Kirchoff and L. Mayfield for laboratory support and assistance.

Author contribution KBK and MMP designed the study. KBK, MMP, and ERH collected pika samples from the field and accessed MVZ specimens. KBK and ERH conducted the laboratory analyses. KBK and MMP analyzed the data. LBN provided all census data. KBK and MMP wrote the manuscript. All authors read and approved the final manuscript.

Funding Funds for travel, sample collection, laboratory analysis and genotyping was provided by a Bodie Foundation grant (Bodie Foundation is a 501 (c) (3) non-profit corporation), the University of Nevada (UNR) Ecology, Evolution and Conservation Biology graduate program summer stipend, and the UNR Diana Hadley-Lynch Scholarship awards to KBK.

Data availability Microsatellite genotypes for this paper are available at Dryad, Dataset, https://doi.org/10.5061/dryad.51c59zwbd http://www

Code availability Not applicable

Conflict of Interest The authors declare that they have no conflicts of interest.

Ethical approval and consent to participate Institutional Animal Care and Use Committees (IACUC) at the University of Nevada, Reno (protocol No. 00557). Scientific collecting permits were obtained from the California Department of Parks and Recreation (No. 2014-

03); California Department of Fish and Wildlife (SC-012184); US Department of

Agriculture (USDA) Forest Service and USDA Forest Service, Inyo National

Forest (No. LVD14021 and No. RMT143)

Consent for publication Not applicable

AghaKouchak A, Cheng L, Mazdiyasni O, Farahmand A (2014) Global warming and changes in risk of concurrent climate extremes: insights from the 2014 California drought. Geophys Res Lett 41:8847–8852. https://doi.org/10.1002/2014GL062308
Anderson TM, White S, Davis B, Erhardt R, Palmer M, Swanson A, Kosmala M, Packer C (2016) The spatial distribution of African savannah herbivores: species associations and habitat occupancy in a landscape context. Phil Trans R Soc B 371:20150314. https://doi.org/10.1098/rstb.2015.0314
Bálint M, Domisch S, Engelhardt CHM, Haase P, Lehrian S, Sauer J, Theissinger K, Pauls SU, Nowak C (2011) Cryptic biodiversity loss linked to global climate change. Nat Clim Chang 1:313–318. https://doi.org/10.1038/nclimate1191
Beever EA, Wilkening JL, McIvor DE, Weber SS, Brussard PF (2008) American pikas (Ochotona princeps) in northwestern Nevada: a newly discovered population at a low-elevation site. West North Am Nat 68:8–14. https://www.jstor.org/stable/41717651
Beever EA, Ray C, Mote PW, Wilkening JL (2010) Testing alternative models of climate-mediated extirpations. Ecol Appl 20:164–178. https://doi.org/10.1890/08-1011.1
Beever E, Ray C, Wilkening J (2011) Contemporary climate change alters the pace and drivers of extinction. Glob Chang 17:2054–2070. https://doi.org/10.1111/j.1365-2486.2010.02389.x
Beever EA, Perrine JD, Rickman T, Flores M, Clark JP, Waters C, Weber SS, Yardley B, Thoma D, Chesley-Preston T, Goehring KE, Magnuson M, Nordensten N, Nelson M, Collins GH (2016) Pika (Ochotona princeps) losses from two isolated regions reflect temperature and water balance, but reflect habitat area in a mainland region. J Mammal 97:1495–1511. https://doi.org/10.1093/jmammal/gyw128
Benson JF, Mahoney PJ, Vickers TW, Sikich JA, Beier P, Riley SPD, Ernest HB, Boyce WM (2019) Extinction vortex dynamics of top predators isolated by urbanization. Ecol Appl 29:e01868. https://doi.org/10.1002/eap.1868
Blomqvist D, Pauliny A, Larsson M, Flodin L-A (2010) Trapped in the extinction vortex? Strong genetic effects in a declining vertebrate population. BMC Evol Biol 10:33. https://doi.org/10.1186/1471-2148-10-33
Booy G, Hendriks RJ, Smulders MJM, Van Groenendael JM, Vosman B (2000) Genetic biodiversity and the survival of populations. Plant Bio 2:379–395. http://doi.org/10.1055/s-2000-5958
Bouchard JR, Fernando DD, Bailey SW, Weber-Townsend J, Leopold DJ (2017) Contrasting patterns of genetic variation in central and peripheral populations of Dryopteris fragrans (Fragrant wood fern) and implications for colonization dynamics and conservation. Int J Plant Sci 178:607–617. http://doi.org/10.1086/693109607
Calkins MT, Beever EA, Boykin KG, Frey JK, Andersen MC (2012) Not-so-splendid isolation: modeling climate-mediated range collapse of a montane mammal Ochotona princeps across numerous ecoregions. Ecography 35:780–791. https://doi.org/10.1111/j.1600-0587.2011.07227.x
Castillo JA, Epps CW, Davis AR, Cushman SA (2014) Landscape effects on gene flow for a climate-sensitive montane species, the American pika. Mol Ecol 23:843–856. https://doi.org/10.1111/mec.12650
Cayan DR, Maurer EP, Dettinger MD, Tyree M, Hayhoe K (2008) Climate change scenarios for the California region. Clim Change 87:S21–S42. https://doi.org/10.1007/s10584-007-9377-6
Chapuis MP, Estoup A (2007) Microsatellite null alleles and estimation of population differentiation. Mol Biol Evol 24:621–631. https://doi.org/10.1093/molbev/msl191
Conner DA (1982) Geographic-variation in short calls of pikas (Ochotona princeps). J Mammal 63:48–52. https://doi.org/10.2307/1380670
Cornuet JM, Luikart G (1996) Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics 144(4):2001–2014. https://doi.org/10.1093/genetics/144.4.2001
Creel S, Spong G, Sands JL, Rotella J, Zeigle J, Joe L, Murphy KM, Smith D (2003) Population size estimation in Yellowstone wolves with error-prone noninvasive microsatellite genotypes. Mol Ecol 12:2003–2009. https://doi.org/10.1046/j.1365-294X.2003.01868.x
David SR, Wright JJ (2017) Genetic variation and biogeography of the spotted gar Lepisosteus oculatus from core and peripheral populations. J Exp Zool Part B Mol Dev Evol 328:596–606. https://doi.org/10.1002/jez.b.22772
Dearing M (1997) The function of haypiles of pikas (Ochotona princeps). J Mammal 78:1156–1163. https://doi.org/10.2307/1383058
Diaz HF, Wahl ER (2015) Recent California water year precipitation deficits: A 440-year perspective. J Clim 28:4637–4652. https://doi.org/10.1175/JCLI-D-14-00774.1
Diffenbaugh NS, Swain DL, Touma D (2015) Anthropogenic warming has increased drought risk in California. Proc Natl Acad Sci USA 112:3931–3936. https://doi.org/10.1073/pnas.1422385112
Erb LP, Ray C, Guralnick R (2011) On the generality of a climate-mediated shift in the distribution of the American pika (Ochotona princeps). Ecology 92:1730–1735. https://doi.org/10.1890/11-0175.1
Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE, a simulation study. Mol Ecol 14:2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x
Foley P, Roth T, Foley J, Ray C (2017) Rodent–pika parasite spillover in western North America. J Med Entomol 54:1251–1257. https://doi.org/10.1093/jme/tjx085
Galbreath KE, Hafner DJ, Zamudio KR (2009) When cold is better: climate-driven elevation shifts yield complex patterns of diversification and demography in an alpine specialist (American pika, Ochotona princeps). Evolution 63:2848–2863. https://doi.org/10.1111/j.1558-5646.2009.00803.x
Galbreath KE, Hafner DJ, Zamudio KR, Agnew K (2010) Isolation and introgression in the intermountain west: contrasting gene genealogies reveal the complex biogeographic history of the American pika (Ochotona princeps). J Biogeogr 37:344–362. https://doi.org/10.1111/j.1365-2699.2009.02201.x
Gibson SY, Van Der Marel RC, Starzomski BM (2009) Climate change and conservation of leading-edge peripheral populations. Conserv Biol 23:1369–1373. https://doi.org/10.1111/j.1523-1739.2009.01375.x
Gilpin M (1991) The genetic effective size of a metapopulation. Biol J Linn Soc 42(1–2):165–175. https://doi.org/10.1111/j.1095-8312.1991.tb00558.x
Goudet J, Perrin N, Waser P (2002) Tests for sex-biased dispersal using bi-parentally inherited genetic markers. Mol Ecol 11:1103–1114. https://doi.org/10.1046/j.1365-294X.2002.01496.x
Grayson DK (2005) A brief history of Great Basin pikas. J Biogeogr 32:2103–2111. https://doi.org/10.1111/j.1365-2699.2005.01341.x
Griffin D, Anchukaitis KJ (2014) How unusual is the 2012–2014 California drought? Geophys Res Lett 41:9017–9023. https://doi.org/10.1002/2014GL062433
Hafner DJ (1994) Pikas and permafrost: post-Wisconsin historical zoogeography of Ochotona in the southern Rocky Mountains, USA. Arct Alp Res 26:375–382. https://doi.org/10.1080/00040851.1994.12003082
Hafner DJ, Smith AT (2010) Revision of the subspecies of the American pika, Ochotona princeps (Lagomorpha: Ochotonidae). J Mammal 91:401–417. https://doi.org/10.1644/09-MAMM-A-277.1
Hekkala E, Shirley MH, Amato G, Austin JD, Charter S, Thorbjarnarson J, Vliet KA, Houck ML, DeSalle R, Blum MJ (2011) An ancient icon reveals new mysteries: mummy DNA resurrects a cryptic species within the Nile crocodile. Mol Ecol 20:4199–4215. https://doi.org/10.1111/j.1365-294X.2011.05245.x
Jeffress MR, Rodhouse TJ, Ray C, Wolff S, Epps CW (2013) The idiosyncrasies of place: geographic variation in the climate-distribution relationships of the American pika. Ecol Appl 23:864–878. https://doi.org/10.1890/12-0979.1
Jeffress MR, Van Gunst KJ, Millar CI (2017) A surprising discovery of American pika sites in the northwestern Great Basin. West N Am Naturalist 77(2):252–268. https://doi.org/10.3398/064.077.0213
Jepsen SM, Harmon TC, Meadows MW, Hunsaker CT (2016) Hydrogeologic influence on changes in snowmelt runoff with climate warming: numerical experiments on a mid-elevation catchment in the Sierra Nevada, USA. J Hydrol 533:332–342. http://dx.doi.org/10.1016/j.jhydrol.2015.12.010
Johnson PCD, Haydon DT (2007) Maximum-likelihood estimation of allelic dropout and false allele error rates from microsatellite genotypes in the absence of reference data. Genetics 175:827–842. https://doi.org/10.1534/genetics.106.064618
Kalinowski ST (2005) HP-RARE 1.0: A computer program for performing rarefaction on measures of allelic richness. Mol Ecol Notes 5:187–189. https://doi.org/10.1111/j.1471-8286.2004.00845.x
Klingler KB, Jahner JP, Parchman TL, Ray C, Peacock MM (2021) Genomic variation in the American pika: signatures of geographic isolation and implications for conservation. BMC Ecol Evo 21:2. https://doi.org/10.1186/s12862-020-01739-9
Lande R (1993) Risks of population extinction from demographic and environmental stochasticity and random catastrophes. Am Nat 142:911–927. https://doi.org/10.1086/285580
Lesica P, Allendorf FW (1995) When are peripheral populations valuable for conservation? Conserv Biol 9:753–760. https://doi.org/10.1046/j.1523-1739.1995.09040753.x
Lowell JL, Antolin MF, Andersen GL, Hu P, Stokowski RP, Gage KL (2015) Single-nucleotide polymorphisms reveal spatial diversity among clones of Yersinia pestis during plague outbreaks in Colorado and the western United States. Vector-Borne Zoonotic Dis 15:291–302. https://doi.org/10.1089/vbz.2014.1714
MacArthur R, Wang L (1973) Physiology of thermoregulation in the pika, Ochotona princeps. Can J Zool 51:11–16. https://doi.org/10.1139/z73-002
Manning T, Hagar JC (2011) Use of nonalpine anthropogenic habitats by American pikas (Ochotona princeps) in western Oregon. West N Am Naturalist 71:106–112. https://doi.org/10.3398/064.071.0114
Mathewson PD, Moyer-Horner L, Beever EA, Briscoe NJ, Kearney M, Yahn JM, Porter WP (2017) Mechanistic variables can enhance predictive models of endotherm distributions: the American pika under current, past, and future climates. Glob Chang Biol 23:1048–1064. https://doi.org/10.1111/gcb.13454
McKelvey K, Schwartz M (2005) Dropout: a program to identify problem loci and samples for noninvasive genetic samples in a capture-mark‐recapture framework. Mol Ecol Notes 5:716–718. https://doi.org/10.1111/j.1471-8286.2005.01038.x
Moilanen A, Smith A, Hanski I (1998) Long-term dynamics in a metapopulation of the American pika. Am Nat 152:530–542. https://doi.org/10.1086/286188
Moyer-Horner L, Mathewson PD, Jones GM, Kearney MR, Porter WP (2015) Modeling behavioral thermoregulation in a climate change sentinel. Ecol Evol 5:5810–5822. https://doi.org/10.1002/ece3.1848
Neam KD, Lacher TE Jr (2018) Multi-scale effects of habitat structure and landscape context on a vertebrate with limited dispersal ability (the brown-throated sloth, Bradypus variegatus). Biotropica 50:684–693. https://doi.org/10.1111/btp.12540
Nichols L (2011) Fecal pellets of American pikas (Ochotona princeps) provide a crude chronometer for dating patch occupancy. West North Am Nat 70:500–507. https://doi.org/10.3398/064.070.0410
Nichols LB, Klingler KB, Peacock MM (2016) American pikas (Ochotona princeps) extirpated from the historic Masonic mining district of eastern California. West North Am Nat 76:163–171. https://doi.org/10.3398/064.076.0203
Park S (2001) Microsatellite toolkit for Excel, v.3.1. Trinity College. University of Dublin, Dublin
Pauls SU, Nowak C, Bálint M, Pfenninger M (2013) The impact of global climate change on genetic diversity within populations and species. Mol Ecol 22:925–946. https://doi.org/10.1111/mec.12152
Peacock MM (1997) Determining natal dispersal patterns in a population of North American pikas (Ochotona princeps) using direct mark-resight and indirect genetic methods. Behav Ecol 8:340–350. https://doi.org/10.1093/beheco/8.3.340
Peacock MM, Smith AT (1997) The effect of habitat fragmentation on dispersal patterns, mating behavior, and genetic variation in a pika (Ochotona princeps) metapopulation. Oecologia 112:524–533. https://doi.org/10.1007/s004420050341
Peacock MM, Kirchoff VS, Merideth SJ (2002) Identification and characterization of nine polymorphic microsatellite loci in the North American pika, Ochotona princeps. Mol Ecol Notes 2:360–362. https://doi.org/10.1046/j.1471-8286.2002.00249.x
Peacock MM, Kirchoff VS (2009) Identification and characterization of 19 polymorphic microsatellite loci in the North American pika, Ochotona princeps, Unpublished Genbank Accession numbers OCP10–28(GQ461705.1-GQ461722.1)
Peacock MM, Hekkala ER, Kirchoff VS, Heki LG (2017) Return of a giant: DNA from archival museum samples helps to identify a unique cutthroat trout lineage formerly thought to be extinct. R Soc Open Sci 4:171253. https://doi.org/10.1098/rsos.171253
Peacock M (2022) Life on the edge - a changing genetic landscape within an iconic American pika metapopulation over the last half century. Dryad, Dataset. https://doi.org/10.5061/dryad.51c59zwbd
Peakall R, Smouse PE (2006) GENALEX 6: Genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes 6:288–295. https://doi.org/10.1111/j.1471-8286.2005.01155.x
Peakall R, Smouse PE (2012) GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research–an update. Bioinformatics 28:2537–2539. https://doi.org/10.1111/j.1471-8286.2005.01155.x
Peel D, Ovenden JR, Peel SL (2004) NeEstimator: Software for estimating effective population size, version 1. 3. Department of Primary Industries and Fisheries, Queensland Government. City East, Queensland, Australia
Pinto SM, MacDougall AS (2010) Dispersal limitation and environmental structure interact to restrict the occupation of optimal habitat. Am Nat 175:675–686. https://doi.org/10.1086/652467
Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155:945–959. https://doi.org/10.1534/genetics.116.195164
Queller DC, Goodnight KF (1989) Estimating relatedness using genetic markers. Evolution 42:258–275. https://doi.org/10.1111/j.1558-5646.1989.tb04226.x
Ray C, Beever EA, Rodhouse TJ (2016) Distribution of a climate-sensitive species at an interior range margin. Ecosphere 7:e01379. https://doi.org/10.1002/ecs2.1379
Reed DH, Frankham R (2003) Correlation between fitness and genetic diversity. Conserv Biol 17:230–237. https://doi.org/10.1046/j.1523-1739.2003.01236.x
Robson KM, Lamb CT, Russello MA (2016) Low genetic diversity, restricted dispersal, and elevation-specific patterns of population decline in American pikas in an atypical environment. J Mammal 97:464–472. https://doi.org/10.1093/jmammal/gyv191
Rodhouse TJ, Beever EA, Garrett LK, Irvine KM, Munts M, Ray C, Shardlow MR (2010) Distribution of the Lava Beds pika (Ochotona princeps Goldmani): conservation implications from the range periphery. J Mammal 91:1287–1299. https://doi.org/10.1644/09-MAMM-A-334.1
Rogers SD, Peacock MM (2012) The disappearing northern leopard frog (Lithobates pipiens): conservation genetics and implications for remnant populations in western Nevada. Ecol Evol 2:2040–2056. https://doi.org/10.1002/ece3.308
Rousset F (2008) Genepop’007: a complete re-implementation of the genepop software for Windows and Linux. Mol Ecol 8:103–106. https://doi.org/10.1111/j.1471-8286.2007.01931.x
Schwalm D, Epps CW, Rodhouse TJ, Monahan WB, Castillo JA, Ray C, Jeffress MR (2016) Habitat availability and gene flow influence diverging local population trajectories under scenarios of climate change: a place-based approach. Glob Chang Biol 22:1572–1584. https://doi.org/10.1111/gcb.13189
Schmidt DA, Waterhouse MD, Sjodin BMF, Russello MA (2021) Genome-wide analysis reveals associations between climate and regional patterns of adaptive divergence and dispersal in American pikas. Heredity 127:443–454. https://doi.org/10.1038/s41437-021-00472-3
Severaid JH (1955) The Natural History of the Pikas (mammalian Genus Ochotona) h (Doctoral dissertation, University of California, Berkeley)
Shen G, He F, Waagepetersen R, Sun IF, Hao Z, Chen ZS, Yu M (2013) Quantifying effects of habitat heterogeneity and other clustering processes on spatial distributions of tree species. Ecology 94:2436–2443. https://doi.org/10.1890/12-1983.1
Smith AT (1974a) The distribution and dispersal of pikas: consequences of insular population structure. Ecology 55:1112–1119. https://doi.org/10.2307/1940361
Smith AT (1974b) The distribution and dispersal of pikas: influences of behavior and climate. Ecology 55:1368–1376. https://doi.org/10.2307/1935464
Smith AT (1978) Comparative demography of pikas (Ochotona) – effect of spatial and temporal age-specific mortality. Ecology 59:133–139. https://doi.org/10.2307/1936639
Smith AT, Ivins BL (1983) Colonization in a pika population: dispersal vs philopatry. Behav Ecol Sociobiol 13:37–47. https://doi.org/10.1007/BF00295074
Smith AT, Weston ML (1990) Ochotona princeps. Mammal Species 352:1–8. https://doi.org/10.2307/3504319
Stewart JAE, Perrine JD, Nichols LB, Thorne JH, Millar CI, Goehring KE, Massing CP, Wright DH (2015) Revisiting the past to foretell the future: summer temperature and habitat area predict pika extirpations in California. J Biogeogr 42:880–890. https://doi.org/10.1111/jbi.12466
Stewart JAE, Wright DH, Heckman KA (2017) Apparent climate-mediated loss and fragmentation of core habitat of the American pika in the northern Sierra Nevada, California, USA. PLoS ONE 12:e0181834. https://doi.org/10.1371/journal.pone.0181834
Stewart JAE, Wright DH (2012) Assessing persistence of the American pika at historic localities in California's Northern Sierra Nevada. Wildl Soc Bull 36:759–764. https://doi.org/10.1002/wsb.220
Trefry SA, Hik DS (2010) Variation in pika (Ochotona collaris, O. princeps) vocalizations within and between populations. Ecography 33:784–795. https://doi.org/10.1111/j.1600-0587.2009.05589.x
Van Oosterhout C, Hutchinson W (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes 4:535–538. https://doi.org/10.1111/j.1471-8286.2004.00684.x
Varner JM, Dearing D (2014) Dietary plasticity in pikas as a strategy for atypical resource landscapes. J Mammal 95:72–81. https://doi.org/10.1644/13-MAMM-A-099.1
Waples RS, Do C (2010) Linkage disequilibrium estimates of contemporary Ne using highly variable genetic markers: a largely untapped resource for applied conservation and evolution. Ecol Appl 3:244–262. https://doi.org/10.1111/j.1752-4571.2009.00104.x
Waterhouse MD, Blair C, Larsen KW, Russello MA (2017) Genetic variation and fine-scale population structure in American pikas across a human-modified landscape. Conserv Genet 18:825–835. https://doi.org/10.1007/s10592-017-0930-1
Wauters LA, Verbeylen G, Preatoni D, Martinoli A, Matthysen E (2010) Dispersal and habitat cuing of Eurasian red squirrels in fragmented habitats. Popul Ecol 52:527–536. https://doi.org/10.1007/s10144-010-0203-z
Whiteley AR, Fitzpatrick SW, Funk WC, Tallmon DA (2015) Genetic rescue to the rescue. Trends Ecol Evol 30:42–49. https://doi.org/10.1016/j.tree.2014.10.009
Wilkening JL, Ray C, Varner J (2015) Relating sub-surface ice features to physiological stress in a climate sensitive mammal, the American pika (Ochotona princeps). PLoS ONE 10:e0119327. https://doi.org/10.1371/journal.pone.0119327
Williams AP, Seager R, Abatzoglou JT, Cook BI, Smerdon JE, Cook ER (2015) Contribution of anthropogenic warming to California drought during 2012–2014. Geophys Res Lett 42:6819–6828. https://doi.org/10.1002/2015GL064924
Williams AP, Seager R, Abatzoglou JT, Cook BI, Smerdon JE, Cook ER (2015) Contribution of anthropogenic warming to California drought during 2012–2014. Geophys Res Lett 42:6819–6828. https://doi.org/10.1002/2015GL064924

Download PDF

Journal Publication

published 27 Sep, 2023

Read the published version in PeerJ →

Version 1

posted

You are reading this latest preprint version

Life on the edge - a changing genetic landscape within an iconic American pika metapopulation over the last half century

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Methods

Study species

BSHP Population census surveys

Genetic sample collection

DNA isolation and PCR amplification

Allelic dropout and genotyping error rates for fecal and museum samples

Population genetic analyses

Results

Patch occupancy

Allelic dropout, null alleles and genotyping error rates

Genetic diversity

Population genetic analyses

Discussion

Population resilience and the role of neutral genetic variation in species’ persistence

Summary

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1