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

Ethical approval and consent to participate

This research was approved by the 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).

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 & 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 & 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 typically 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 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 & 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 & Smith, 1997).

Climate

To examine trends in summer temperature at the study site from 1940 to present day, we used monthly historical PRISM data (Daly, Neilson & Phillips, 1994) extracted, using bilinear interpolation, at the mean location of the 82 habitat patches surveyed in annual censuses. Mean summer temperatures were calculated as the means of daily high and low temperatures for June, July, and August of each year. We used linear regression to assess trend over times and a two-sided Welch Two sample t-test to compare temperatures from 1940–1960 with temperature from 2000–2030. Analyses were conducted in R version 4.3 (R Core Team, 2023).

BSHP population census surveys

Eighty-two individual pika habitat patches, representing the majority of the total suitable habitat at BSHP (~100 patches) were surveyed in 1972, 1977, 1989, 1991–2001, 2003–2006, and 2008–2020 (Smith, 1974a; Moilanen, Smith & Hanski, 1998; Nichols, this study). The presence of fresh fecal pellets and/or active haypiles within each pika territory was used to determine whether each patch was occupied and to estimate the number of pikas occupying each patch (Dearing, 1997; Moilanen, Smith & Hanski, 1998). While all 82 patches were surveyed during later survey years (i.e., after 2008), earlier years have some idiosyncratically distributed missing data (i.e., ~2% missing data, overall). To ameliorate the influence of missing data, we linearly imputed missing occupancy and count data over time (Fig. S1). We used nonparametric tests (i.e., two-sided Spearman’s rank correlation) to evaluate the strength and significance of trends over time for patch occupancy, population count, and number of pikas per occupied patch. We used quasibinomial logistic and negative binomial regressions, respectively for occupancy and population count data, to visualize trends over time. Population count trends were evaluated 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 (see Fig. 2). Patch occupancy and number of pikas per occupied patch were evaluated for both the full period of record (1972–2022) and for the subset of years that excludes the first two survey years (1989–2022). Analyses were conducted in R version 4.3 (R Core Team, 2023).

Genetic sample collection

Mid-20^th century museum samples (1947–1949). Approximately 200 pika museum specimens collected during the mid-20^th century at BSHP by 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 (Supplemental 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) collected by S.B. Benson and O.P. Pearson in 1947 (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 to test for a pattern of increasing regional isolation of the BSHP population over the time frame from 1948–2015.

Late 20^th century samples (1988–1991). Pika ear tissue was collected from live-trapped individuals in BSHP as part of a metapopulation dynamics study of 24 occupied patches in the northern portion of the study area (1988–1990; Peacock & 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 21^st 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 vs 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. The forceps were wiped down between samples, submerged in 70% to 100% ethanol and allowed to air dry between sample collections to remove cellular debris. Two to five fecal pellets collected per individual animal were immediately transferred to a 2.0 ml vial with 500 µl of ATL tissue lysis buffer (DNeasy96 Blood and Tissue kits; QIAGEN, Hilden, Germany) after collection in the field. At the end of each collection period (~2–4 h), 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. In order to rule out multiple individuals in DNA isolated from multiple fecal pellets and to also control for possible contamination during pellet collection, 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 also 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).

DNA isolation and polymerase chain reaction (PCR) amplification

Total genomic DNA was isolated according to the manufacturer’s protocol (DNeasy96 Blood and Tissue kits; QIAGEN, Hilden, Germany) with modifications specific to each sample type. Details regarding modifications to DNA extraction protocol for fecal and museum genetic samples can be found in Hekkala et al. (2011) and Peacock et al. (2017). All DNA extractions and polymerase chain reactions (PCR) for mid-20^th 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, Kirchoff & Merideth, 2002; OCP10, OCP17, OCP21; GQ461705.1–GQ461722.1, 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.2 G thermal cycler (Thermo Fisher Scientific, Waltham, MA, USA), and included a 15 min hot start at 94 °C, followed by 41 cycles at 94 °C, 30 s each for denaturing, a touch down annealing temperature for 90 s (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 s (modified from Peacock, Kirchoff & Merideth, 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, Waltham, MA, USA). Fragment size analysis was carried out on an Applied Biosystems 3730 DNA Genetic Analyzer at the UNR Nevada Genomics Center (https://www.unr.edu/genomics/services). 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 & 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.

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 & 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 eight 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 & Haydon, 2007). In addition, the programs MicroChecker (version 2.2.3; Van Oosterhout & Hutchinson, 2004) and FreeNA (Chapuis & 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, Perrin & Waser, 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 (version 8.0, 1998) together with post hoc Tukey tests. 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, Ovenden & Peel, 2004; Waples & 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.

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 a Bayesian genotype clustering approach (STRUCTURE version 2.3.4; Pritchard, Stephens & Donnelly, 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, Regnaut & Goudet, 2005).

Results from an earlier study which tracked tagged individuals and used variable number of tandem repeat minisatellite nuclear DNA data (Peacock & 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. To assess changes in movement dynamics and population genetic structure we compared the patch level microsatellite datasets for 1988–1991 and 2013–2015 using STRUCTURE (version 2.3.4, Pritchard, Stephens & Donnelly, 2000). To further characterize genetic changes within and among patches between the two time periods, 1988–1990 and 2013–2015, we calculated population level relatedness (r) for each time period using the Queller & Goodnight (1989) method 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 made pairwise comparisons of genetic similarity between patches sampled in 1988–1991 and again in 2013–2015 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 for matched patch pairs 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.

Climate

From 1940 to 2022, mean summer temperature increased by 0.38 °C/decade (95% CI [0.27–0.50], P = 4.7 × 10⁻⁹) (Fig. 3). Mean summer temperatures from 2000–2020 were 2.07 °C (95% CI [1.51–2.63]) higher than mean summer temperatures from 1940–1960 (P = 9.1 × 10⁻⁹).

Patch occupancy

Across the entire patch network, amidst interannual fluctuations, the percentage of patches occupied has declined over time (1972–2022: r_s= −0.47, P = 0.0062; 1989–2022: r_s = −0.36, P = 0.049; Fig. 4A) while the population count has remained relatively unchanged (r_s = −0.15, P = 0.42; Fig. 4C), resulting in an increased concentration of pikas on remaining occupied patches (1972–2022: r_s = 0.45, P = 0.0092; 1989–2022: r_s = 0.39, P = 0.032; Fig. 4D). Data from census surveys show the collapse of the middle and southern pika subpopulations occurring sometime during the 1980s followed by intermittent extirpations and recolonizations from 2008–2022. Population counts in these regions declined significantly over time (middle: r_s = −0.67, P = 0.000017; southern: r_s = −0.52, P = 0.0021). The southern and middle portions of the study area are dominated by small (1–4 territories) and medium (5–9 territories) patches (Fig. 4B) and have never been fully reoccupied and remain largely vacant (Fig. 4C). After collapse in the 1980s, subpopulations in the southern patch network (all patches south of Silver Hill, see Fig. 2) experienced a small increase in the number of patches occupied in 2003 and 2004 (4–5 patches occupied), but occupancy declined to one patch in 2005 before complete extirpation by 2008 (Fig. 4C). Since then, there have been a few recolonization attempts in the southern patch network, but no patch has remained consistently occupied; from 2008–2022 three patches have experienced intermittent occupancy, with up to two patches occupied at a time. The talus patches in the middle patch network (patches on Silver Hill south of Mono Shaft, see Fig. 2) were all mostly unoccupied by the early 2000s and have remained so (Fig. 4C). In contrast, the northern patch network (all patches north of and including Mono Shaft) has maintained consistently occupied patches over the same period and although population size has fluctuated on these patches, the northern network has not experienced a population decline (r_s = 0.06, P = 0.75; Fig 4C).

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.

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. 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 largest 26 patches sampled consistent with a signature of genetic bottleneck (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.0001). 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 × OCP21, P = 0.0001, 7,280 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 × OCP9) and Hobart (OCP4 × OCP8, OCP8 × 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_2,18 = 10.228, df = 2, P = 0.001), and Pa (F_2,18 = 20.89, P = 0.000), and He (F_2,18 = 8.774, P = 0.006), but not Ho (F_2,18 = 0.087, P = 0.917) over the three time periods at the BSHP site (1948–2015; Figs. 5A–5D ). 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), 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). Effective population size also 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. 5E).

The small sample size from 1948–1949 precludes an accurate assessment of 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 to 0.128]) and 1988–1991 (r = 0.024, 95% CI [−0.060 to 0.060]) but an r value consistent with an average relatedness among first cousins (r = 0.158, 9% CI [−0.028 to 0.021]) in 2013–2015 (Fig. 5F). 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 = 1,780.085, df = 2, P = 0.000; Fig. 6A). For the subset of patches for which we had sufficient data from both 1988–1991 and 2013–2015, within patch comparisons among individuals also showed a significant increase in genetic similarity between time periods (Mann-Whitney U, 10.5–97595.5, df = 1, P = 0.000–0.009; Fig. 6B), as well as increases in average relatedness for all patch pair comparisons across years 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. S2). 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, 1,560 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).

Bayesian genotype clustering (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 (Figs. 7A and 7B). 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 7C). These results were also supported by the PCoA analysis, which largely separated the 2013–2015 dataset from the 1947–1949 and 1988–1991 datasets on coordinate axis 1 (Fig. 7D). Whereas the BSHP 1948–1949 and 1988–1991 samples showed considerable overlap in PCoA space and the BIH 1947 samples overlapped completely with 1948–1949 BSHP samples.

Although an assignment test approach, such as Bayesian genotype clustering analysis, is not well suited to capture rapidly shifting population dynamics, STRUCTURE results did reveal an increase in population genetic structure among BSHP subpopulations in 2013–2015 compared to 1988–1991. There was no statistical support for population subdivision in the 1988–1991 dataset. However, in the 2013–2015 dataset, there was statistical support for both k = 2 and 4 (k = 2, average LnP(D) = −2492.92, SD ± 1.540, Δk = 33.59457; k = 4, average LnP(D) = −2,466.82, SD ± 9.583, Δk = 31.943; Fig. S3).