Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Biological consequences of rapid environmental change in the American pika Ochotona princeps Waterhouse, Matthew David 2017

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2017_november_Waterhouse_Matthew.pdf [ 7.13MB ]
JSON: 24-1.0355850.json
JSON-LD: 24-1.0355850-ld.json
RDF/XML (Pretty): 24-1.0355850-rdf.xml
RDF/JSON: 24-1.0355850-rdf.json
Turtle: 24-1.0355850-turtle.txt
N-Triples: 24-1.0355850-rdf-ntriples.txt
Original Record: 24-1.0355850-source.json
Full Text

Full Text

BIOLOGICAL CONSEQUENCES OF RAPID ENVIRONMENTAL CHANGE IN THE AMERICAN PIKA OCHOTONA PRINCEPS  by  Matthew David Waterhouse  B.Sc., University of Maine at Farmington, 2005 M.Sc., University of Wisconsin at Stevens Point, 2011  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  DOCTOR OF PHILOSOPHY  in  THE COLLEGE OF GRADUATE STUDIES  (Biology)     THE UNIVERSITY OF BRITISH COLUMBIA  (Okanagan)  September 2017  © Matthew David Waterhouse, 2017  ii  The undersigned certify that they have read, and recommend to the College of Graduate Studies for acceptance, a thesis entitled:   BIOLOGICAL CONSEQUENCES OF RAPID ENVIRONMENTAL CHANGE IN THE AMERICAN PIKA OCHOTONA PRINCEPS  submitted by Matthew David Waterhouse in partial fulfillment of the requirements of the degree of Doctor of Philosophy.    Michael Russello, Department of Biology, University of British Columbia Okanagan Supervisor, Professor   Karen E. Hodges, Department of Biology, University of British Columbia Okanagan Supervisory Committee Member, Professor   Jason Pither, Department of Biology, University of British Columbia Okanagan Supervisory Committee Member, Professor   Erik Beever, Department of Ecology, Montana State University Supervisory Committee Member, Professor   Kevin Hanna, Department of Earth and Environmental Sciences, University of British Columbia Okanagan University Examiner, Professor   Mary Peacock, Department of Biology, University of Nevada External Examiner, Professor   September 28, 2017 (Date Submitted to Grad Studies)   iii Abstract Species are often confronted with rapid environmental change that require an adaptive response to maintain viability. The source of this environmental change can be natural, but is frequently anthropogenic in nature. In this thesis, I document the biological consequences of rapid environmental change in a sensitive mammal, the American pika (Ochotona princeps). In Chapter 2, I used microsatellite genetic makers to investigate the consequences of landscape modifications for pika populations relative to two major developments in British Columbia, Canada: a large open-pit copper mine (Highland Valley Copper) under partial reclamation and a bisecting major highway (97C). I found evidence for restricted movement associated with human-modified landscapes and showed evidence for fragmentation potentially resulting from the highway. Rapid environmental change often elicits a stress response in animals involving elevated glucocorticoids. In Chapter 3, I applied a novel technique to measure chronic stress in pikas by measuring extracted corticosterone from hair samples. Applying this method along two elevational transects of populations in North Cascades National Park, WA, allowed me to assess individual factors influencing stress in the American pika. My results showed elevated stress levels associated with smaller body sizes, female pikas, and with low early-spring ambient temperatures. This study provided direct physiological evidence for thermal stress in a climate-sensitive mammal and provides a useful tool for assessing climate stress in the American pika. Selective pressure from such thermal stress can result in local adaptation via natural selection. In Chapter 4, I used genotyping-by-sequencing to assess genotype-environment associations along these same elevational transects. The results showed a consistent genomic response to climate conditions across both transects and provide preliminary evidence for cold stress and hypoxia in the  iv American pika. Additionally, I resolved consistent evidence for a downslope bias in gene flow, which may interfere with the upslope movement of individuals that is predicted with ongoing climate change. Taken together, these chapters show consistent patterns of restricted gene flow and thermal stress in the American pika and show the potential that climate-mediated natural selection is resulting in the development of local climate adaptations.  v Preface Chapters 2, 3, and 4 of this thesis are or will be published with the contributions of several co-authors.   Chapter 2 is based on work conducted in collaboration with Cheryl Blair and Karl Larsen at Thompson Rivers University. Cheryl Blair collected all samples and Karl Larsen was responsible for developing the Highland Valley Copper mine as a study site. Michael Russello and I designed the study. I was responsible for genotyping all samples, performing the statistical analysis, and writing the manuscript. A version of Chapter 2 has been published:   Waterhouse MD, Blair C, Larsen K, Russello MA (2017) Genetic Variation and Fine- Scale Population Structure in American pika across a Human-Modified Landscape.  Conservation Genetics.  Michael Russello and I designed the study for Chapter 3. Chris Ray and Jenifer Wilkening contributed data and reference samples for this project. Liesl Erb helped develop study sites in the North Cascades National Park. I collected samples, Bryson Sjodin extracted and quantified corticosterone samples, I conducted all statistical analysis, and wrote the manuscript. A version of Chapter 3 has been published:   Waterhouse MD, Sjodin B, Ray C, Erb L, Wilkening J, Russello MA (2017) Individual-based analysis of hair corticosterone reveals factors influencing chronic stress in the American pika. Ecology and Evolution, 12, 4099–4108.  Michael Russello, Liesl Erb, Erik Beever, and I designed the study for Chapter 4. I collected all samples with the help of field technicians, I conducted all statistical analysis and wrote the chapter. A version of this chapter has been submitted to Molecular Ecology.   vi Table of Contents Abstract ................................................................................................................................... iii	Preface ...................................................................................................................................... v	Table of Contents ................................................................................................................... vi	List of Tables .......................................................................................................................... ix	List of Figures .......................................................................................................................... x	Acknowledgements ................................................................................................................ xi	Chapter 1: Introduction ......................................................................................................... 1	1.1 Climate change .................................................................................................................... 1	1.2 Local adaptation .................................................................................................................. 3	1.3 Tests for selection ............................................................................................................... 4	1.4 Study system ....................................................................................................................... 6	1.5 Thesis objectives ................................................................................................................. 7	Chapter 2: Genetic variation and fine-scale population structure in American pikas across a human-modified landscape ................................................................................... 10	2.1 Background ....................................................................................................................... 10	2.2 Methods............................................................................................................................. 12	2.2.1 Study sites ........................................................................................................................... 12	2.2.2 Sampling ............................................................................................................................. 13	2.2.3 Genetic data collection ........................................................................................................ 14	2.2.4 Site-level analysis ................................................................................................................ 15	2.2.5 Landscape-level analysis ..................................................................................................... 16	2.3 Results ............................................................................................................................... 18	2.3.1 Data quality ......................................................................................................................... 18	2.3.2 Site-level genetic analysis ................................................................................................... 18	2.3.3 Landscape-level genetic analysis ........................................................................................ 19	2.4 Discussion ......................................................................................................................... 20	Chapter 3: Individual-based analysis of hair corticosterone reveals factors influencing chronic stress in the American pika .................................................................................... 32	 vii 3.1 Background ....................................................................................................................... 32	3.2 Methods............................................................................................................................. 35	3.2.1 Sample site and sample collection ...................................................................................... 35	3.2.2 Microclimate measurements ............................................................................................... 36	3.2.3 Molecular sexing ................................................................................................................. 36	3.2.4 Extraction and immunoassay of corticosterone .................................................................. 37	3.3 Results ............................................................................................................................... 40	3.3.1 Laboratory validation .......................................................................................................... 40	3.3.2 Equivalence of corticosterone estimates from different sample sources ............................ 40	3.3.3 Pika stress analysis along elevational gradients .................................................................. 41	3.4 Discussion ......................................................................................................................... 42	3.4.1 Laboratory validation .......................................................................................................... 42	3.4.2 Field study ........................................................................................................................... 45	Chapter 4: Adaptation to climate warming from environmentally-mediated selection may be impeded due to directional gene flow in a thermal-sensitive mammal ............... 58	4.1 Background ....................................................................................................................... 58	4.2 Materials and methods ...................................................................................................... 61	4.2.1 Sample site and sample collection ...................................................................................... 61	4.2.2 Climate measurement .......................................................................................................... 63	4.2.3 DNA extraction and RADseq genotyping ........................................................................... 64	4.2.4 Reference assembly and SNP discovery ............................................................................. 65	4.2.5 Outlier detection and annotation ......................................................................................... 66	4.2.6 Population genetic analysis ................................................................................................. 69	4.2.7 Gene flow analysis .............................................................................................................. 70	4.2.8 Stress hormone analysis ...................................................................................................... 71	4.3 Results ............................................................................................................................... 71	4.3.1 Sample and environmental data collection ......................................................................... 71	4.3.2 Reference assembly and SNP discovery ............................................................................. 72	4.3.3 Outlier detection .................................................................................................................. 73	4.3.4 Population genetic analysis ................................................................................................. 75	4.3.5 Gene flow analysis .............................................................................................................. 76	4.3.6 Stress hormone analysis ...................................................................................................... 76	 viii 4.4 Discussion ......................................................................................................................... 77	4.4.1 Outlier detection .................................................................................................................. 78	4.4.2 Outlier annotation ................................................................................................................ 80	4.4.3 Population genetic analysis ................................................................................................. 81	4.4.4 Directional migration .......................................................................................................... 83	4.4.5 Stress hormone analysis ...................................................................................................... 84	4.4.6 Summary ............................................................................................................................. 85	Chapter 5: Conclusion ........................................................................................................ 102	5.1 Research findings ............................................................................................................ 102	5.1.1 Evidence for restricted dispersal ....................................................................................... 102	5.1.2 Evidence for cold stress .................................................................................................... 104	5.1.3 Metabolic implications ...................................................................................................... 105	5.2 Limitations ...................................................................................................................... 106	5.2.1 Sampling limitations ......................................................................................................... 106	5.2.2 Genetic limitations ............................................................................................................ 107	5.2.3 Geographic limitations ...................................................................................................... 108	5.3 Management implications ............................................................................................... 109	5.4 Significance and future directions .................................................................................. 110	References ............................................................................................................................ 112	   ix List of Tables  Table 2.1. Sample locations, approximate patch size (m2), sample sizes (n), and genetic diversity metrics for the 15 sampling sites ................................................................. 26 Table 2.2. AMOVA results showing the distribution of genetic variation explained by different hierarchical classification schemes .............................................................. 27	Table 2.3. Pairwise site comparisons of genetic differentiation (θ) for American pika within and around Highland Valley Copper and highway 97C. ............................................ 28	Table 3.1. Site description of the Pyramid Peak (PP) and Thornton Lakes (TL) sampling transects in North Cascades National Park, WA. ....................................................... 50	Table 3.2. Correlation matrix of independent variables ......................................................... 51	Table 3.3. Information-theoretic analysis of mixed-effects models explaining variation in corticosterone estimates of American pika samples ................................................... 52	Table 4.1. Summary of sample sizes (n), observed heterozygosity (Ho), total corrected heterozygosity (He), effective number of alleles (Ae), and inbreeding estimates (Fis) for sample sites along the Pyramid Peak (PP) and Thornton Lakes (TL) transects in the North Cascades National Park, WA. ..................................................................... 87 Table 4.2. Number of SNPs retained after various filtering procedures. ................................ 88	Table 4.3. Summary of outlier detection analyses. ................................................................. 89	Table 4.4. Summary of top BLASTN search results for outlier loci detected by BAYESCAN (BS), LFMM (LF), and BAYPASS (BP) ...................................................................... 90	Table 4.5. Pairwise comparisons of Fst (Top diagonal; Weir & Cockerham 1984) and associated p-value from a permutation analysis ......................................................... 91	Table 4.6. Directional migration analysis for pairwise movement between sites along the PP and TL transects using either the neutral or outlier datasets. ...................................... 92	Table S4.1. Climate variables from i-Button data obtained along the elevational transects established in the North Cascades National Park, WA. .............................................. 97 Table S4.2. ClimateWNA data for each site along the elevational transects established in the North Cascades National Park, WA. ........................................................................... 98	Table S4.3. Gene Ontology terms associated with significant outlier loci Blast hits. ............ 99	   x List of Figures Figure 2.1. Map of study area indicating sampling locations within and around the Highland Valley Copper mine in south central British Columbia .............................................. 29	Figure 2.2. Photograph depicting typical artificial habitat (a) contrasted with natural habitat (b) in the Highland Valley Copper mine. .................................................................... 30	Figure 2.3. STRUCTURE bar plots averaged over 25 iterations showing the genetic division (a) between sites north and south of highway 97C (ΔK = 507.9) and (b) in the south only when analyzed independently (ΔK = 73.7) ................................................................. 31	Figure 3.1. Photograph of an American pika and reference hair sample (inset) weighing approximately 10 mg .................................................................................................. 53	Figure 3.2. Sample sites in North Cascades National Park, Washington, USA. Thornton Lake (TL) and Pyramid Peak (PP) sampling sites shown as circles .................................... 54	Figure 3.3. Top: Parallelism between the standard curve (solid line with circles) and serial dilutions of one sample (squares, no line) .................................................................. 55	Figure 3.4. Relationship between sample mass and estimated corticosterone concentration using NOCA samples (squares), hair samples from paired plasma and fecal samples (triangles), and differing masses from PP04T08 (circles) .......................................... 56	Figure 3.5. Box and whisker plot showing average corticosterone per site after correcting for extraction efficiency .................................................................................................... 57  Figure 4.1. Map of Pyramid Peak (PP) and Thornton Lakes (TL) transects in the North Cascades National Park, Washington, USA ............................................................... 94	Figure 4.2. Summary plots for genotype-environment associations in American pika .......... 95	Figure 4.3. Admixture proportions for American pika along Pyramid Peak (panels A and C) and Thornton Lakes (panels B and D) transects using only neutral (panels A and B) or outlier loci (panels C and D) with inset cross-entropy (CE) estimations for each value of K ............................................................................................................................. 96 Figure S4.1. Climate variables and loadings for the first (PC1) and second principal component (PC2) of the environmental PCA for the Pyramid Peak and Thornton Lakes transects .......................................................................................................... 100	Figure S4.2. Admixture plot for all the sites sampled in NOCA .......................................... 101	  xi Acknowledgements This research would not have been possible without the support, suggestions, critical feedback, field assistance, analytical advice, and software assistance from my advisory committee, fellow graduate students, and researchers at the University of British Columbia Okanagan and staff at the North Cascades National Park.   I am particularly indebted to my advisor, Mike Russello. He has taught me the results of putting research first and the benefits of considering the evidence behind any claim. Mike’s desire to contribute to our understand of ecological genetics for the betterment of conservation is inspirational. I am thankful to the rest of my committee members for their assistance, feedback and support. Karen Hodges has been a mentor and role model for me since I started my program. Jason Pither has been a great source for academic support. Erik Beever’s passion for ecology and life continue to be motivational.   I am thankful to the current and past members of the Ecological and Conservation Genomics Lab, especially Evelyn Jensen, Bryson Sjodin, Brett Ford and Andrew Veale for their support, assistance and understanding. I am thankful to all the wonderful friends I have met during the course of my program.   I am grateful to my family for their love and support. I have always found the love and devotion my sister, Melinda, shows her family a source of wonder and inspiration. My sister, Christine, always goes the extra mile to help and has been a constant source of strength for me. I get my predisposition for analysis from my father. My mother has always been my beacon for compassion and understanding.   This research was funded by an NSERC Discovery Grant and Seattle City Light. I was partially supported by University Graduate Fellowships from UBC. xii Dedication   For my nephews, Zachary; your natural curiosity is amazing, never stop asking why; Owen, you’re the strongest little man I know; and my niece, Allison, you’re a firecracker, the world needs more of your character. My hope is that through diligent effort we can leave a bright future for kids like you.   1 Chapter 1: Introduction 1.1 Climate change Anthropogenic climate change is causing profound and ubiquitous alterations to natural environments (Stocker et al. 2013). An estimated 350 billion tons of carbon dioxide have been emitted by human activities since 1959 (Ballantyne et al. 2012), resulting in a 0.72°C increase in global temperature (Stocker et al. 2013). Long-term climate modeling indicates that the past three decades were the warmest of the past millennium and that the current rate of warming is far greater than previous natural rates (Jones et al. 2001; Pacifici et al. 2017). While some action has been taken, the moral complexity of the problem (Gardiner 2012) and necessity of a global consensus on action makes abatement of climate-change unlikely. Furthermore, the repercussions of additional carbon already added to the atmosphere will last for at least a millennium (Solomon et al. 2009). The resulting environmental perturbations from climate change are extensive and include sea-level rise, ocean acidification, increased ocean and land surface temperatures, changes in precipitation, and an increase in extreme weather events (Stocker et al. 2013).    These rapid alterations to the environment are resulting in dramatic ecological and physiological stresses on natural populations. The direct ecological impacts of climate change have been documented on every continent, in every ocean, and in most major taxonomic groups (Parmesan 2006). One of the major influences of climate change will be the direct physiological stress of rising temperature (Somero 2010; but see Cahill et al 2013). Many species are adapted to a narrow range of temperatures (stenothermy), and exposure to elevated temperatures can reduce metabolic efficiency or, in extreme cases, mortality can result from acute heat shock. To maintain viability, impacted species must either disperse,  2 adapt, or acclimate (Lovejoy & Hannah 2005; Beever et al. 2017).  Species are expected to undergo poleward or altitudinal range shifts in the face of climate change to maintain optimal thermal conditions (Root et al. 2003). Parmesan & Yohe (2003) documented this expected range shift in 80% of the 434 species surveyed in a meta-analysis, with an average poleward range shift of 6.1 km per decade or an increase in elevation of 6.1 m per decade. A more recent meta-analysis (Chen et al. 2011), reported median shifts nearly 2 to 3 times greater than these estimates, indicating species responses to climate change are much greater than previously expected. Additionally, the focus on unidirectional movement ignores complex interactions among temperature, precipitation and species-specific tolerances, and therefore may underestimate the influence of climate change (VanDerWal et al. 2012). Consequently, some species are disproportionally affected, such as European butterflies, which showed a 200 km northward range shift (Parmesan et al. 1999). The majority of these range shifts likely result from the extirpation of low-elevation or equatorward populations and expansion of poleward or high-elevation populations (Davis et al. 2005).  The ability of a species’ range to shift in response to climate change may be significantly impaired by biotic and abiotic factors. Interactions with native species may limit the ability of dispersing species to occupy new niches (Urban et al. 2012). Natural geographic boundaries can limit the poleward or altitudinal dispersal of many species, especially alpine specialists. Anthropogenic impediments to dispersal may impose the most serious threat to the range shift of species in response to climate change (Hoffmann & Sgrò 2011). Habitat fragmentation, in particular, is a significant threat to biodiversity (Wilcox & Murphy 1985; Opdam & Wascher 2004), and  can reduce dispersal, depending on species- 3 specific behavior (Baguette et al. 2003; van Dyck & Baguette 2005). These natural and anthropogenic limitations imposed on dispersal and colonization may severely limit range shifts in response to climate change, making local adaptation increasingly important.     1.2 Local adaptation Biotic response to the local environment can come from either phenotypic plasticity or local adaptation. Phenotypic plasticity is the non-genetic alteration of an organism’s behavior, morphology, physiology, or development in response to environmental change. Some studies suggest that phenotypic plasticity could play an important role in population acclimatization to elevated temperatures (Hoffmann & Sgrò 2011). Phenotypic change can happen rapidly in natural populations, especially as a result of anthropogenic disturbances (Hendry et al. 2008). For example, warmer temperatures have led to earlier spring emergence in yellow-bellied marmots, Marmota flaviventris, and earlier weaning of young (Ozgul et al. 2010). The resulting lengthening of the growing season for the marmots led to increased body masses and a decline in adult mortality. This type of phenotypic plasticity may provide a short-term mechanism to cope with climate change; however, this response is unlikely to provide long-term solutions to populations experiencing continued directional selection from climate change (Gienapp et al. 2008).  Climate change is likely to increase selective pressures on populations, which may lead to a genetic response in impacted species (Hoffmann & Sgrò 2011). Directional selection from increased drought or thermal stress may drive adaptation in affected populations. Selection experiments have provided evidence for the rapid evolution of climate-related characteristics (Reusch & Wood 2007). For example, Drosophila  4 melanogaster evolves increased heat shock tolerance when subjected to artificially elevated temperatures (Cavicchi et al. 1995). Equally compelling, heritable morphological adaptations induced in D. melanogaster afforded a twofold increase in desiccation resistance after selective pressures were applied to 26 generations (Telonis-Scott et al. 2006). While these studies provide examples of artificially induced climate adaptations, documentation of genetic adaptations to climate change in natural settings has been more difficult (Merilä & Hendry 2014). Rigorous observations have revealed that the pitcher-plant mosquito, Wyeomyia smithii, has undergone genetic adaptation to alter its phenology with altered seasonality in response to recent climate change (Bradshaw & Holzapfel 2001). However, other studies show either a lack of genetic response to climate change (Gienapp et al. 2008) or the lack of adaptability of the species. For example, populations of Drosophila birchii appear incapable of adapting to drier conditions due to low levels of genetic variability for desiccation resistance and will likely go extinct in the drier conditions predicted to occur with climate change (Hoffmann et al. 2003). This lack of adaptability in D. birchii indicates that impacted species will likely require preexisting genetic variation for evolutionary responses to match climate change (Davis et al. 2005; Hoffmann & Sgrò 2011).   1.3 Tests for selection Current genomic technologies allow for unparalleled insight into evolutionary processes occurring in natural populations (Ellegren 2014). The application of next-generation sequencing to non-model organisms continues to shed light on the inner workings of adaptive evolution. Although the precipitous drop of sequencing costs has facilitated the application of direct sequencing to conservation-related questions, whole genome sequencing is still  5 prohibitively expensive and generally superfluous for addressing ecological questions. Reduced-representation techniques have facilitated the application of next-generation sequencing technology for ecological applications by providing high-quality genomic datasets that are still computationally feasible to analyze. One such technique, restriction-site associated DNA sequencing (RADseq), has become increasingly popular due to the ease of application and cost effectiveness for identifying and genotyping of thousands of single nucleotide polymorphisms (SNPs) randomly distributed throughout the genome of organisms with or without a reference genome (Baird et al. 2008).  Analytical methods to detect candidate loci under selection from this type of genomic data generally scan for excessive divergence (FST-based analysis) at specific loci or attempt to resolve genomic correlations to environmental parameters (spatially explicit analysis; Rellstab et al. 2015). FST-based analyses scan for allelic frequencies that diverge from neutral patterns as a potential result of natural selection (Beaumont & Balding 2004; Excoffier & Lischer 2010; Helyar et al. 2011; Günther & Coop 2013). Spatially explicit analysis methods can identify informative correlations between allele frequency and environmental data (Joost et al. 2007; Frichot et al. 2013; Jones et al. 2013). The resulting “outlier” SNPs from either of these two methods may not be directly under selection but instead may be linked to adaptive loci that are providing a selective advantage (i.e., genetic hitch-hiking; Smith and Haigh 1974). These approaches are powerful for teasing apart locus-specific and genome-wide patterns in natural populations (Luikart et al. 2003). Furthermore, combining these two independent analytical approaches has the benefit of reducing the rate of false positives when tests occur over thousands of loci (Rellstab et al. 2015).    6 1.4 Study system The family Ochotonidae contains 30 species of pikas, two of which are native to North America: the Collared pika, Ochotona collaris, and the American pika, Ochotona princeps (Hoffman and Smith 2005). The American pika is an alpine specialist distributed from central British Columbia, Canada south to the Sierra Nevada Mountains in California, USA and east to the eastern edge of the Rocky Mountains from New Mexico to Alberta (Smith & Weston 1990). Five discrete lineages of American pikas associated with separate mountain ranges were reconstructed using mitochondrial DNA from the cytochrome b/D-loop (Galbreath et al. 2009). The Cascade lineage was likely the first to diverge and extends along the Cascade Mountain range from southern British Columbia to northern California.  The geographic range and habitat requirements of pikas are influenced by the species’ thermal sensitivity (MacArthur & Wang 1974; Smith 1974a). Pikas inhabit talus slopes ranging in elevation from sea level to 3,000 m in the northern range and are rarely found below 2,500 m near the southern limits of their distribution (Smith and Weston 1990; but see Millar and Westfall 2010 and Varner and Dearing 2014). Talus provides a cool, moist microclimate that acts as a refuge from the summer heat (Beever et al. 2008); for example, temperatures can be as much as 32°C cooler than ambient temperatures just 1 m into the talus in the Columbia River Gorge, OR (Varner & Dearing 2014). Winter temperatures are also moderated by the insulation of talus and snowpack (Simpson 2009). American pikas use these microclimates to thermoregulate behaviorally (Jeffress et al. 2013; Rodhouse et al. 2017). MacArthur & Wang (1974) found that pika strictly regulate the intensity and duration of activity to maintain their body temperatures just 2-3°C below upper lethal temperatures, albeit this study may have suffered from a low sample size (6 pikas).  7 Although the direct physiological evidence for thermal sensitivity in pikas is limited (Smith 1974), several lines of evidence indicate that climate change has already had a negative impact on American pika populations. Paleontological analysis of pikas in the Great Basin shows that the minimum inhabitable elevation has risen 150 m over the past century (Grayson 2005; Beever et al. 2011). Furthermore, recent extirpations of pika populations in the Great Basin (Beever et al. 2010) and marked declines in abundance (Beever et al. 2013) have been associated with regional climate change. Likewise, Galbreath et al. (2009) found genetic evidence of recent demographic declines in all five lineages of the American pika. Niche modeling projects that American pikas in the USA could disappear from 55% of their current range even under low-carbon emission scenarios (Ray et al. 2012). Concerns over the impacts of climate change on pika populations have spurred consideration of the species for listing as threatened or endangered at both federal and state levels (USFWS 2010; Osborn and Applebee 2011). Moreover, the thermal sensitivity of O. princeps makes it a useful mammalian model for testing the ecological and genetic impacts of climate change. In this light, the America pika may serve as an early warning sign of climate change (Ray et al. 2012) and may inform conservation efforts seeking to abate the loss of biodiversity in changing environments.   1.5 Thesis objectives This thesis consists of three complementary approaches to assessing the biotic consequences of rapid environmental change. In Chapter 2, I apply population-genetic analyses to investigate the response of pikas to major anthropogenic habitat alterations in the form of open-pit mining and potential fragmentation induced from road construction. I demonstrate  8 that the resulting genetic structure of pikas colonizing the newly created habitat from mining activity is likely influenced by human-made impediments to gene flow across the landscape. I show several lines of evidence that a major highway (97C) bisecting the sample site has likely reduced connectivity to the point where populations delineated by the highway represent distinct genetic units.   In Chapter 3, I validate and apply a technique for measuring corticosterone from pika hair samples and show the ecological application of this method by conducting an individual-based analysis of chronic stress. This chapter represents the first application of this technique to pika and potentially the first direct physiological evidence of climate-induced physiological stress in this species. I show how smaller body size increases stress in pikas likely resulting from an increase in the mass-specific metabolic rate of smaller individuals. Additionally, I show higher stress levels in female pikas and demonstrate the potential for cold stress during the early spring molting period by establishing a relationship between increased stress levels and lower ambient temperatures.  Chapter 4 focuses on the robust identification of genomic correlates to environmental parameters along elevational gradients. In this chapter, I use three distinct analyses for detecting outlier loci along two independent elevational transects to identify potential targets of climate-induced natural selection in the American pika. The results of these analyses provided consistent evidence for numerous loci that have been influenced by climatic patterns. Furthermore, I highlight several potential genomic regions and functional mechanisms that may be a source of thermal adaptation in the American pika.   These findings highlight several important biotic responses to rapid environmental change in a species of conservation interest. More specifically, I show that anthropogenic  9 habitat modification can alter the connectivity of American pika populations, potentially limiting the species’ range shift in response to climate change and increasing the importance of local adaptation to climate. Next, I show that climate can have measurable impacts on physiological stress in the American pika highlighting the role climate could be playing as a selective force in pikas. Finally, I reveal direct genomic correlates to microclimate conditions, which may represent thermal adaptations resulting from climate-induced natural selection. This body of research presents evidence for restricted gene flow, environmental stressors, natural selection and adaptive genetic diversity in the American pika. Taken together, these patterns represent the necessary conditions for local adaptation to occur, which may be especially critical during the current period of rapid environmental change.  10 Chapter 2: Genetic variation and fine-scale population structure in American pikas across a human-modified landscape  2.1 Background There is a diverse array of research and knowledge on how habitat alterations influence wildlife populations (Saunders et al. 1991; Keyghobadi 2007). Fragmentation of habitat impedes dispersal (Baguette et al. 2003; Buchmann et al. 2013), having disproportionate effects on some taxa, such as small mammals (Sauvajot et al. 1998). Reduced dispersal can hinder metapopulation dynamics (Fischer & Lindenmayer 2007), leading to reduced resilience in the face of ecological stress. Furthermore, barriers to dispersal may limit the potential of a species to shift its range in response to climate change (Parmesan & Yohe 2003). Habitat modifications leading to reductions in population size can also negatively impact genetic diversity (Frankham 1996), which may influence a species’ ability to adapt to changing environments (Sgrò et al. 2011).  Reclamation activities are increasingly applied to disturbed landscapes in an attempt to improve the suitability of the land for wildlife, with or without a specific target species in mind (Ruiz-Jaen & Aide 2005). For example, many mining operations are now heavily involved in reclamation, often with the end goal of mitigating the initial disturbance and improving habitat recolonization by plants and animals (Eaton et al. 2014). Such activities could represent an important mechanism for reducing impacts to biodiversity on multiple scales. However, we have little direct evidence of the potential consequences colonizing artificial habitat may have on the resulting demographic structure of wildlife populations on these landscapes. In particular, the juxtaposition of local, artificial habitat and neighboring  11 natural habitat may influence the extent and direction of gene flow. Moreover, anthropogenic landscape features (e.g., roadways) commonly associated with these areas may lead to further habitat fragmentation. Understanding how ecological processes act within reclaimed landscapes is an important step forward in our efforts to improve the resiliency of wildlife populations.  Species with narrow habitat requirements that are able to colonize human-modified environments may provide good opportunities for investigating the relationship between demography, population genetics, and landscape modifications. The American pika (Ochotona princeps) is a small lagomorph that has shown a narrow tolerance range for ambient temperature (Smith 1974b; Hafner & Sullivan 1995; Beever et al. 2010; Stewart et al. 2015). This species is patchily distributed in rocky, talus-type habitats across mountainous areas throughout western North America from central British Columbia and Alberta, Canada, south to the Sierra Nevada in California and east to New Mexico, USA. In some instances, American pikas have colonized reclaimed mining landscapes, most notably in Bodie, California, USA (Peacock and Smith 1997a). The fragmented nature of their habitat and limited dispersal ability (Henry et al. 2012; Castillo et al. 2014; Robson et al. 2016) has increased focus on the American pika for studies of metapopulation dynamics, island biogeography, and source-sink dynamics (Peacock & Smith 1997a; Moilanen et al. 1998; Kreuzer & Huntly 2003; Beever et al. 2013).  In this study, we investigated the genetic consequences of landscape modifications on the American pika relative to two major developments in British Columbia, Canada, a large open-pit copper mine under partial reclamation and a bisecting major highway. We collected microsatellite genotypic data for individuals sampled at sites within and adjacent to the mine  12 both north and south of the highway, and employed site- and landscape-level analyses to quantify levels of variation and connectivity across this human-modified landscape.   2.2 Methods 2.2.1 Study sites Highland Valley Copper (HVC) is located approximately 54 km southwest of Kamloops, BC (Figure 2.1a). The original mine was commissioned in 1962, although mineral explorations in the area date back to 1954. Originally, three mines operated in the Highland Valley. In 1986, they were amalgamated into one mine, which is now one of the world’s largest open-pit copper mines. Currently, mining operations occupy approximately 6,200 ha (Freberg & Gizikoff 1999). Surface developments from the operation of the mine include open pits, waste rock dumps, tailings, infrastructure, water diversions, and roads. Several natural and anthropogenic features punctuate the landscape. The Highland Valley runs east to west through the study site representing an approximate 300 m change in elevation with a seasonal stream (Witches Brook) at the bottom. Additionally, Highway 97C was completed in 1990, and runs along the bottom of this valley and bisects the HVC mine. Approximately 1,320 vehicles use this highway daily, with peak hours between 5-8 am and 4-8 pm (British Columbia Ministry of Transportation and Highways 2009). Major mining operations in the section of HVC north of the highway ceased in 1982 while the southern section remains active. Extensive reclamation activities at HVC were initiated in the late 1980s (Freberg & Gizikoff 1999) and have largely been structured around revegetation and lake remediation. Revegetation goals include: the establishment of forage for cattle, native shrubs and trees for wildlife browse, and conifers for wildlife corridors (Bloodgood et al. 1998). More recently, reclamation plans have been modified to increase  13 biodiversity of the mine (Teck Resources Limited 2012). These plans have never specifically identified American pikas as a target species, but the presence of pikas in the reclaimed landscape was first observed by mine workers around 2005, and then formally documented by Howie (2007). The closest natural population of American pikas is adjacent to the mine property within 0.5 km. The elevation of pika-occupied sites in the mine (1350-1550 MSL) is comparable to occupied surrounding sites (1350-1850 MSL).    2.2.2 Sampling  This study is part of a larger project investigating the population ecology of pikas both within and near the HVC operating area (Blair & Larsen unpubl.). Sampling locations were selected on both natural (n = 7) and artificial (n = 8) habitat north and south of Highway 97C based on site occupancy. Pika-occupied sites were initially located by using aerial mapping and local knowledge to identify rocky areas of both natural and artificial habitat. Intense searching on the ground then was conducted to identify individual pika territories through direct observation of the animals and/or their hay piles. Pika territories were considered to be in artificial habitat when in waste rock dumps or riprap from road construction (Figure 2.2a). All other sites were considered ‘natural’ and represented talus patches in a relatively undisturbed state (Figure 2.2b). Sites were exhaustively trapped between April and October 2012 and between April and November 2013, allowing for the identification of individual pika territories. Pikas were captured using Tomahawk model 202 (Hazelhurst, WI) collapsible live traps in accordance with BC Ministry of Forests, Lands and Natural Resource Operations wildlife permits KA12-78714 and KA13-86652, and Animal Ethics Protocol #100102 (Thompson Rivers University, BC). All habitat patches within a 500 m radius were  14 considered part of the sample site and were measured using a range finder and tape measure to approximate the total area of habitat for each site. A GPS coordinate was taken at each individual territory and all GPS points within the 500 m radius were averaged to obtain site coordinates.  To determine the age class of captured animals, we combined observational data with estimates of mass using a spring scale and cranial diameter using calipers. Individuals with a mass under 150 g and cranial diameter under 2.5 cm were categorized as juveniles (informed by Smith and Weston 1990); all such individuals were generally trapped emerging from their natal nest and were considered young of the year. To eliminate resampling individuals, each pika was marked with two Monel #1 ear tags (unique number combinations) and one unique color tag. A small tuft of hair was plucked and stored in a coin envelope for subsequent genetic analysis.   2.2.3 Genetic data collection  DNA was extracted using the Macherey-Nagel NucleoSpin Tissue kit (Macherey-Nagel GmbH & Co. KG, Duren, Germany) and manufacturer’s protocols. Eleven polymorphic microsatellite loci were used to genotype each sample (Supplemental Table 2.1; Peacock et al. 2002; Peacock and Kirchoff, unpublished report). Conditions for PCR amplification followed Henry et al. (2012) including one additional locus (Ocp 10) that was not previously used, but was amplified under the same touchdown PCR protocol. PCR products were co-loaded and run on an ABI 3130XL DNA automated sequencer (Applied Biosystems, Foster City, CA) with GENESCANTM 500 LIZ® size standard. Genotype calls were made using GENEMAPPER 4.0 (Applied Biosystems, Foster City, CA). To assess allele-scoring error, 40%  15 of the samples were re-amplified and re-genotyped independently and compared to the original scores.  Sex was determined for each sample by the selective co-amplification of an allosomal-linked locus (SRY) and an autosomal control locus (Ocp 10) as described by Lamb et al. (2013). Scoring was conducted by running the PCR product on a 1.5% agarose gel containing 2.5% SYBR Safe (Invitrogen, Carlsbad, California). The genotypic data were examined for evidence of large allele dropout and null alleles using MICROCHECKER (van Oosterhout et al. 2004). All loci were tested for deviations from Hardy-Weinberg expectations (HWE) in each sample site using an exact test implemented in GENEPOP 4.0 (Raymond and Rousset 1995; Rousset 2008). Linkage disequilibrium was tested between all pairs of loci in each site using the exact test of Guo and Thompson (1992) as implemented in GENEPOP 4.0. Type I error rates for tests of linkage disequilibrium and departure from HWE were corrected for multiple comparisons using the sequential Bonferroni method (Rice 1989).  2.2.4 Site-level analysis Sex ratios (M:F) were calculated for each site based on the molecular sexing data. We tested for even sex ratios using a chi-squared (χ2) goodness-of-fit test and the chisq.test function in R version 3.3.1 (R Core Team 2015) between the number of males and females for all artificial and natural sites, respectively. Unbiased expected heterozygosity (He) was calculated for each population using ARLEQUIN 3.5 (Excoffier & Lischer 2010). Allelic richness (AR) was estimated using a rarefaction method described by Leberg (2002) to account for biases caused by unequal sample sizes as implemented in HP-RARE v1.0  16 (Kalinowski 2005). The inbreeding coefficient (Fis) was calculated for each site and tested for statistical deviations from zero using 10,000 permutations in GENETIX (Belkhir et al. 2004). Pairwise relatedness was calculated between all samples at each site using the estimator developed by Queller and Goodnight (1989) and tested for significance using a permutation test with 1,000 replicates in GENALEX (Peakall & Smouse 2006). Heterozygosity, rarefied allelic richness, site-level relatedness, and Fis were compared between natural and artificial sites using a two-tailed t-test assuming unequal variances using the t.test function in R.  2.2.5 Landscape-level analysis The presence of discrete genetic units was assessed using a Bayesian clustering method implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000). An admixture model with correlated allele frequencies was used with a run length of 1,000,000 MCMC replicates after a burn-in period of 500,000. The most likely number of clusters (K) was determined by varying K from 1 to 17 with 25 iterations per value of K and implementing the ΔK method (Evanno et al. 2005) using STRUCTURE HARVESTER (Earl & VonHoldt 2011). Additional population structure was assessed by re-running the STRUCTURE analysis for each resolved genetic unit separately using the same parameters but varying K from 1 to the number of sampling sites in the resolved genetic unit plus 2. An Analysis of Molecular Variance (AMOVA) was performed in GENALEX using the resolved genetic groups from the STRUCTURE analysis in addition to the natural/artificial classification to determine the degree of genetic differentiation explained by these groupings and tested for significance using a permutation test with 1,000 replicates.  The level of genetic differentiation between pairwise comparisons of sites was  17 estimated using θ (Weir & Cockerham 1984) and tested for significance using 10,000 permutations as calculated in GENETIX (Belkhir et al. 2004) and corrected for multiple comparisons using the false discovery rate correction (Benjamini and Hochberg 1995). To test for a pattern of isolation-by-distance (IBD), a matrix of genetic distances (θ) was compared to a matrix of Euclidean distances using a Mantel test implemented in the ISOLATION BY DISTANCE WEB SERVICE (Jensen et al. 2005) with default parameters and tested for significance using 1,000 permutations. This analysis was repeated using the same parameters, but using only sites north or south of highway 97C, respectively. Additionally, to examine the possibility of spatially variable IBD patterns, we first constructed a genetic similarity matrix of (1- θ) between all pairwise comparisons of sites. We then used this genetic similarity matrix to construct a non-stationary genetic friction map displaying the relative genetic divergence per unit of geographic distance as implemented in LOCALDIFF (Duforet-Frebourg & Blum 2014) using 4 simulated neighbors at a distance of 0.1 and 100 posterior replicates.  The direction and magnitude of contemporary gene flow were assessed using a non-equilibrium Bayesian method implemented in BAYESASS v. 3 (Wilson & Rannala 2003). A run length of 10,000,000 MCMC replicates with a burn-in period of 1,000,000 replicates was used, sampling the chain every 100 iterations. To evaluate consistency, the program was run five times with a different random seed. Given the recovered structure in the dataset (see Results below), we grouped sites for this analysis as follows: north natural, north artificial, south natural, and south artificial. Significance was assumed when the 95% credibility set [mean ± 1.96 × standard deviation (sd)] did not encompass zero, as recommended by the developers (Wilson & Rannala 2003).  18 2.3 Results 2.3.1 Data quality A total of 109 pikas were sampled from 15 sites, 8 artificial and 7 natural, with an average sample size of 7.3 animals per site (Table 2.1). Our habitat assessment showed an average of 4,117 m2 of habitat per site with a mean pika territory of approximately 590 m2. These parameters agree with previous reports (Smith and Weston 1990) indicating our sampling was likely representative of the total population at these sites. Of these individuals, 37 were unambiguously assigned as juveniles, 66 were assigned as adults, and the remaining 6 were indeterminate. The overall dataset contained 1.0% missing genotypic data, with no sample having missing data at more than 3 loci. Independently genotyping 47 random samples showed a 1.4% allelic scoring error rate, which is within reasonable expectations for the use of hair as a genetic source material and is not expected to skew population genetic analyses (Smith & Wang 2014). There was no evidence for large allele dropout or null alleles at any sampling sites with the exception of Ocp 15, which showed evidence of null alleles at four sites (HFGR, Relic_1, SGS, and FHF). Following sequential Bonferroni corrections, no locus deviated from HWE or showed evidence of linkage disequilibrium. Given the lack of a consistent trend in the evidence for null alleles across sites, all loci were retained for further analysis.  2.3.2 Site-level genetic analysis The average sex ratio (M:F) was 1.9 (sd = 1.5) for artificial sites and 1.3 (sd = 1.0) for natural sites (Table 2.1). There was no significant difference between the number of males and females on either natural sites (χ 2 = 5.43, df = 7, p = 0.61) or on artificial sites (χ2 = 5.00, df = 6, p = 0.54). Site-level heterozygosity of pikas ranged from 0.530 (EC) to 0.707 (BRCH),  19 while AR ranged from 2.14 (SGS) to 2.58 (BRCH). Site-level relatedness ranged from 0.028 to 0.394 and was significantly greater than expected under random mating for 11 of 15 sites (Table 2.1). Inbreeding estimates ranged from -0.179 to 0.270; 2 of the 15 sites exhibited inbreeding estimates significantly greater than 0. There was no difference in mean genetic diversity estimates between natural and artificial sites for He (t = 0.965, df = 7.694, p = 0.82), AR (t = 1.816, df = 9.996, p = 0.95), or Fis (t = 1.034, df = 11.628, p = 0.84). Relatedness of pikas was significantly higher at artificial sites (mean = 0.251, sd = 0.098) than natural sites (mean = 0.146, sd = 0.071; t = -2.396, df = 12.580, p = 0.016).  2.3.3 Landscape-level genetic analysis  The STRUCTURE analysis revealed evidence for K = 2, corresponding to clusters of sites north and south of highway 97C (ΔK = 507.9; Figure 2.3a). Further analysis did not reveal additional genetic units in the north, but resolved three southern genetic units (ΔK = 73.7; Figure 2.3b). The two natural sites (Relic_1 and Relic_2) grouped together with an artificial, admixed site (SGE). Two additional artificial sites (HNR and HFGR) grouped together, while one artificial site (SGS) formed a largely distinct cluster from the other southern sites. The AMOVA results showed the north/south divide of highway 97C consistently explained the greatest amount of genetic variation, and the natural/artificial demarcation explained the least (Table 2.2). While all grouping scenarios were significant, the four genetic clusters resolved by the STRUCTURE analysis were the most explanatory followed by the natural/artificial grouping when the north/south divide was added to the hierarchical structure (K = 4).     Significant genetic differentiation was found for 35 of 105 of the pairwise site comparisons of θ (Table 2.3). A weak but significant pattern of IBD was also detected (r2 =  20 0.080, p = 0.009) where genetic distance increased with geographic distance when considering all sites. This pattern did not hold when only considering sites north (r2 = 0.003, p = 0.402) or south (r2 = 0.087, p = 0.122) of highway 97C. The genetic friction map resolved a localized area of disproportionately high genetic differentiation in the central region of the study area largely corresponding with landscape modification associated with HVC and highway 97C (indicated in red; Figure 2.1b). There was no genetic evidence of migration between the northern and southern genetic units. Within these genetic units, 7.6% (sd = 3.8%) and 24.8% (sd = 6.7) of pikas residing on natural sites in the north and south, respectively, were estimated to be recent migrants from adjacent artificial sites. No significant migration from natural to artificial sites was detected in either region.  2.4 Discussion  In this study, we investigated genetic variation and connectivity within and among sites occupied by American pikas across a human-modified landscape. We detected evidence that American pikas are influenced by habitat modification at both site- and landscape-level spatial scales, the nature of which may have implications for metapopulation dynamics of wildlife populations inhabiting such landscapes.   Previous research of American pikas inhabiting mining locations have shown that these populations can be highly variable and smaller than expected given habitat availability (Smith 1980; Moilanen et al. 1998), which may have potential genetic consequences. For example, in ore dumps in Bodie, California, American pika averaged 2.49 alleles per locus (Klinger and Peacock, in prep), similar to the levels reported for artificial sites here (AR = 2.33). These values are substantially lower than those reported by studies that used partially  21 overlapping loci within natural habitat in the range core in Nevada (AR = 4.4; Meredith 2002), and Oregon (AR = 5.7; Castillo et al. 2014). It is important to note, however, that our study occurred towards the northern range margin of the American pika (Smith & Weston 1990); theory predicts that levels of within-population genetic diversity declines towards range peripheries (Lesica & Allendorf 1995; Durka 1999; Eckert et al. 2008). As a case in point, levels of allelic richness and heterozygosity detected in the current study were similar to those reported at natural sites at the northern range margin in Tweedsmuir South Provincial Park in British Columbia (AR = 2.8, He = 0.62; Henry et al. 2012). Consequently, we cannot disentangle the relative impacts of fine-scale landscape modification from broader-scale range-wide patterns in interpreting the low levels of within-site genetic variation in and around HVC.  Additionally, we saw no difference in either heterozygosity or allelic richness between pika inhabiting artificial and natural sites possibly owing to the limited sample sizes associated with such a fine-scale assessment or a lack of significant demographic perturbation associated with development. There was, however, a significant increase in relatedness of pika on artificial sites. Artificial sites were originally formed by mining activities (1962 or newer). Given their relatively contemporary origin, these sites were likely colonized much more recently than surrounding natural sites, and are therefore potentially subject to founder effects (Mayr 1963; Nei et al. 1975). Moreover, the artificial sites show some evidence of isolation, exhibiting both detectable levels of genetic divergence from and unidirectional migration towards natural sites, which may have contributed to the elevated levels of relatedness. The evidence of directional migration from artificial to natural sites also has  22 implications for metapopulation dynamics in this system. Peacock (1997) found that dispersal in American pikas is resource dependent, where the primary resource is available habitat, and dispersing individuals generally settle on the first available territory. Moreover, immigration patterns in American pikas are largely a function of local demographic processes of birth rates and habitat saturation (Kreuzer & Huntly 2003). Habitat saturation can be highly variable, but can occur even in artificial habitat in a mine setting (Smith 1980). In this context, the artificial sites studied here may have a lower carrying capacity, spurring directional movement towards more natural settings. As a case in point, preliminary analyses indicate significant differences of both thermal and vegetative characteristics between our artificial and natural sites. Both surface and subsurface temperatures at American pika territories on artificial habitat were significantly more variable than their natural counterparts, and ambient temperatures tended to be higher on artificial habitat (Spilker, unpublished data). Additionally, there were marked differences in the plant (forage) communities between the artificial and natural territories, due in part to the types of species used in the reclamation process. However, nutritional (i.e., nitrogen) composition did not notably differ in the plants appearing in haypiles at the two territory types (Leung, unpublished data). These thermal and vegetative differences could alter habitat quality for American pikas on artificial sites and, in turn, influence local metapopulation dynamics and patterns of gene flow as has been found in other well-studied pika populations (Moilanen et al. 1998). Ongoing ecological assessment of American pikas at our study site could further elucidate metapopulation dynamics in the region and help determine the degree to which variable habitat quality may play a role.   On a broader scale, we resolved extensive genetic structure associated with landscape features. We found evidence for a significant genetic break in this system,  23 corresponding to north and south of the highway, respectively (Figure 2.3). Moreover, the central region of the study system bisected by the highway also constitutes an area of high genetic friction (Figure 2.1). An increase in genetic structure from reduced connectivity is a central prediction of the genetic effects of roads on wildlife (Balkenhol & Waits 2009), and can occur over relatively short timespans (Martínez-Cruz et al. 2007). However, the degree of genetic impact is species- and context-specific as exemplified by the lack of genetic structure detected in the pygmy rabbit (Brachylagus idahoensis), another small bodied lagomorph with a presumed limited dispersal ability (Estes-Zumpf et al. 2010). While this study was conducted over a similar geographic scale and across comparable landscape impediments such as highways, creeks, and reservoirs, the study area contained no mining activity or associated reclaimed habitat.  At a finer level, three genetic units were detected south of the highway, one of which was comprised of a single site, SGS, that formed a unique genetic unit despite close proximity to SGE (550 m). This distance is well within the American pika dispersal capacity reported elsewhere across the range (maximum distances between 2 km-10 km; Hafner and Sullivan 1995; Peacock 1997). Interestingly, SGS is the only site completely surrounded by intense mining activity (Cheryl Blair, personal observation). Although direct mining activities may have contributed to the isolation of SGS, the additional structure detected south of the highway may be the result of differing sources or timing of colonization of newly created habitat with landscape modifications subsequently limiting gene flow. This hypothesis could be tested in the future with broader sampling of potential source populations to the south. There are several natural geographic barriers to gene flow that may account, in part,  24 for the genetic structure observed around the mining site. The highway lies at the bottom of a valley representing an approximate 300 m change in elevation, with a small seasonal creek at the bottom. Previous research indicates that both topographic relief (Henry et al. 2012) and water bodies (Castillo et al. 2014) can significantly inhibit pika movement, making natural geographic boundaries a possible alternative explanation for the north/south genetic division. However, natural geographic barriers would only account for the north-south genetic division and not the degree of genetic structure observed in the south nor the pattern of genetic friction across the landscape since no other natural barriers to gene flow were observed. Additionally, the degree of topographic relief previously shown to inhibit American pika dispersal was far more extreme than anything found around HVC (Henry et al. 2012; Robson et al. 2016). Future research could potentially disentangle the influence of natural and anthropogenic barriers to gene flow by using a larger sampling of the American pika genome and coalescence-based genetic analyses to determine if the development of observed genetic structure was concurrent with human modifications of the landscape.   In summary, we found evidence that landscape modifications have likely influenced the distribution of genetic variation within this study system, documenting several of the expected patterns of fragmentation on small mammals (Gaines et al. 1997). Specifically, we detected site-level changes in genetic characteristics of pika, a slight but significant degree of genetic differentiation of American pikas inhabiting artificial sites, and significant impacts on genetic structuring and migration that were likely associated with landscape modifications. These alterations could influence metapopulation dynamics, including responses to future environmental stressors. In general, it appears that inhabiting artificial habitat might predispose some species to develop fine-scale genetic structure due, in part, to the  25 colonization patterns of the newly available area. Additionally, by its nature, artificial habitat is generally in close proximity to other landscape modifications; in this study, the artificial habitat sites were bisected by a major highway. These additional landscape modifications could act to further reinforce the development of fine-scale genetic structure.  Overall, this area of reclamation appears successful in promoting occupancy for American pikas within HVC, even though the species was not specifically targeted; however, barriers to gene flow likely associated with resource extraction and road construction may limit connectivity across the landscape. Mitigation strategies for promoting connectivity may be limited for American pikas given their thermal sensitivity and habitat requirements. However, American pikas have been documented inhabiting riprap around a small bridge (Henry et al. 2012), indicating habitat corridors and highway bypasses may be effective in this species, but additional study is required. Furthermore, awareness of the potential demographic and genetic consequences of similar landscape alterations may help encourage the integration of mitigation promoting connectivity directly into management planning in order to benefit other wildlife species in the affected areas. Additionally, this study may serve as a reference point for fine-scale genetic analysis across a human-modified landscape enabling contrast between natural and anthropogenically-induced genetic structure in the American pika.   26 Table 2.1 Sample locations, approximate patch size (m2), sample sizes (n), and genetic diversity metrics for the 15 sampling sites. Sex ratios (males per female; M:F), unbiased expected heterozygosity (He), rarified allelic richness (AR), inbreeding coefficient (Fis), and relatedness (rxy) are shown for each site, type indicates either natural talus or artificial habitat. Significance (α< 0.05) is shown by an asterisk for Fis and rxy.   Site Area n M:F Type He AR Fis rxy BBB 2120 6 1 Artificial 0.582 2.21 -0.007 0.363* BLD 9800 6 5 Artificial 0.603 2.32 -0.179 0.394* BRCC 1200 4 3 Natural 0.675 2.48 0.254* 0.072 BRCH 2200 6 0.5 Natural 0.707 2.58 0.018 0.028 BSDE 195 3 2 Artificial 0.618 2.35 0.200 0.111 BSDG 7000 5 1.5 Artificial 0.604 2.27 0.008 0.283* EC 1800 3 2 Natural 0.530 2.20 0.270* 0.210 FHF 10960 13 1.2 Natural 0.654 2.45 0.083 0.143* HFGR 1740 14 1.3 Artificial 0.590 2.24 -0.038 0.248* HNR 1285 6 0.5 Artificial 0.601 2.30 0.021 0.206* Relic_1 8000 16 1.7 Natural 0.596 2.31 0.073 0.209* Relic_2 1800 4 0.3 Natural 0.601 2.34 0.018 0.180* SGE 5800 7 0.8 Artificial 0.638 2.40 -0.025 0.143* SGS 4720 12 3 Artificial 0.557 2.14 0.024 0.262* TG 3135 4 0.3 Natural 0.594 2.31 0.022 0.180*   27 Table 2.2 AMOVA results showing the distribution of genetic variation explained by different hierarchical classification schemes. Numbers represent percentages of genetic variation explained by each grouping. The STRUCTURE clusters (K = 2) represent the north/south divide of highway 97C, while the STRUCTURE clusters (K = 4) takes into account the observed substructure in the south. Natural/artificial (K = 2) groups sites by habitat type, while natural/artificial (K = 4) includes groupings by habitat type divided into north and south of highway 97C.       Among groups Among sites within groups Within sites Significance  STRUCTURE clusters (K = 2) 5.7 10.5 83.8 <0.001 STRUCTURE clusters (K = 4) 6.7 8.7 84.6 <0.001 Natural/artificial (K = 2) 0.9 13.4 85.7 <0.005 Natural/artificial (K = 4) 6.0 9.0 85.0 <0.001     28 Table 2.3. Pairwise site comparisons of genetic differentiation (θ) for American pika within and around Highland Valley Copper and Highway 97C.    BLD BRCC BRCH BSDE BSDG EC FHF TG HFGR HNR Relic_1 Relic_2 SGE SGS BBB 0.234 0.160 0.115 0.117 0.169 0.178 0.163* 0.152 0.253* 0.224 0.210* 0.177 0.165* 0.191* BLD   0.131 0.168 0.116 0.080 0.253 0.134* 0.194 0.200* 0.209 0.228* 0.255 0.215* 0.268* BRCC     0.048 0.132 0.064 -0.005 0.070 0.059 0.118 0.131 0.097 0.107 0.107 0.216* BRCH       0.092 0.109 0.080 0.037 0.080 0.096* 0.053 0.097* 0.110 0.043 0.102* BSDE         0.084 0.174 0.131 0.103 0.147 0.110 0.124 0.099 0.103 0.210 BSDG           0.175 0.121* 0.137 0.161* 0.174 0.195* 0.188 0.163 0.211* EC             0.121 0.127 0.188 0.194 0.149 0.167 0.114 0.217 FHF               0.046 0.142* 0.105* 0.141* 0.156* 0.092* 0.202* TG                 0.186* 0.152 0.181* 0.143 0.152 0.261* HFGR                   0.024 0.108* 0.055 0.092* 0.151* HNR                     0.088* 0.023 0.063 0.094* Relic_1                       0.042 0.017 0.178* Relic_2                         0.039 0.157* SGE                           0.115* * indicates values that are statistically significant after correction for false discovery rate, Pcritical < 0.015  29  Figure 2.1. Map of study area indicating sampling locations within and around the Highland Valley Copper mine in south central British Columbia (a). See Table 2.1 for site descriptions. Green diamonds represent sites with natural habitat while red squares represent sites with artificial habitat. Topographic lines show 100 m changes in elevation and highway 97C is shown. Gridlines show 0.05° changes in latitude and longitude. Two large open pit mines are visible at (-121.04, 50.59) and (-121.04, 50.45). Darker areas represent wooded areas whereas lighter areas represent cleared areas from mining or other human activities; talus is not visible at this scale. Aerial photograph obtained from Google™Earth. Genetic friction map computed across American pika sampling locations (b). The degree of genetic friction is indicated by the inset contour lines and color (red indicates a relative increase in genetic differentiation per unit of geographic distance). Points indicate relative locations of sampling sites as shown on the site map.   −121.00 −120.9650.4450.4650.4850.5050.5250.5450.56V2V30.1207840.1207860.1207880.120790a. b.  30   Figure 2.2. Photograph depicting typical artificial habitat (a) contrasted with natural habitat (b) in the Highland Valley Copper mine. Artificial habitat consisted of rock dumps generated from mining activities or road construction. Relatively undisturbed talus patches were considered natural habitat.    a.   b.   31   Figure 2.3. STRUCTURE bar plots averaged over 25 iterations showing the genetic division (a) of pikas between sites north (BBB through TG) and south of highway 97C (HFGR through SGS; ΔK = 507.9) and (b) in the south only when analyzed independently (ΔK = 73.7). No further genetic subdivisions were resolved within the northern genetic unit. a.   b.  0 0.5 1 BBB BLD BRCC BRCH BSDE BSDG EC FHF TG HFGR HNR RELIC_1 RELIC_2 SGE SGS 0 0.5 1 HFGR HNR RELIC_1 RELIC_2 SGE SGS  32 Chapter 3: Individual-based analysis of hair corticosterone reveals factors influencing chronic stress in the American pika  3.1 Background Rapid environmental change represents a potential stressor and selective force on wildlife populations (Reeder & Kramer 2005; Wingfield & Romero 2011). The main physiological response to long-term environmental stress is the activation of the hypothalamic-pituitary-adrenal axis (HPA) resulting in the release of glucocorticoids (GC), in the form of corticosterone or cortisol, into the bloodstream (Sapolsky et al. 2000; Sheriff et al. 2011). This increase in GC facilitates a suite of adaptive responses to stressful stimuli, such as behavioral changes and energy mobilization via gluconeogenesis, which can enhance short-term survival (Wingfield & Romero 2011). However, long-term activation of the HPA signifies chronic stress and can have detrimental physiological consequences including: suppressed immune response and growth, severe protein loss, fat deposition and hypertension, as well as undesirable behavioral changes including decreased cognitive functioning, inhibition of reproductive behavior, and depression (Wingfield et al. 1998; Sapolsky et al. 2000). For these reasons, the relative levels of GC often reflect overall health and fitness (Blas et al. 2007; Bonier et al. 2009), and measurement of GC is increasingly incorporated into ecological and conservation studies (Busch & Hayward 2009). The relative strengths and feasibility of methodologies to assess stress in wildlife has therefore been a major recent focus (Sheriff et al. 2011).  Several techniques have been developed for measuring stress in wildlife, including measuring GC levels in hair, blood, or saliva, and measuring glucocorticoid metabolites  33 (GCM) in fecal samples. Levels of GC in both saliva and blood respond rapidly to stress and therefore require capture techniques that allow sampling before the activation of the HPA in response to capture stress (generally 2-5 minutes; Sheriff et al. 2011). Where rapid sampling has been possible, this technique has revealed fundamental insights into factors governing GC levels within and across mammalian species. For instance, in a recent meta-analysis, Haase et al. (2016) found a surprisingly strong connection between plasma cortisol levels and mass-specific metabolic rate across a wide variety of taxa, providing a predictive framework for GC levels within and among species. However, obtaining timely blood samples may not be feasible for all species. In such cases, measuring fecal GCM offers a less invasive technique for assessing stress in wildlife (Touma & Palme 2005). Fecal samples accumulate the metabolic byproducts of stress hormones only during gut passage, and therefore primarily reflect chronic stress experienced over a number of hours or days. Additionally, most species exhibit diurnal and seasonal shifts in GC (Reeder & Kramer 2005; Sheriff et al. 2012), making it necessary to obtain a time-series of samples to effectively assess long-term chronic stress using blood, saliva or fecal samples.  The measurement of GC incorporated into hair is a relatively new approach to assess stress in wildlife (Koren et al. 2002). While the direct mechanism by which GC is incorporated into hair is still unknown (Gow et al. 2010), this sample source offers the potential to measure relative levels of GC over the duration of time the hair was grown, which typically encompasses several weeks or months. This longer-term record makes hair analysis a powerful approach for assessing long-term chronic stress (Russell et al. 2012). Despite this major advantage, only a limited number of wildlife studies have used hair as opposed to more established alternatives (Sheriff et al. 2011). While it is likely that using  34 hair samples for stress analysis would provide deeper insights into climate-induced stress in wildlife, another study cautioned that this approach may be more appropriate for detecting population rather than individual stress responses (Mastromonaco et al. 2014). This suggestion came after analyzing long-term trends in fecal GCM and hair GC in eastern chipmunks, Tamias striatus, there was a significant increase in GCM associated with logging but no change in hair GC. The authors concluded the time period measured by hair samples was too long to reflect individual differences in stress. However, this critique would depend on the research objective and species examined. If this methodology can detect individual responses to long-term chronic stress, then it could afford important insights into physiological responses to environmental stressors in climate-sensitive species.   The American pika, Ochotona princeps, is a small lagomorph generally considered to be a thermally sensitive, cold-adapted specialist (Figure 3.1; MacArthur & Wang 1974; Smith 1974). Pikas have an exceptionally high metabolic rate (Lovegrove 2003) and low thermal conductance (MacArthur & Wang 1973), which allows them to survive in an alpine climate without hibernating. However, these features also result in the pika having a resting body temperature only a few degrees below its lethal threshold (MacArthur & Wang 1974; Smith 1974). Pikas require access to cool microclimates to behaviorally thermoregulate (MacArthur & Wang 1974; Hafner 1993). It is thought that this thermal sensitivity may predispose the American pika to the negative ramifications of climate change, casting them as a sentinel species for detecting the ecological consequences of climate change (Beever et al. 2003; Wilkening et al. 2011; Jeffress et al. 2013; Schwalm et al. 2016). Recent analysis of pika populations in the Great Basin supports this view; the minimum elevation inhabited by pikas in the region has risen by 150 m in the past century (Grayson 2005), and climate has  35 been implicated in local extirpations (Beever et al. 2011). Therefore, assessing the biotic response to rapid environmental change in the American pika has become increasingly important as an early warning sign.  Here, we evaluated the utility of hair samples for measuring long-term chronic stress in the American pika. First, we demonstrate the sensitivity of the assay protocol through validation, and we compare estimates of GC from hair with estimates of plasma GC and fecal GCM from the same individuals. Next, we apply this method to investigate relationships between hair GC, microclimate, body size and sex over two elevational gradients to assess whether hair samples can provide direct insights related to climate-mediated stress responses.  3.2 Methods 3.2.1 Sample site and sample collection Sample sites were located along two previously established transects representing elevational gradients in North Cascades National Park, WA (NOCA; Russello et al. 2015). Four sites were located along each transect, spanning ~1,000 m (Figure 3.2, Table 3.1). The two transects, Thornton Lakes (TL) and Pyramid Peak (PP), ran roughly southwest and northeast, respectively. Each transect included the highest (subalpine) and lowest occupied sites identified in its respective region.   Transects were sampled between July 20th and August 29th of 2014; pikas were live-trapped using Tomahawk (Hazelhurst, WI) model 202 collapsible traps following University of British Columbia Animal Care Protocol # A11-0371-006 and U.S. National Park Service Permit # NOCA-2014-SCI-0022. Trapping was generally conducted between 0700 and 1100 or 1600 and 2000 when temperatures were between 5°C and 18°C. After capture, each  36 animal was transferred to a handling bag. A small tuft of fur (about 5 mm2; Figure 3.1 inset) was plucked from a hind limb, as shaving was not practical in this species. Hair samples were stored in individually labeled coin envelopes within a container of silica desiccant. Two small samples of ear tissue (3 mm in diameter) were collected and stored in ethanol for molecular sexing. Cranial diameter was measured with digital calipers to the nearest millimeter, and body mass was measured to the nearest 5 grams using a Pesola scale.   3.2.2 Microclimate measurements Microclimate measurements were taken at each site by deploying four temperature sensors (DS1921G Thermochron i-Button, Maxim Integrated Products, Sunnyvale, CA). Following sampling, sensors were deployed in weather-proof housing; two “ambient sensors” were placed 1.5 m above the talus, each under a white plastic shade in neighboring trees, while two “talus sensors” were deployed in a central location at each site approximately 0.8 m below the talus surface. Temperature was recorded every four hours (starting at 0200) during August 24-31, 2014 and June 1-August 15, 2015. To represent summer microclimatic differences among sites, mean daily maximum and mean daytime (1000 to 1800) temperatures were averaged for the two talus (Tal_max and Tal_day, respectively) and the two ambient sensors at each site (Amb_max and Amb_day, respectively). Additionally, mean nighttime (1800 to 1000) talus temperatures were calculated for each site (Tal_night).  3.2.3 Molecular sexing Morphological differences between male and female genitalia are poorly defined in pika (Duke 1951); therefore, sex was determined using the molecular protocol described by Lamb  37 et al. (2013). DNA was extracted from tissue samples using the Macherey-Nagel NucleoSpin Tissue kit (Macherey-Nagel GmbH & Co. KG, Duren, Germany) following the manufacturer’s protocols. Sex was determined by the selective co-amplification of an allosomal-linked locus (SRY) and an autosomal control locus (Ocp 10). Scoring was conducted by running the PCR product on a 1.5% agarose gel containing 2.5% SYBR Safe (Invitrogen, Carlsbad, California). To ensure accuracy, 50% of the samples were sexed independently a second time and assigned sexes were compared.  3.2.4 Extraction and immunoassay of corticosterone Corticosterone extraction from hair samples followed Meyer et al. (2014) using the DetectX® Corticosterone Enzyme Immunoassay (EIA) kit (Arbor Assays Design, Inc., catalogue no. K014-H1). All hair follicles were removed with a razor blade to avoid the addition of skin tissue (Gow et al. 2010). The remaining hair was added to a 15 mL tube, washed twice with 3 mL 99.7% high-performance liquid chromatography (HPLC) grade isopropanol by rotating for 3 minutes, then decanted to remove external contaminants, and dried under a fume hood. Dried samples were weighed to the nearest 0.1 mg and transferred to a reinforced 2.0 mL tube with three 3.2 mm chrome-steel beads. Samples were then pulverized in 3-minute intervals for 3-18 minutes at 30 hertz on a MM301 Mixer Mill (Retsch®, Newtown, PA). Once samples were uniformly pulverized, 1.5 mL of HPLC-grade methanol was added then samples were rotated for 24 hours at room temperature. Samples were then centrifuged at 12,000 rpm for 10 minutes and 1 mL of the supernatant was transferred to a 1.5 mL microcentrifuge tube without disturbing the hair pellet. This extract was dried under a gentle stream of air in a fume hood for approximately 1-3 days until all methanol had evaporated.  38 The extract was reconstituted using 200 µL of the EIA buffer supplied with the kit, vortexed vigorously, and then immediately frozen at -20°C until analyzed.   Each sample was run in duplicate along with six standard concentrations of corticosterone and two nonspecific binding (NSB) and two maximum binding (BO) wells for each plate. Absorption values at 450 nm were recorded using a Syngery HT microplate reader (Biotek, Winooski, VT). Final concentrations of GC were expressed as picogram (pg) of corticosterone per milligram (mg) of washed, dried hair.  Methods used to measure hair corticosterone concentrations were validated by: 1) demonstrating parallelism between the standard curve and serial dilutions of hair extract; and 2) determining the recoverability of exogenous corticosterone added to hair extracts prior to analysis. For the addition of exogenous corticosterone, six samples were diluted 1:1 each with one of the standard curve solutions. Extraction efficiency likely varies with initial sample mass yielding proportionally higher estimates of glucocorticoids in smaller samples (Millspaugh & Washburn 2004; Tempel & Gutierrez 2004). All corticosterone estimates were plotted against sample mass to identify potential relationships, and, if present, a nonlinear model was fitted using the nls function and used to account for the influence of extraction efficiency. For further comparison, nine additional hair samples were obtained from pikas previously analyzed for plasma corticosterone (Wilkening & Ray 2016) and baseline fecal GCM (i.e., stress levels before capture; Wilkening et al. 2013). Plasma samples were not collected within 3 minutes of capture and thus measured an acute stress response, while fecal samples were collected prior to the stress signature documented for pikas (GCM increases 11-15 hours after capture, Wilkening et al. 2013) and represented a chronic stress response. Relationships between hair, plasma, and fecal GCM were assessed with a linear  39 regression using the lm function.  All analyses (unless otherwise noted) were conducted using the stats package in R version 3.1.3 (R Core Team 2015).   3.2.5 Data analysis The distribution of each independent variable was assessed for normality according to the Shapiro Wilk test using the shapiro.test function, and by inspecting a normal probability plot using the qqnorm function. Any variable that deviated from normality was transformed using a natural log transformation and retested for normality. Since elevation can have an overriding influence on most environmental parameters, we next tested for collinearity among all independent variables by using the Pearson correlation coefficient and corr.test function in R. To reduce collinearity, where significant (alpha = 0.05) and strong (|r| ≥ 0.80) correlations were found, the more physiologically relevant variable was retained based on previous studies of pika ecology (see Results). Mixed-effects models were used to compare corticosterone estimates to all combinations of remaining independent variables using the lmer function of the lme4 package (Bates et al. 2015), after setting REML = FALSE to allow for model selection via AICc. Model fit was assessed by calculating the marginal R2 as suggested by Nakagawa & Schielzeth (2013) using the sem.model.fits function of the piecewiseSEM package (Lefcheck, in press). The top model was tested for all the basic assumptions of linear regression including: linearity, homoscedasticity, and normality of residuals. Additionally, excessively influential data points were assessed using the influence.ME package (Nieuwenhuis et al. 2012).     40 3.3 Results 3.3.1 Laboratory validation Serial dilutions showed a parallel response of samples across the entire standard curve, but the GC concentration from the highest dilution (1:24) was lower than expected relative to the standard curve (Figure 3.3). Due to the lower than expected GC concentrations, all remaining samples were analyzed undiluted. The addition curve showed that 99.8% of exogenous corticosterone was recoverable in the sample matrix. The mean intra-assay variation was 2.01 and 1.88%, and the inter-assay variation was 2.02 and 4.58% for the BO and NSB standards, respectively. There was a negative nonlinear relationship between sample mass and corticosterone concentration, indicating a decrease in extraction efficiency with increasing sample mass (Figure 3.4). To verify this trend, varying initial quantities of hair (1.2, 3.0, 5.2, 10.2, and 18.3 mg) were extracted and analyzed from one sample (PP04T08) resulting in the same nonlinear relationship. A power model was fitted to the data and all subsequent analyses were conducted on the residuals from this model to correct for extraction efficiency. The residuals of this model showed heteroscedasticity, indicating sample masses less than approximately 2.3 mg could lead to inaccurate estimates of corticosterone.     3.3.2 Equivalence of corticosterone estimates from different sample sources The sample masses for the nine hair samples used to contrast hair GC, plasma GC, and fecal GCM were low (5.8 mg ± 4.7 SD) and the mean coefficient of variation (CV) of corticosterone estimates was high (29.85%) between replicates. There was a substantial skew in plasma GC estimates; to facilitate comparison, plasma GC estimates were natural log transformed. There was a non-significant positive relationship between hair and plasma GC  41 (F = 3.57, df = 7, R2 = 0.338, P = 0.101) and no relationship between hair GC and fecal GCM (F = 1.09, df = 7, R2 = 0.134, P = 0.332). Additionally, there was no relationship between plasma GC and fecal GCM (F = 0.830, df = 7, R2 = 0.106, P = 0.393).  3.3.3 Pika stress analysis along elevational gradients  Hair samples were obtained from a total of 49 pikas (23 females and 26 males; Table 3.1). All pikas were unambiguously sexed with no replicate returning a different sex. A mean of 19.4 mg (±7.3 SD) of washed, trimmed hair was obtained from each sample, and the minimum sample weight was 6.5 mg. Resulting corticosterone estimates had a mean CV of 7.19% between replicates. A one-way ANOVA showed significant deviation in corticosterone levels among the sample sites (Figure 3.5).  Only cranial diameter and body mass were non-normally distributed (respectively: W = 0.948, P = 0.032; W = 0.911, P = 0.001). A natural log transformation did not establish normality nor approximate a normal distribution; therefore, non-transformed data were used in subsequent analysis. Due to significant collinearity among all temperature metrics and elevation (Table 3.2), elevation was eliminated in favor of a more direct assessment of microclimate variation. Similarly, mean ambient daily temperature was eliminated in favor of mean maximum daily temperature (Amb_max), which may be a better metric of thermal stress (Beever et al. 2011). All talus temperature metrics were collinear with Amb_max; however, we included Tal_night as this metric represented the mean nighttime temperature pika were likely subjected to, in contrast to Amb_max, which represented the mean maximum daytime temperature. Finally, body mass was eliminated in favor of cranial diameter since body mass is likely to fluctuate on a seasonal basis and cranial diameter was more accurately  42 measured in the field (personal observation).  A total of 13 mixed-effects models were assessed using Amb_max, Tal_night, sex, and cranial diameter as independent variables, including all possible models except those based on the highly collinear variables Amb_max and Tal_night. The top model incorporated all of these variables except Tal_night (Table 3.3), with lower corticosterone estimates at sites with higher ambient temperatures and for larger, male pika. The marginal R2 for this model showed these variables explained about 36.8% of the variation in corticosterone estimates. The residuals of this model were normally distributed (W = 0.9628, P = 0.1231), showed no signs of heteroscedasticity, nor a linear relationship with the fitted value (F = 0.251, df = 47, R2 = 0.005, P = 0.617) or with the main predictor (cranial estimates; F = 1.411E-26, df = 47, R2 = 3.00E-28, P = 1). No excessively influential data points were identified (Cook’s distance < 0.85 for all sites and < 0.5 for all samples).   3.4 Discussion  3.4.1 Laboratory validation While we agree with other authors regarding the power of biological validation of GC measurements (Touma & Palme 2005; Sheriff et al. 2011), the limited feasibility of validation must also be acknowledged for some species. The logistical difficulties of a long-term stress trial in wild animals can lead to mixed results when conducting a biological validation of hair-based stress analyses (Koren et al. 2002; Mastromonaco et al. 2014). For example, Mastromonaco et al. (2014) was only able to resample 12 of the original 23 eastern chipmunks used in a biological validation of hair samples (ACTH challenge), and only 3 of the 5 samples in their experimental group exhibited elevated glucocorticoid levels.  43 Additionally, some species of conservation interest, such as the American pika, would likely exhibit high mortality during a rigorous stress trial (MacArthur & Wang 1973). However, hair-based measurements of GC have been validated in other lagomorph species; for example, Peric et al. (2017) documented a significant increase in cortisol incorporated into hair samples from New Zealand white rabbits after stressful events.   Using individual-based comparisons, we documented a limited connection between plasma- and hair-based estimates and no connection between hair GC and fecal GCM. This lack of correspondence may be attributable to the different time periods over which these sample sources are sensitive. Levels of GC in the bloodstream can be significantly elevated in just a few minutes after a stressful stimulus (Sheriff et al. 2011), and plasma measurements in our study reflect an acute stress response. GCM measurements reflect GC levels on a time scale of several hours or days prior to collection; however, levels of GC in hair represent the accumulation of GC during the relatively long period of hair growth (Yang et al. 1998; Koren et al. 2002). Accordingly, our hair samples likely measured long-term chronic stress following the summer molt, whereas plasma or fecal samples reflected chronic stress experienced during the few hours or days before the time of capture. Additionally, plasma GC, fecal GCM and hair GC measure slightly different hormonal signatures, and other lagomorph studies have documented a lack of correlation among alternative stress metrics (Monclús et al. 2006; Cabezas et al. 2007). For example, only free GCs (those not bound to the carrier protein, corticosterone-binding globulin) circulating in the bloodstream are metabolized by the liver and converted into GCMs; thus, fecal GCM levels mirror free GC levels, but not total GC levels in plasma (Sheriff et al. 2010). The manner in which GCs are incorporated into hair is largely unknown, so questions remain about whether circulating  44 free GC concentration in the blood is proportionately reflected in hair and the influence of confounding factors such as GC contributions from saliva or scents (Sheriff et al. 2011). These temporal and measurement differences are likely responsible for the weak correlations previously observed when hair hormone levels have been compared to those of plasma (Yang et al. 1998) and fecal samples (Mastromonaco et al. 2014).   One of the potential difficulties of using hair is the apparent decrease in extraction efficiency with higher sample masses. Interestingly, this same pattern was reported in fecal samples for both mourning doves, Zenaida macroura (Millspaugh & Washburn 2004) and California spotted owls, Strix occidentalis (Tempel & Gutierrez 2004), reinforcing the need to correct for extraction efficiency. Our approach was to establish a nonlinear relationship to account for this influence. This approach may be preferable when the mass or number of samples is low, as it obviates the need to standardize sample sizes by eliminating smaller samples or truncating larger ones. Of course, this nonlinear relationship suggests that estimates based on low sample masses are less precise (another reason not to standardize samples to the lowest sample mass). Of particular note, the nine hair samples used here to contrast with estimates of fecal GCM and plasma GC were generally low in mass, potentially contributing to the weak relationship observed. We agree with Macbeth et al. (2010) who recommended a minimum sample weight of 5 g when analyzing GC from hair. In our analysis, the relationship between sample mass and GC estimates was approximately linear for samples larger than 5 g and the residuals of our model explaining extraction efficiency were disproportionally high for the low sample masses. We further recommend researchers consider the influence of sample mass in GC extraction efficiencies even for larger samples; our data suggests such effects continue even at higher sample masses.  45  To summarize, we recommend researchers employ a traditional validation technique if possible when applying a novel stress analysis protocol. When traditional validation is not feasible, as in our case, we suggest a cautious approach to cross-validation as each sample source may measure unique temporal and physiological elements of stress. Additionally, we emphasize the importance of considering and correcting for the influence of extraction efficiency. Sample mass was demonstrated here and elsewhere to have a strong influence on estimated corticosterone concentrations. We recommend standardizing sample masses, when possible, above 5 mg; however, the relationship between extraction efficiency and sample mass may be species-specific, so this cutoff may need to be reassessed, especially in species with hair that is coarser. When logistical considerations make standardizing sample mass impractical, we recommend correcting for extraction efficiency by assessing for a nonlinear correlation between sample mass and GC estimates. The resulting equation can then be used to correct for the influence of sample mass.   3.4.2 Field study We demonstrated the utility of hair samples by directly investigating factors influencing long-term chronic stress at the individual level. The sensitivity of this analysis allowed us to evaluate the American pika for several patterns of stress hormone activity well-documented in other mammals. Our results showed that hair GC was mainly influenced by body size, a pattern perhaps mediated by individual differences in mass-specific metabolic rates. In mammals, there is a negative relationship between body mass and GC concentration, underpinned by a relative increase in mass-specific metabolic rate with decreasing body mass (Haase et al. 2016). Since the production but not the degradation of GC is a metabolic  46 function, smaller pikas with higher metabolic rates would be more prone to accumulate GC. Furthermore, smaller animals generally lose heat faster due to their higher surface area to volume ratios, and would need to elevate their baseline metabolism disproportionally to compensate. Additionally, it may be possible that larger pikas would have longer hair, better insulating them from cold stress or influencing the incorporation of GC into the hair. However, an analysis of New Zealand White rabbits found no influence of hair length or body location on GC estimates using a similar protocol (Comin et al. 2012). To our knowledge, this is the first time that a relationship between GC and body size has been reported within a single species; however, pikas may be exceptional given their relatively high metabolic rate (Lovegrove 2003), and further investigation is needed to determine whether this pattern is prevalent within other mammals. Here, we report perhaps the first direct connection between chronic stress and microclimate variation in the American pika, a species with a reputation for narrow thermal tolerance (Smith 1974; Moyer-Horner et al. 2015). Our results further support the potential for the negative effects of chronic cold stress in this species (Beever et al. 2010, 2011; Ray et al. 2012; Jeffress et al. 2013; Schwalm et al. 2016). The increase in GC observed at colder sites could be a function of when our hair samples were grown. The American pika molts twice each year, during summer and fall (Smith & Weston 1990). While we cannot determine the exact time period over which stress was measured, our samples likely captured the GC profile of pikas just after the summer molt, which typically occurs around June to mid-July (Krear 1965). An increase in GC associated with lower ambient temperatures may indicate the necessity of a higher metabolic rate to maintain homeostasis in colder conditions (Lovegrove 2003), particularly during a molt. Being a small alpine mammal that does not  47 hibernate, the American pika may be especially dependent on a fine-tuned metabolic rate, given that smaller animals would be disproportionately affected by low temperatures (Moyer-Horner et al. 2015). As a case in point, Boratyński et al. (2016) found that both the basal metabolic rate and non-shivering thermogenesis in Siberian hamsters, Phodopus sungorus, was highly plastic during the summer months to meet local thermal conditions. If such patterns generalize to the current study, the elevational pattern of GC reported here may represent the metabolic plasticity of pikas to local thermal conditions. We should note that these data do not refute the potential for heat stress in pikas, as the record of GC in our hair samples would likely have been from early summer when the risk of heat stress was minimal. Additionally, it was the mid-elevational sites that had the lowest stress levels along each of the respective transects, a pattern indicating these sites may have been thermally optimal for pikas, with the potential for stress at lower or higher temperatures. The lower stress levels reported here for male pikas match the general pattern observed in most mammalian species (Reeder & Kramer 2005). As both male and female pikas are highly territorial (Smith & Weston 1990), this aspect of behaviour is unlikely to contribute to sex-specific differences in stress. Our samples likely represented the post-breeding period and thus would not capture the increase in stress associated with mating found in other small mammals (Koren et al. 2008). However, our sampling period coincided with gestation and lactation. The costly metabolic demands of rearing offspring may be responsible for elevated stress hormone levels in female mammals (Gittleman & Thompson 1988; Wade & Schneider 1992). Interestingly, female pikas possess a larger adrenal gland than males (Smith & Weston 1990), potentially to meet these physiological demands. However, the relative stress level of each sex may fluctuate seasonally as males and females  48 perform differing tasks, which could decouple acute and long-term stress measurements. For instance, Wilkening et al. (2013) reported higher GCM levels in male pikas, but a longer duration of GCM response to an acute stressor in females.  One of the known limitations of the elevational transect experimental design is the high degree of covariation among microclimate variables typically observed, which can preclude the identification of specific climate influences (Sundqvist et al. 2013). The high degree of covariation within our microclimate estimates was indicative of the overarching influence of elevation on climate within our sample area. As such, our measurements likely represent relative microclimate differences between our sites, independent of time period. Fittingly, our direct measurement of Amb_max was highly related to mean annual temperature at our sites for 2014 (R2 = 0.878, F = 43.2, df = 6, P < 0.001) using downscaled weather station data from the CLIMATEWNA model (Wang et al. 2012). While this addresses our concern over using microclimate measurements taken subsequent to our hair samples and over a short period of time, it does render identifying more specific climate influences challenging with this dataset. This limitation could potentially be addressed by careful selection of additional elevational transects in the future. In conclusion, we suggest a cautionary approach when attempting GC measurements in a species without the ability to validate the methodology. Identifying biologically relevant and well-supported relationships such as GC covariance with body size can assist in the development of novel measurement protocols. In addition, cross-referencing GC metrics among analysis methods may support novel applications in some cases. However, we urge careful consideration of sample type in addressing physiological stress in wildlife, as sample sources vary in the time periods over which they actively measure stress. Finally, we report  49 the only known correlation between directly measured physiological stress and climate variation in the American pika. Our results add to the recent evidence of cold stress in pikas (Beever et al. 2010, 2011; Ray et al. 2012; Jeffress et al. 2013; Schwalm et al. 2016). We suspect that the elevated metabolic rate needed to endure colder ambient conditions as a small bodied, non-hibernating mammal may be responsible for the elevated GC levels reported here. Further research assessing physiological stress in the American pika may assist in conservation and monitoring efforts as we enter a period of rapid environmental change.      50 Table 3.1. Site description of the Pyramid Peak (PP) and Thornton Lakes (TL) sampling transects in North Cascades National Park, WA. Mean daily temperatures for ambient maximum (Amb_max), ambient daytime (Amb_day), talus maximum (Tal_max), talus daytime (Tal_day), and talus nighttime temperatures (Tal_night) are reported in °C along with elevation (m), sample size (n), mean cranial diameter (mm), mean weight (g), and sex ratio (M/F) of pikas at each site.  Site Elevation Amb_max Amb_day Tal_max Tal_day Tal_night n Cranial Weight M/F PP01 480 25.2 22.8 16.6 16.2 15.8 5 24.5 153.0 0.67 PP02 810 23.8 20.3 14.2 13.7 13.7 6 25.6 138.3 1.00 PP03 1327 20.0 17.7 14.7 13.9 12.6 3 23.5 123.3 0.50 PP04 1550 17.8 15.6 13.1 12.4 11.4 9 23.9 152.2 2.00 TL01 504 23.6 21.2 19.2 18.3 17.9 7 27.7 182.1 0.75 TL02 760 23.3 20.6 16.3 15.8 14.6 8 29.2 182.5 1.67 TL03 1409 17.6 15.7 11.0 10.4 10.4 7 26.9 162.6 1.33 TL04 1665 15.4 14.2 12.4 11.6 10.9 4 23.8 131.3 1.00   51 Table 3.2. Correlation matrix of independent variables including pika morphological variables and site microclimate variables (see Table 3.1 for definitions).      Cranial Weight Elevation Amb_max Amb_day Tal_max Tal_day Weight 0.78*  	 	 	 	 	Elevation -0.35 -0.23  	 	 	 	Amb_max 0.29 0.17 -0.97*  	 	 	Amb_day 0.29 0.19 -0.98* 0.99*  	 	Tal_max 0.28 0.24 -0.86* 0.80* 0.84*  	Tal_day 0.29 0.24 -0.88* 0.83* 0.86* 1.00*  Tal_night 0.30 0.24 -0.93* 0.86* 0.89* 0.98* 0.98*   * Significance after sequential Bonferroni correction (p ≤ 0.007) 52 Table 3.3. Information-theoretic analysis of mixed-effects models explaining variation in corticosterone estimates of American pika samples (see Table 3.1 for definitions and sample sizes). All variables demonstrated negative relationships with corticosterone estimates. The negative slope for sex indicates males had lower corticosterone estimates. Site was used as a random effect in all models and the null model included only the random effect.    Model K AICc ∆ AICC Weight Evidence Ratio Marginal R2 - Cranial - Sex - Amb_max 6 264.60 - 0.248 1.00 0.368 - Amb_max - Cranial 5 264.80 0.20 0.225 1.10 0.326 - Cranial 4 265.31 0.71 0.174 1.43 0.208 - Sex - Cranial 5 265.60 1.00 0.151 1.64 0.239 - Tal_night - Cranial 5 266.80 2.20 0.083 3.00 0.261 - Cranial - Sex - Tal_night 6 266.90 2.30 0.079 3.16 0.299 - Amb_max 4 271.01 6.41 0.010 24.64 0.152 - Amb_max - Sex 5 271.10 6.50 0.010 25.73 0.186 Null 3 271.83 7.23 0.007 37.21 0.000 - Sex 4 272.21 7.61 0.006 44.90 0.027 - Tal_night 4 272.61 8.01 0.005 54.85 0.091 - Tal_night - Sex 5 272.80 8.20 0.004 60.20 0.122    53   Figure 3.1. Photograph of an American pika and reference hair sample (inset) weighing approximately 10 mg. Photo curtesy of Andrew Veale.      54   Figure 3.2. Sample sites in North Cascades National Park, Washington, USA. Thornton Lake (TL) and Pyramid Peak (PP) sampling sites shown as circles. Inset map shows approximate location in Washington state. Topographic lines represent 100-m intervals of elevation.   55  Figure 3.3. Top: Parallelism between the standard curve (solid line with circles) and serial dilutions of one sample (squares, no line). Bottom: Addition curve showing a linear  relationship (P < 0.001) between observed and expected GC when samples were mixed 1:1 with standard concentrations of corticosterone from the standard curve.      2040608010050 500 5000%Bound/total boundCorticosterone concentraction (pg/mL)y = 1.009x - 29.53R² = 0.9980400800120016000 400 800 1200 1600Observed concentraction (pg/mL)Expected concentration (pg/mL) 56  Figure 3.4. Relationship between sample mass and estimated corticosterone concentration using NOCA samples (squares), hair samples from paired plasma and fecal samples (triangles), and differing masses from PP04T08 (circles). Inset exponential relationship (solid  line) was developed using all samples. Dashed line shows our suggested 5-mg minimum sample weight cutoff.   0204060801001201401600 5 10 15 20 25 30 35Corticosterone (pg/mg)Sample Mass (mg)y = 118.09x-0.83P < 0.001 57   Figure 3.5. Box and whisker plot showing average corticosterone per site after correcting for extraction efficiency (see Table 3.1 and Figure 3.2 for site descriptions). Boxes represent medians, 25% and 75% quartiles while whiskers extend through 95% interquartile range. A one-way ANOVA showed significant deviation among sites (F = 5.028, df = 7, P < 0.001). Sites are numbered to reflect relative elevation, where 01 = lowest.  PP01T PP02T PP03T PP04T TL01T TL02T TL03T TL04T−10−50510Corrected corticosterone concentration (pg/mg) 58 Chapter 4: Adaptation to climate warming from environmentally-mediated selection may be impeded due to directional gene flow in a thermal-sensitive mammal   4.1 Background Global climate change has profound biological ramifications for nearly all species, many of which have already responded through alterations in morphology, phenology, behavior, abundance, and/or range shifts (Walther et al. 2002; Parmesan & Yohe 2003; Brown et al. 2016; Pacifici et al. 2017). However, in many cases these ecological responses are constrained by species interactions (Van der Putten et al. 2010; Urban et al. 2012), natural limits to dispersal (Pearson & Dawson 2003), or anthropogenic barriers (Opdam & Wascher 2004). As an alternative or complementary outcome, many species have been able to respond to climate change due to phenotypic plasticity (Merilä and Hendry 2014), but this mechanism alone is unlikely to provide a long-term solution when populations are subjected to consistent directional selection (Visser 2008). Adaptation to new climate conditions offers a potential long-term mechanism for impacted species to maintain viability in rapidly changing environments (Savolainen et al. 2013). Therefore, the identification of adaptations to climate and thermal stress has become an important focus in conservation biology (Reusch & Wood 2007; Merilä & Hendry 2014).   Identifying genetic adaptations requires the initial detection of genetic changes in a population over time in response to a selective pressure and the subsequent demonstration that this genetic change leads to increased fitness. Documenting evolutionary responses to climatic pressure in wild populations at this level has presented many methodological difficulties. Although induced thermal tolerance has been documented in laboratory settings  59 (Cavicchi et al. 1995), few studies have attempted to document longitudinal genetic changes in natural populations associated with climate pressure (Merilä 2012). The high degree of phenotypic plasticity in natural populations and difficulty of longitudinal genetic analysis has largely confounded the direct identification of genetic adaptations. Gienapp et al. (2007) concluded that the direct evidence for evolutionary responses to climate change is limited and insufficient to draw inferences on the adaptability of species to climate change. Even given these difficulties, several studies have documented local adaptation in response to climate change, albeit using common-garden experiments to indirectly assess genetic changes (Bradshaw & Holzapfel 2001; Franks et al. 2007). For example, an investigation of pitcher plant mosquitoes (Wyeomyia smithii) revealed genetic shifts in critical photoperiod related to seasonal shifts induced by climate change (Bradshaw & Holzapfel 2001). Additionally, an experiment that contrasted flowering phenology in Brassica rapa before and after a severe drought showed a rapid genetic response to climate-induced stress (Franks et al. 2007). Still, major challenges exist in identifying specific genomic responses to climate change in a broader array of taxa, including highlighting gene regions that may be influenced by natural selection and determining the pace at which impacted species will respond to environmental change.  Population genomics provides a framework for detecting signatures of natural selection associated with local adaptation (Luikart et al. 2003) and for predicting the long-term viability of species impacted by climate change (Visser 2008; Gienapp & Lof 2013; Beever et al. 2016). In particular, genotype-environment association (GEA) analyses provide a powerful approach for differentiating the locus-specific effects of natural selection from genome-wide patterns of evolution caused by processes such as gene flow, genetic drift, and  60 inbreeding (Balding & Nichols 1995; Vitalis et al. 2001; Oleksyk et al. 2010; Vitti et al. 2013; Rellstab et al. 2015). GEA methods seek to detect natural selection by developing functional association between allelic diversity and environmental pressures. For example, Guo et al. (2016) found strong genomic evidence for adaptive divergence in Andrew’s toads (Bufo andrewsi) by demonstrating elevational patterns in genes associated with binding and metabolic processes while controlling for neutral divergence. Such applications of GEA analyses to identify potential climate adaptations are still very limited in mammals. A recent review of phenotypic and evolutionary responses in mammals found no evidence for adaptation in any of the reviewed studies, concluding that it was unclear whether this lack of evidence indicates a lack of adaptability or simply a paucity of data (Boutin & Lane 2014).   The American pika, Ochotona princeps, has been identified as a sentinel mammalian species for the ecological impacts of climate change (Beever et al. 2003; Wilkening & Ray 2016). Recent population extirpations (Grayson 2005; Beever et al. 2011), genetic evidence of range-wide declines (Galbreath et al. 2009), and direct physiological evidence (Waterhouse et al. 2017) all suggest that American pikas are thermally sensitive. Due to this sensitivity, numerous niche models predict precipitous declines in the range of the America pika in response to climate change (Galbreath et al. 2009; Wilkening et al. 2011; Jeffress et al. 2013; Stewart et al. 2015). A draft reference genome and relatively close phylogenetic relationship with a model organism, the European rabbit (Oryctolagus cuniculus; Lanier and Olson 2009), make the American pika a useful mammalian system to assess genomic responses to climate change (Lemay et al. 2013; Henry & Russello 2013; Russello et al. 2015; Robson et al. 2016). For instance, Henry and Russello (2013) identified genomic signatures of environmentally-induced natural selection using anonymous AFLP markers.  61 Similarly, Lemay et al. (2013) identified fixed differences in transcriptome sequences between populations living in different climate zones. Moreover, several elevationally-associated loci were identified in the first proof-of-concept studies to pair non-invasive sampling and genotyping-by-sequencing (Russello et al. 2015). These studies were useful for providing preliminary genomic insights and resources for the American pika. However, the characteristics of the employed markers and limitations in sampling populations or the genome precluded a thorough investigation of climate adaptation in these previous studies.  In this study, we used a space-for-time design and restriction site-associated DNA sequencing (RADseq; Baird et al. 2008; Etter et al. 2011) to genotype tissue samples of American pika in order to: (1) investigate evidence for local adaptation across two elevational transects established in North Cascades National Park, Washington, USA; (2) annotate candidate gene regions exhibiting robust signatures of selection and explore biological implications relative to previously published results; (3) assess elevational patterns of genomic diversity and gene flow to infer both the extent and directionality of pika movement across a sharp climate gradient; and (4) investigate if individual levels of inbreeding have a functional impact on stress hormone levels associated with climate stress.   4.2 Materials and methods  4.2.1 Sample site and sample collection Two previously established elevational transects in North Cascades National Park (NOCA), WA, USA were resampled for this study (Waterhouse et al. 2017; Russello et al. 2015). Four sites were located along each transect, representing an approximately span of 1,000 m in elevation, an expected 6.5°C change in mean annual temperature (Briggs 1997), and an  62 approximate 10% reduction in available oxygen (Peacock 1998). The Thornton Lakes transect (TL) spanned 6.1 km from the low (TL01) to highest site (TL04) and all sites were approximately linear along this south- to southeast-facing transect (Figure 4.1). The Pyramid Peak transect (PP) was 18.8 km east (from TL01 to PP01), spanned 5.4 km, and exhibited a largely northeast aspect; the limited availability of low-elevation pika populations in the area necessitated the low site be slightly offset from the general aspect of this transect. Both transects started at the lowest elevation at which we found pika occupied sites large enough for our sampling requirements and terminated in the subalpine region. Additionally, we sampled pikas at supplemental low-elevation sites where one or several pikas were identified.   Pikas were live-trapped between July 20th and August 29th of 2014 using Tomahawk  model 202 (Hazelhurst, WI) collapsible traps and following University of British Columbia Animal Care Protocol #A11-0371-006 and U.S. National Park Service Permit # NOCA-2014-SCI-0022. Between 15 and 22 traps were set at each main sampling site around signs of pika activity (e.g., hay piles and latrine sites) and approximately 50 m apart to avoid re-trapping the same individual. Traps were set one week in advance and baited with alfalfa cubes every other day to acclimate the animals to the presence of the traps. Trapping was generally conducted between 0700 and 1100 or 1600 and 2000 when temperatures were between 5°C and 18°C and traps were checked every three hours to avoid heat stress. After capture, each pika was transferred to a handling bag and two small samples of ear tissue (each 3 mm diameter) were collected and stored in ethanol for subsequent DNA extraction. Cranial diameter (zygomatic width) was measured with digital calipers to the nearest millimeter and body mass was measured to the nearest 5 grams using a Pesola scale. A small tuft of fur was removed from each individual for corticosterone analysis (Waterhouse et al.  63 2017). We sampled each of the 8 main sites until a target of 8 pikas per site were captured, or until three consecutive trap sessions failed to capture any new pikas.  4.2.2 Climate measurement  We obtained climate data from each site using two methods: first, by using an elevational model to downscale existing climate measurements; and second, by obtaining direct thermal measurements at each site using temperature sensors. The ClimateWNA model uses partially derived elevational models to downscale climate data from PRISM (Daly et al. 2002) and ANUSPLIN (Hutchinson 1989) to a high-resolution raster providing both downscaled direct climate measurements and derived variables of biological relevance (Wang et al. 2012). Long-term averages (1981-2010) for the directly calculated annual (n = 8) and seasonal (n = 20) variables provided by the ClimateWNA model were calculated based on the longitude, latitude, and elevation of each site (Table S4.1).  Since local climatic variations driven by small-scale topographic elements and talus specific properties are not captured by ClimateWNA, we also took direct measurements of ambient and below-talus temperatures using four Thermochron iButton temperature sensors (model DS1921G, Maxim Integrated Products, Sunnyvale, CA) deployed in weather-proof housings at each of the eight sites. Two “ambient sensors” were deployed 1.5 m above the talus, each under a separate white plastic shade in neighboring trees, and two “talus sensors” were placed in a central location at each site approximately 0.8 m below the talus surface. Each sensor collected temperature readings every 4 hours starting at 0200 from September 8, 2014 through August 15, 2015. To better approximate site-wide conditions, temperature data were averaged between paired sensors (i.e., between the 2 ambient and between the 2  64 interstitial sensors) at each site.  To represent the fine-scale annual climate conditions at each site, we calculated the following biologically relevant summary variables  for ambient and talus sensors at each site (Beever et al. 2010, 2011, 2013; Wilkening et al. 2011): number of days with a daily minimum below -10° C, number of days with a daily maximum above 28° C, mean annual temperature, mean summer temperature (June – August), and mean winter temperature (December – February).  We expected a high degree of collinearity among climate variables due to the influence of elevation on climate patterns (Sundqvist et al. 2013). To reduce the dimensionality of the original climate dataset, a principal component analysis (PCA) was conducted using all directly calculated annual and seasonal variables from ClimateWNA (n = 28) and summary variables from iButton measurements (n = 8). We ran a PCA for each transect using prcomp in R v.3.3.2 (R Core Team 2015) on normalized climate variables (mean set to zero and standard deviation set to 1). Subsequently, we used the PCA eigenvectors that collectively accounted for >95% of the variation in the climate data as synthetic climate variables.   4.2.3 DNA extraction and RADseq genotyping DNA was extracted from one tissue sample per individual using the NucleoSpin Tissue Kit (Macherey-Nagel) following the manufacturers’ suggestions with the addition of RNase A (Qiagen). The extracted concentration of genomic DNA was measured using a Quant-iT™ PicoGreen® dsDNA quantification assay kit (Thermo Fisher Scientific) following the manufacturers’ recommendations and using a ViiA 7 real-time PCR machine (Thermo Fisher Scientific).   65 We constructed two RADseq libraries using a protocol described by Baird et al. (2008) as modified in Lemay & Russello (2015). Briefly, 500 ng of DNA was digested with the Sbf 1 restriction enzyme for each individual (New England Biolabs). P1 adaptors were ligated, each containing a unique barcode for each sample in the library. Barcodes were six nucleotides in length and differed by at least two bases (Hohenlohe et al. 2010; Narum et al. 2013). Once pooled, samples were sheared to a mean length of approximately 500 bps using a Bioruptor® (Diagenode) and all fragments between 400-600 base pairs in length were separated using a Pippin PrepTM (Sage Science) by following the manufacturers’ recommendations. Library amplification was carried out following the addition of the P2 adaptor to pooled samples. A second round of size selection was conducted targeting fragments between 450 and 650 bps to account for the addition of the P2 adaptor. Libraries were sequenced using one full lane of Illumina HiSeq 2000 per library and targeting 100-bp single-end reads.   4.2.4 Reference assembly and SNP discovery Libraries were demultiplexed and trimmed to 94-bp using the process_radtags program in STACKS v.1.30 (Catchen et al. 2011, 2013). We used BWA v.0.7.12 (Li & Durbin 2009) to index the pika reference genome (Ensembl, release 88, Ochotona_princeps.pika.dna_sm.toplevel.fa) and masked repeated gene scaffolds in order to make a reference for the subsequent assembly of RADseq reads. Each read was aligned to this reference using the aln algorithm and the resulting SAM files were indexed and converted into BAM files using SAMTOOLS v.1.2 and piped into BCFTOOLS v.1.3.1 for SNP discovery using default parameters (Li et al. 2009a). The resulting VCF file was filtered in  66 VCFTOOLS v.1.12b (Danecek et al. 2011) to ensure: a minimum depth of 6, a maximum depth of 100, no greater than 30% missing data, a minor allele frequency greater than 0.05, only one SNP retained per RADtag, presence of a maximum of 2 alleles, and minimum phred-scale quality score of 20, which corresponds to 99% confidence in genotypes. Additionally, we removed any indels and sequences that mapped to multiple locations of the reference genome. We used the –hardy option in VCFTOOLS to check for Hardy-Weinberg Equilibrium (HWE) at each site with greater than 3 individuals. We excluded any locus that was out of HWE (α = 0.05) at greater than 50% of these sites (minimum of 2). The resulting data were split into PP- and TL-specific datasets, which only contained polymorphic loci present in 70% of the individuals along their respective transects.   4.2.5 Outlier detection and annotation The vast number of genomic markers needed to conduct genomic scans makes controlling for false positives in outlier locus detection a major analytical hurdle (François et al. 2016). One method for overcoming this limitation is to combine multiple independent statistical analyses. In a simulation study, De Villemereuil et al. (2014) found this approach greatly reduces error rates when identifying outlier loci in large genomic datasets. Here, we used three outlier detection methods conducted independently for each transect. First, we used the FST-outlier method of Beaumont and Balding (2004) as implemented in BAYESCAN 2.1 (Foll & Gaggiotti 2008) with a prior odds value of 10, using 100,000 iterations and a burn-in of 50,000 iterations. This analysis was conducted 5 independent times for each transect using a different random seed each time and median scores were calculated for all test statistics. Outliers were identified using a q-value of 0.2 to test for statistical significance.   67  We next used latent factor mixed models (LFMM) to identify specific genomic correlates to environmental parameters using the R package LEA v.1.4.0 (Frichot & François 2015). This method introduces hidden latent factors to account for the influence of genetic structure and spatial autocorrelation of genetic data and then tests for nonrandom associations between environmental parameters and allele frequencies (Frichot et al. 2013). The number of latent factors introduced is determined by the genetic structure of the data. Each transect was first analyzed for the presence of genetic subdivisions (K) using the pca and the snmf functions. We inferred the appropriate number of latent factors by assessing the number of significant axes in the PCA; significance was assessed by comparing eigenvalues to a Tracy-Windom distribution using the tracy.windom function (Patterson et al. 2006). The snmf function uses a Bayesian clustering algorithm to estimate individual admixture coefficients similar to the program STRUCTURE (Pritchard et al. 2000). We identified the K value which returned the minimum cross-entropy in the snmf analysis with 10 repetition for each K value between 1 and 6 for each transect. Next, we used this K value to inform the LFMM analysis assessing if allele frequencies were associated with any of the climate variables. We used the lfmm function with 100,000 iterations, 50,000 burin-in cycles, and 10 replicates. The resulting z-scores were combined using the Stouffer method (Whitlock 2005) recommended by the authors of the software and we applied the genomic inflation factor recommended by Devlin & Roeder (1999). A Bonferroni correction was used to control for multiple comparisons (Rice 1989).       Lastly, we used BAYPASS v.2.1.0 as a third independent method for identifying loci under divergent selection (Gautier 2015). BAYPASS employs an updated version of the Bayesian hierarchical model proposed by Coop et al. (2010). This analysis first constructs a  68 population covariance matrix to estimate neutral genetic divergence, then identifies outliers by checking for nonrandom associations between locus-specific divergence (FST) and environmental parameters using a Bayesian framework. BAYPASS introduces modeling modifications to improve the accuracy of the estimated population covariance matrix and a complete reprogramming of the core Markov Chain Monte Carlo (MCMC) algorithm from BAYENV2 (Günther & Coop 2013), which allows for more efficient analysis of genomic datasets. We conducted five independent runs of the standard covariate model with 20 pilot runs of 1,000 iterations and confirmed the proposal distributions for Metropolis and Metropolis-Hastings updates were within an acceptance range between 0.2 and 0.4 to ensure proper convergence (Gilk et al. 1996). A run length of 25,000 iterations was used, sampling every 25 iterations (thinning procedure) after a burn-in of 5,000 iterations. We selected the run with the value of the lowest deviance information criterion and assessed evidence of natural selection at each locus using a cutoff Bayes Factor of 10 decibans as evidence for ‘strong support’ (following Jeffreys 1961).  We constructed Venn diagrams of overlapping outliers that were significant between transects and methods (BAYSCAN, LFMM, and BAYPASS) using an online tool available at All outliers that were either identified in both transects or in multiple methods were separated into an ‘outlier dataset’. We assessed linkage disequilibrium between all pairs of outlier loci in all populations using the exact test of Guo and Thompson (1992) with 10,000 dememorization steps, 100 batches, and 10,000 iterations per batch implemented in GENEPOP 4.3 (Rayman & Rousset 1995; Rousset 2008) and used a Benjamini-Hochberg correction procedure to assess significance (Benjamini & Hochberg 1995). Additionally, the reference DNA sequence from a 200-bp window around  69 each outlier SNP was subjected to a BLASTN search of all sequences in the NCBI non-redundant database (Altschul et al. 1990). We used a word size of 11, mismatch score of 2, -3, and maximum e-value of 10-10.   4.2.6 Population genetic analysis Natural selection can bias estimates of genetic diversity and divergence (Beaumont & Nichols 1996; Luikart et al. 2003). We form a ‘neutral dataset’ by eliminating all loci that were identified as outliers by any analysis in either transect. This dataset was used to estimate observed and expected heterozygosity, effective number of alleles, and Fis inbreeding (Nei 1978) at each site. We tested for differences in diversity metrics between transects using a two-tailed t-test implemented in R v.3.1.3 (R Core Team 2015). Next, the hierarchical portioning of genomic diversity was estimated by conducting an analysis of molecular variance (AMOVA) and levels of genetic differentiation among sites were assessed using pairwise comparisons of θ (Weir & Cockerham 1984). All genetic diversity tests were performed in GENODIVE v.2.0b27 (Meirmans & Van Tienderen 2004) with a permutation test of 1,000 replicates to test Fis, θ, and AMOVA results for statistical significance.   We investigated the genetic structure among our main sites using a Bayesian clustering method implemented in the program ADMIXTURE v1.3 (Alexander et al. 2009). Overall genetic structure was inferred from the neutral dataset using all the sites in both transects and by running the program for each value of K between 1 and the number of sites plus 2. The optimal number of genetic units was assessed by inspecting the ten-fold cross-validation statistic (CV) and selecting the K with the lowest CV. Strong local genetic structure can mask further subdivision; therefore, we repeated the ADMIXTURE analyses for  70 each transect individually using the same parameters as above. We investigated if outlier loci followed a similar spatial structure by generating transect-specific datasets of loci identified either in both transects or by two outlier detection methods within their respective transects, and reran admixture again using the same parameters as above.   4.2.7 Gene flow analysis Evidence of recent migration (i.e., during the past few generations) along each transect was assessed using the Bayesian method implemented in BAYESASS (Wilson & Rannala 2003). This analysis detects migration by extracting information from transient disequilibria observed at individual multilocus genotypes. Importantly, this method makes relatively few assumptions, allowing it to be applied to nonstationary populations and to genotype proportions out of Hardy-Weinberg equilibrium. We estimated migration rates separately for neutral and outlier loci along each transect using 10,000,000,000 iterations after a burn-in of 1,000,000 steps and sampled every 100th iteration. The inbreeding coefficient, allele frequency, and migration rate mixing parameters were adjusted for each dataset to ensure acceptance rates between 0.2 and 0.6 for each variable. We conducted 5 independent runs each using a separate random number and compared results using the mean coefficient of variance for migration estimates between runs to ensure proper convergence. We used the median migration rates from these runs and constructed 95% credible sets by multiplying the mean standard deviation for each migration rate by 1.96 as suggested by Wilson & Rannala (2003). Significance was indicated when the credible set did not encompass zero.   The supplemental sites we sampled along each transect likely represent recently dispersed pikas since no established population was identified. To further elucidate migration in this system, we conducted an assignment test for each individual at these supplemental  71 sites using 5,000 randomly selected neutral SNPs and the estimator of Rannala & Mountain (1997) in GENECLASS version 2.0.b (Piry et al. 2004).    4.2.8 Stress hormone analysis To investigate the potential role of inbreeding in the climate response of the American pika, we first calculated individual estimates of Fis. To do this, we used a method of moments using the –het function in VCFtools (Danecek et al. 2011). Since this estimate requires population-wide gene frequencies we calculated this estimate separately for each transect using neutral loci (See Admixture results below). We used the hair corticosterone estimates of 49 pikas from Waterhouse et al. (2017) as estimates of long-term chronic stress. This study found that body size, as measured by cranial diameter, sex and maximum summer ambient temperature were all important predictors of individual stress levels. We re-ran the linear mixed-effects models from Waterhouse et al. (2017) by using all combinations of these variables with the addition of individual Fis as fixed effects and used site as a random effect and selected the best mixed effect model via AICc.   4.3 Results 4.3.1 Sample and environmental data collection We sampled 59 pikas during the summer of 2014 (Table 4.1). Our main sites had an average sample size of 6.6 pikas, ranging from 3-11 per site. Additionally, three supplemental sites along the PP transect and two along the TL transect were sampled. One pika was captured from each supplemental site with the exception of TL0.5 where two pikas were captured.     All 32 iButton sensors captured temperature data from all sites from September 8, 2014 through August 15, 2015. Despite the 1.5 m-high placement, one pair of ambient  72 sensors (PP04) appeared to be insulated by snow cover, leading to an overestimation of ambient temperatures at that site between January 15 and May 15, 2015. Despite this additional insulation at PP04, there was high correspondence between climate metrics obtained from ClimateWNA and iButtons (Tables S4.1 and S4.2). A Pearson correlation showed these two methods of assessing climate produced highly comparable estimates of mean annual temperature (r = 0.992, n = 8, p < 0.001).   The PCA of ClimateWNA and iButton variables showed the overriding influence of elevation on site-level climate. The first principal component (PC1) explained 90.3% and 94.9% of the variance in climate in the PP and TL transects, respectively. Adding the second principal component (PC2) increased the amount of variation explained to 98.6% and 99.1% for the PP and TL transects, respectively. PC1 had nearly equal loadings of all the climate variables along each transect (Figure S4.1) and was positively correlated to elevation in both the PP transect (r = 0.999, n = 4, p = 0.001) and TL transect (r = 0.999, n = 4, p = 0.001). Along the PP transect, the second principal component primarily reflected high solar radiation, high below-talus winter temperature and low ambient temperature. Along the TL transect, the second principal component primarily reflected low estimated solar radiation, high temperature differential between summer and winter, high talus temperature, and low ambient temperature.   4.3.2 Reference assembly and SNP discovery DNA sequencing resulted in ~317.8 million reads; after demultiplexing and filtering reads with ambiguous barcodes, low quality scores (phred score < 10), and ambiguous RAD-tags, we obtained a mean of 4.1 (SD = 1.2) million reads per sample. Approximately 76.1% (SD =  73 0.44%) of these reads aligned to the reference genome scaffolding during assembly. After various filtering methods, we obtained a final dataset of 30,763 high-quality SNPs for the 59 samples with an average depth of 25.6 reads (SD = 7.5; Table 4.2). Samples had a mean of 6.4% missing genotypes (SD = 4.8%) per individual. There were 3,634 and 3,784 loci that were monomorphic along the PP and TL transects, respectively.     4.3.3 Outlier detection The BAYESCAN analysis identified 67 and 52 outlier loci along the PP and TL transects, respectively; among these, one locus was shared between the two transects (Table 4.3). Outlier loci had high Fst values, with a mean of 0.226 (SD = 0.052) and 0.225 (SD = 0.044) as compared to putatively neutral loci with mean Fst values of 0.056 (SD = 0.005) and 0.070 (SD = 0.006), along the PP and TL transects, respectively. These results were highly consistent; the mean coefficient of variation for q values between runs was 1.01% and 0.95% for PP and TL transects, respectively.  To conduct the LFMM analysis, we first determined the appropriate number of latent factors by summarizing our genomic data along each transect using a PCA and snmf analysis. The PCA indicated no axes were significant for the PP transect (p > 0.20), whereas the first two axes explained 18.7% of the genomic variation along the TL transect (p < 0.01). The snmf analysis for the PP transect returned a minimum cross-entropy for K = 1 of 0.879 with K = 2 being the next lowest of 0.885. The snmf analysis supported K = 2 for the TL transect returning a cross-entropy of 0.763 while the next lowest was K = 4 with a minimum cross-entropy of 0.765. We therefore used one latent factor along the PP transect and two latent factors along the TL transect to account for neutral genetic structure. The LFMM analysis  74 highlighted 991 loci along the PP transect as being correlated with either PC1 or PC2, with 23 being correlated to both principal components. There were 445 loci along the TL transect correlated with either PC1 or PC2, with 21 being correlated to both principal components. Of these, 54 loci were highlighted independently in both transects by the LFMM analysis (Figure 4.2a and 4.2b).   BAYPASS highlighted 292 loci associated with either PC1 or PC2 along the PP transect, with one locus being associated with both environmental principal components. Along the TL transect, 354 loci were associated with either PC1 or PC2, with one locus being associated with both environmental principal components. There were seven outlier loci that were highlighted independently in both transects by the BAYPASS analysis (Figure 4.2c and 4.2d).   To control for type 1 error, we considered only loci highlighted by multiple methods or in both transects for further analysis. This method indicated robust evidence for 173 outlier loci (Table 4.3). There were 70 loci highlighted by multiple methods along the PP transect and 55 along the TL transect. A total of 79 loci were highlighted in both transects as outliers, 7 of which were highlighted by multiple methods in one of the transects as well. There were no indications of linkage between any pairwise comparisons of outlier loci in any sites after applying a Benjamini-Hochberg correction procedure.   Our BLASTN results returned 44 matches to the NCBI nr database, 11 of which were to genes of known function including those with metabolic processes involving acyl-coenzyme A, oxygen transportation, cellular structure and division, and immune function (Table 4.4). Several gene ontology (GO) terms were common: iron ion binding, oxidoreductase activity, metal ion binding, protein binding, and hydrolase activity were all  75 associated with at least two of these genes (Table S4.3).      4.3.4 Population genetic analysis Population genetic analyses based on the neutral dataset (e.g. outliers removed; n = 28,750 loci; Table 2) revealed higher genetic diversity along the PP transect compared to the TL transect across most genetic diversity metrics (Table 1). Both observed (T = 2.63, df = 5.15, p = 0.05) and expected heterozygosity (T = 5.93, df = 5.90, p < 0.01) were significantly higher along the PP transect than along the TL transect. There tended to be a relatively higher effective number of alleles along the PP transect compared to along the TL transect (T = 2.33, df = 4.76, p = 0.07). Both PP01 and TL03 exhibited significant inbreeding rates, while all other sites showed some indications of outbreeding.     Our results revealed significant genetic divergence between transects and among sites within transects. All θ values between pairwise comparisons of sites between transects were significant, ranging from 0.143 to 0.201 (Table 4.5). Pairwise comparisons within transects showed lower θ values ranging from 0.020 to 0.107, apart from PP01–PP03 and PP02–PP03 all comparisons were significant. Likewise, AMOVA results indicated approximately 11.1% divergence between transects (p = 0.029) and 4.7% divergence among sites within transects (p = 0.001), with the remaining genetic diversity being housed within sites.    The ADMIXTURE results clearly signaled that each transect comprised a unique genetic unit with little genetic mixture between transects (CV = 0.524 for K = 2, next closest CV = 0526 for K = 3; Figure S4.2). We resolved no additional substructure when analyzing the PP transect-specific neutral data and identified weak evidence for two genetic units along the TL transect (Figure 4.3). There was support for further subdivisions when assessing the  76 outlier datasets along each transect. For the PP transect, two genetic units were resolved from outlier loci, while there was evidence to support K = 3 among outliers along the TL transect. In both cases, there was clear evidence of a unique low-elevation genetic unit.   4.3.5 Gene flow analysis All runs of BAYESASS returned optimal acceptance rates after adjusting mixing rates except for the allele frequency parameter (-a), which returned acceptance rates higher than our maximum target of 0.6 in all analyses. Elevated acceptance rates can occur when the likelihood surface is relatively flat, as is likely in biallelic data (Wilson & Rannala 2003). However, independent runs of these analyses returned a low coefficient of variation for all analyses (<1%), indicating consistent model convergence. The only evidence for migration along the PP transect was from PP04 to PP02 and this migration was only apparent when analyzing the neutral dataset (Table 4.6). The TL transect revealed more significant dispersal from TL03 to all other sites along the transect in the neutral dataset and from TL03 to TL02 in the outlier dataset. The 95% credibility intervals of all other migration rates encompassed zero and were therefore taken to be non-significant. The quality index of self-assignment was high along both transects (82.6% and 76.7% long PP and TL, respectively) indicating relatively high assignment power. All three pikas captured from the PP supplemental sites assigned to PP04 with 100% probability. Both pikas from TL0.5 assigned to TL03 while the one pika from TL07 assigned to TL02 again with 100% probability in all cases.   4.3.6 Stress hormone analysis Individual inbreeding levels were relatively low, averaging 0.016 but highly variable among  77 individuals (SD = 0.097). The best fit mixed-effects model incorporated cranial diameter, ambient maximum temperature, and individual Fis (Table 4.7). This model indicated that smaller pikas in colder climates with higher inbreeding levels had elevated stress levels and explained 37.6% of the variance in chronic long-term stress. A similar model but without the individual Fis explained 32.6% of the variance, which indicated the incorporation of Fis increased the explanatory value of the model by approximately 5%.    4.4 Discussion Here we investigated evidence for climate-mediated natural selection by identifying genotype-environment associations (GEA) across two independent elevational transects of American pika populations. We employed RADseq genotyping-by-sequencing within a space-for-time design, which enabled us to circumvent the difficulties associated with long-term genetic monitoring and allowed us to highlight several environmental axes potentially driving natural selection (Lotterhos & Whitlock 2015). Additionally, GEA methods have relatively high power to resolve loci under natural selection (De Mita et al. 2013) and do not require specific candidate loci (Hoffmann & Willi 2008). In our case, these were major advantages given the relative paucity of specific candidate loci from other mammalian systems under similar climatic conditions (Boutin & Lane 2014).   While GEA methods have the power to highlight specific environmental drivers to natural selection, it is common for many biotic and abiotic parameters to be geographically autocorrelated, especially along elevational gradients (Sundqvist et al. 2013). We saw strong evidence for this pattern in our environmental datasets. Both long-term climate measurements from ClimateWNA and contemporary climate measurements taken onsite with temperature sensors provided similar estimates of environmental parameters. To reduce the  78 collinearity of this dataset, we developed synthetic environmental parameters by conducting a PCA of all the climate data along each transect. The result of this analysis highlighted the overriding influence of elevation on climate at our sample sites, limiting the interpretation of genomic correlates to specific environmental variables. For instance, both mean summer and winter temperature decrease with increasing elevation, making it difficult to determine whether selection is occurring in response to heat or cold stress. Another limitation with this approach is an inability to estimate the rate of adaptive evolution precisely. For instance, our environmental gradients covered approximately a 5.5 °C change in mean annual temperature; while we provide evidence of local adaptation over this geographic area, it is unclear over how many generations these adaptations developed. The lack of temporal scale makes it difficult to assess the speed at which species may adapt to future climate conditions, a critical question to assess the continued viability of species in a rapidly changing world (Visser 2008). Even with these caveats, GEA provides a powerful approach for detecting climate-mediated natural selection and can highlight specific genomic targets for further investigation.   4.4.1 Outlier detection It has long been suspected that the thermal sensitivity of the American pika will threaten the species’ long-term viability in the face of climate change, which has already been implicated in local extirpations (Beever et al. 2003; Wilkening et al. 2011; Jeffress et al. 2013). In fact, we recently provided physiological evidence of climate-induced stress in pika populations at our sample sites (Waterhouse et al. 2017). What is less clear is whether these environmental stressors will result in an adaptive response. In general, outlier detection methods can suffer  79 from false-positives (De Villemereuil et al. 2014; François et al. 2016). We addressed this concern by considering outliers to be significant only if they were identified by multiple detection methods or across both elevational transects. This approach may have been overly conservative, potentially missing some significant outliers. However, we felt this was an acceptable tradeoff for the increased robustness to our downstream interpretations, especially as they relate to preliminary evidence for climate-mediated natural selection in the American pika.     Our combined outlier detection method highlighted 173 loci with robust evidence of divergent selection. Close to half of these outliers (45.7%) were identified in both transects and a significant portion were highlighted in both transects and by two independent analyses in one of the transects (12.1%). This congruence was especially meaningful given the strong neutral genetic divergence documented between transects, indicating that PP and TL are effectively independent replicates. This congruent genomic response in two independent transects further suggests that pikas across these elevational transects exhibit a consistent microevolutionary response to climate. The congruent patterns documented here are in contrast to other studies that looked at evolutionary responses to similar environmental change across multiple, independent locations (Muir et al. 2014; Franks et al. 2016). For example, Franks et al. (2016) found only 0.025% of outliers were identified in both of two independent populations of field mustard, Brassica rapa, concluding that these populations were largely on separate evolutionary trajectories for dealing with drought conditions.  The majority of the outliers with robust evidence for natural selection were associated with the first environmental principal component in at least one of the transects (n = 115). This ordination axis was highly correlated to a number of climate variables apparently  80 influenced by elevation. Two of the principal environmental conditions that change with elevation are temperature and oxygen concentration. In general, pika species are inclined to thermal and hypoxic stress (Wang et al. 2006; Yang et al. 2008; Li et al. 2009b). Our analysis indicates that these environmental stressors may be leading to an adaptive response in numerous areas of the American pika genome. These findings complement previous population genetic and transcriptomic studies of pikas along elevational gradients that found evidence for climate-induced natural selection (Lemay et al. 2013; Henry & Russello 2013).          4.4.2 Outlier annotation  A small proportion of our outlier loci returned significant annotations to known genes (6.4%). This low annotation rate is not unexpected, as the challenges related to associating DNA sequence data and gene function in non-model organisms are well-documented (Pop & Salzberg 2008; Yandell & Ence 2012). Nevertheless, several interesting patterns emerged including multiple genes associated with metabolic processes and oxygen transport. For example, scaffold_32008_8033 mapped to a gene encoding for 2-acylglycerol O-acyltransferase 3-like and scaffold_29985_7301 mapped to a gene encoding for acyl-coenzyme A thioesterase 1-like (ACOT1; Table 4.4, Table S4.3). Both of these genes are involved in metabolic processes that use Acyl-CoA, which is known to play a role in maintaining the fluidity of lipid bilayers during cold acclimatization (Nozawa 2011). American pika are known to be thermally sensitive and inclined to cold stress (Beever et al. 2010, 2011; Ray et al. 2012; Jeffress et al. 2013; Schwalm et al. 2016), a trait shared with other pika species. In fact, the Plateau pikas (Ochotona curzoniae) show a number of seasonal regulatory responses to cold including increased activity of cytochrome c oxidase  81 (Wang et al. 2006). Additionally, a study of six pika species revealed evidence for positive selection acting on leptin, an important hormone involved with energy homeostasis, highlighting the role of metabolic processes in pika species to combat colder temperatures (Yang et al. 2008).  Mammalian species are often faced with hypoxic stress from diving, hibernation, and living in burrows or at high elevation (Ramirez et al. 2007). Most pika species live at high elevations and are therefore subjected to selective pressure from hypoxic stress (Hoffman & Smith 2005). Li et al. (2009b) found high expression of HIF-1a in Plateau pikas, which is a key transcription factor involved in a number of cellular and systemic adaptions to hypoxia. This metabolic pathway is mediated by dioxygenase proteins, one of which exhibited putative signature of divergent selection here in the American pika (scaffold_9501_20274). Likewise, another locus aligned with a hemeprotein with known roles in oxygen transport (scaffold_15992_4066; Zaphiropoulos 1997). Taken together, these loci may represent a genomic response in the American pika to living at high elevation brought about by cold and hypoxia stress. These genes and pathways are excellent targets for future validation studies, including elevational analysis of gene expression.   4.4.3 Population genetic analysis There is a consistent relationship between population size and genetic diversity in wildlife populations (Frankham 1996). The American pika is a habitat specialist (Smith & Weston 1990); the higher genetic diversity along the PP transect may indicate more favorable environmental conditions supporting larger populations of pikas. This difference could reflect more favorable thermal conditions arising from the general northern aspect of the PP  82 transect. For instance, both lower sites along the TL transects had higher mean annual temperatures and higher talus summer temperatures, two factors expected to induce thermal stress in American pikas (Varner & Dearing 2014; Stewart et al. 2015). We attempted to calculate effective population sizes from our genomic data to investigate this possibility. However, many of the estimates encompassed infinity, a known issue when using single sample linkage disequilibrium methods to calculate effective population sizes (Luikart et al. 2010). Future studies could test the interaction between population size, genetic diversity, and habitat quality in the American pika.   We resolved significant genetic structure between transects using the neutral dataset and extensive genetic structure among sites within transect using the outlier dataset. Adaptive genetic diversity has a higher power to resolve fine-scale genetic structure due to local adaptation (Funk et al. 2012). Moreover, pikas exhibit limited dispersal (Henry et al. 2012; Castillo et al. 2014; Robson et al. 2016). This pattern of limited dispersal was beneficial to our experimental design demonstrating that the two elevational transects represent largely independent comparisons. However, it remains possible that, despite marked contemporary structure among transects, congruent patterns could be due to ancestral polymorphism. For instance, these shared outliers could potentially represent climate adaptations that occurred lineage-wide during the American pika’s range expansion, preceding the peak of the last glacial maximum approximately 21,000 years ago (Galbreath et al. 2009). A landscape genomic study across the entire North American range of the species is currently in progress that will explicitly test this hypothesis (Russello, pers. com.).     83 4.4.4 Directional migration Dispersal in American pikas is largely resource-dependent, where the primary resource is habitat (Peacock 1997). Additionally, multiple analyses have shown that dispersal capacity is lowest in warm and dry areas that are most physiologically stressful (Castillo et al. 2014, 2016; Schwalm et al. 2016). Smith (1974) reported that at higher elevation, pikas typically occupy a greater proportion of potential territories, a pattern we observed while in the field (Waterhouse, pers. obs.). Our data are consistent with the hypothesis that pika dispersal is habitat-driven; specifically, our gene flow analysis revealed one high elevation site along each transect that acted as the major source of migrants. Along the TL transect, TL03 was a source of migrants to all other sites; while TL03 was not the highest site, we believe this site did represent the highest population density along the TL transect. Along the PP transect, PP04 was the only site with significant emigration and represented the largest population we sampled in NOCA. Additionally, 5 of the 6 pikas sampled at supplemental sites were from the main sources of emigration identified from our directional migration analysis suggesting that even at lower elevations the high-elevation sites are still the major source of immigrants. Overall, these patterns suggest that the upward retreat of American pika populations with climate change is likely to occur through the extirpation of low elevational populations rather than their movement upslope (Beever et al. 2003; Grayson 2005). Likewise, we observed relatively minimal gene flow in the outlier dataset. One of the strengths of the migration analysis employed is that it does not require loci to conform to Hardy-Weinberg expectations and can be applied to nonstationary populations (Wilson & Rannala 2003). However, since this analysis is based on detecting disequilibrium in individual multilocus genotypes, it is sensitive to the number of genetic markers used; the  84 lack of migration we documented in the outlier dataset could therefore have been a function of the lower number of markers (transect-specific datasets: outlier n = 123-129; neutral n = 25,211-25,297). Alternatively, the lower level of gene flow could be biologically meaningful,  representing resistance to genetic introgression at locally adapted loci (Tigano & Friesen 2016).  The reduced levels of gene flow at putatively adaptive loci and the general downward direction of migration indicate potential impediments for the movement of thermal adaptations across elevational gradients of pika populations. Adaptations resulting from higher temperatures occurring at lower sites may not be able to move upward as the climate warms. This concern may be present in other alpine species as well, as lower elevation populations are expected to suffer reduced viability while higher elevation populations may persist at high densities (Parmesan & Yohe 2003). Under this scenario, we find it likely that source-sink dynamics in American pikas and other alpine species would naturally exhibit a down-slope bias, where high elevation populations act as sources of migrants to generally smaller, lower elevation populations.   4.4.5 Stress hormone analysis Inbreeding is a common occurrence in wild populations and can have deleterious effects through inbreeding depression (Keller & Waller 2002). If sever enough, these impacts can threaten the continued viability of the species (O’Grady et al. 2006). The American pika is known to exhibit mate choice based on intermediate levels of relatedness (Peacock & Smith 1997b) and generally dispersal is dictated by competition for territory rather than inbreeding avoidance (Peacock 1997). We investigated the potential consequences of inbreeding in pika  85 by estimating individual levels of inbreeding and assessing if these levels had an impact on long-term chronic stress as measured by extracted corticosterone from hair samples. We documented a positive correlation between inbreeding and corticosterone levels, which indicates the potential for inbreeding depression in this species especially as it relates to climate stress. Climate change is expected to cause continued reductions in populations sizes of American pikas (Beever et al. 2013; Yandow et al. 2015) which could further exacerbate inbreeding trends (Frankham 1995). With this in mind, it may be important to continue to monitor inbreeding trends in smaller pika populations.      4.4.6 Summary  In this study, we assessed the potential for climate adaptation in a sentinel mammalian species, the American pika, by identifying genotype-environment associations in populations occurring over an environmental gradient. We resolved robust evidence of natural selection in this system and identified several gene regions potentially associated with local adaptation. Our analysis adds to the growing body of evidence for climate-induced natural selection in the American pika (Lemay et al. 2013; Henry and Russello 2013) and provides a relatively rare mammalian example of genetic adaptation to contemporary climate conditions (Boutin & Lane 2014). Future work could test the potential role of the genes identified here as outliers through broader-scale sampling or analyses of gene expression along elevational gradients. Investigation of the American pika transcriptome has already been completed (Lemay et al. 2013). Additionally, we resolved evidence for consistent directional migration, which may hold important ramifications for the movement of low elevation thermal adaptations. More broadly, the results from this study provide insights useful in applying  86 assisted gene flow for mitigating the negative effects of climate change (Aitken & Whitlock 2013). Continued research is needed to determine if the rate of adaptation in the American pika will be able to keep pace with a rapidly changing environment. If adaptation cannot keep pace with climate change, future conservation efforts could consider the translocation of American pika from populations highlighted as having thermal adaptations to populations suspected of undergoing reduced viability due to thermal stress. Given their thermal sensitivity and the fact that their rocky talus slope habitat is largely not impacted by direct human activities, the American pika may represent a rare mammalian model system for evaluating and implementing such conservation strategies for mitigating the deleterious consequences of contemporary climate change.    87 Table 4.1. Site characteristics of American pika populations sampled in the North Cascades National Park, WA. Observed heterozygosity (Ho), total corrected heterozygosity (He), effective number of alleles (Ae), and inbreeding estimates (Fis) shown for each site with a minimum sample size (n) of 2 along the Pyramid Peak (PP) and Thornton Lakes (TL) in North Cascades National Park, WA, USA.  Site Elevation Northing Easting n Ho He Ae Fis PPTrl 330 636449 5396874 1 - - - - PPHwy 415 640203 5395358 1 - - - - PP01 450 640367 5393862 5 0.27 0.28 1.426 0.038* PP1.5 615 637260 5396080 1 - - - - PP02 820 638227 5395281 6 0.316 0.292 1.459 -0.08 PP03 1330 638082 5393892 3 0.305 0.288 1.415 -0.061 PP04 1580 637027 5392854 9 0.302 0.298 1.475 -0.014 TL0.5 150 625978 5391339 2 0.254 0.264 1.345 0.039* TL0.7 265 624461 5389771 1 - - - - TL01 490 622400 5388108 7 0.286 0.255 1.406 -0.12 TL02 780 622818 5390984 11 0.26 0.26 1.415 -0.003 TL03 1390 623134 5393254 8 0.258 0.269 1.42 0.039* TL04 1700 623614 5393933 4 0.265 0.255 1.384 -0.039 *significant at the 0.001 level   88 Table 4.2. Number of SNPs retained in the American pika genomic dataset after various filtering procedures. Transect specific SNP counts shown for the Pyramid Peak (PP) and Thornton Lakes (TL) transects in the North Cascades National Park, WA, USA (NOCA).     From reads to SNPs SNP count Aligned to reference (min depth 6X) 607,468 Max missing data of 30% 189,696 Minor allele frequency > 0.05 53,162 One SNP/tag, indels and duplicates removed 31,298 Max of 2 alleles and 100x depth, passed HWE, minimum phred quality of 20 30,763 PP max missing data of 30% 30,128 PP Polymorphic  27,129 TL max missing data of 30% 30,411 TL Polymorphic 26,979 Outliers (any technique) 2,013 NOCA Neutral dataset 28,750 PP outlier dataset 123 PP neutral dataset 25,297 TL outlier dataset 129 TL neutral dataset 25,211     89 Table 4.3. Summary of outlier detection analyses in American pika. Number of putative outliers is summarized for each analysis; PC1 and PC2 refer to analyses conducted using the first and second environmental principal component, respectively. Numbers in parentheses correspond to outliers significantly correlated to both principal components.   Method   Analysis Number of loci Percent of loci BAYESCAN    PP 67 0.25   TL 52 0.19   Both  0 0.00 LFMM   PP_PC1 493 1.82   PP_PC2 521 (23) 1.92   TL_PC1 305 1.13   TL_PC2 161 (21) 0.60   Both  54 0.23 BAYPASS   PP_PC1 169 0.62   PP_PC2 124 (1) 0.46   TL_PC1 180 0.67   TL_PC2 175 (1) 0.65   Both  7 0.03 Any   Both  79 0.34  90 Table 4.4. Summary of top BLASTN search results for outlier loci in American pika detected by BAYESCAN (BS), LFMM (LF), and BAYPASS (BP). The number after the analysis in the comparison column indicates which environmental principal components was used in the correlative analysis.   Locus Comparison SNP Top blast hit (accession) Abbreviated description scaffold_34750_45136 PP_BP_2, PP_BS, TL_LF_2 A/G XM_012930919.1 Ochotona princeps coiled-coil domain containing 77 (CCDC77) scaffold_32008_8033 PP_BP_2, TL_BP_1, TL_BS C/T XM_004587187.1 Ochotona princeps 2-acylglycerol O-acyltransferase 3-like scaffold_15992_4066 PP_BP_2, PP_BS, PP_LF_2 T/A XM_012925448.1 Ochotona princeps cytochrome P450 2C18-like scaffold_332_200703 PP_LF_2, TL_LF_1 G/T LT160000.1 Macaca fascicularis complete genome, chromosome chr1 scaffold_35953_8025 PP_LF_2, TL_LF_1 C/T XM_004598126.1 Ochotona princeps MAD1 mitotic arrest deficient-like 1 scaffold_9501_20274 PP_LF_1, TL_LF_1 G/T XM_004595405.1 Ochotona princeps 2-oxoglutarate and iron-dependent oxygenase domain containing 2 scaffold_34209_3344 PP_LF_2, TL_LF_1 C/A XM_002711101.3 Oryctolagus cuniculus tubulin alpha-1A chain scaffold_29985_7301 PP_LF_2, TL_LF_1 G/A XM_004584411.2 Ochotona princeps acyl-coenzyme A thioesterase 1-like (ACOT1) scaffold_3620_102829 TL_BP_2, TL_BS T/C NG_008098.1 TNF receptor superfamily member 11a (TNFRSF11A) scaffold_30_208287 TL_BP_2, TL_BS A/G AL391986.12 Human DNA sequence from clone RP11-426E5 on chromosome 10 scaffold_590_131441 PP_BP_1, PP_BS C/A XM_004578932.1 Ochotona princeps interleukin 20 (IL20) scaffold_168552_899 PP_BP_1, PP_BS C/T XM_004597033.1 Ochotona princeps DEAH (Asp-Glu-Ala-His) box polypeptide 34 (DHX34) scaffold_4242_61990 TL_BP_2, TL_LF_1 A/T XM_012755825.2 Microcebus murinus adenosine monophosphate deaminase 3 (AMPD3) transcript variant X6 scaffold_7931_46314 PP_LF_2, TL_LF_2 G/A CP011890.1 Ovis canadensis canadensis isolate 43U chromosome 5 sequence * All e-scores < 10-10  91 Table 4.5. Pairwise comparisons of Fst (Top diagonal; Weir & Cockerham 1984) between pika populations and associated p-value from a permutation analysis (1,000 permutations; bottom diagonal). Pairwise Fst values from comparisons between transects shaded in grey. See Table 4.1 for site descriptions.     PP01 PP02 PP03 PP04 TL01 TL02 TL03 TL04 PP01 -- 0.072* 0.066 0.050* 0.197* 0.201* 0.173* 0.188* PP02 0.002 -- 0.049 0.043* 0.181* 0.186* 0.158* 0.175* PP03 0.058 0.080 -- 0.036* 0.189* 0.185* 0.154* 0.179* PP04 0.004 0.004 0.044 -- 0.163* 0.172* 0.143* 0.156* TL01 0.003 0.001 0.008 0.001 -- 0.104* 0.075* 0.107* TL02 0.001 0.001 0.003 0.001 0.001 -- 0.031* 0.061* TL03 0.001 0.001 0.004 0.001 0.001 0.004 -- 0.020* TL04 0.014 0.006 0.032 0.004 0.003 0.004 0.032 -- *significant at the 0.05 level  92 Table 4.6. Directional migration analysis of American pika along the Pyramid Peak (PP) and Thornton Lakes (TL) transects using either the neutral or outlier datasets. Migration values are shown as proportion of the population expected to be recent migrants from the source population, value of the 95% credibility interval shown in parentheses, significant migration values shaded in grey (i.e., credibility interval does not encompass zero).     Dataset    Direction Migration Direction Migration Direction Migration Direction Migration PP Outlier   PP01àPP01 0.854 (0.104) PP02àPP01 0.072 (0.087) PP03àPP01 0.037 (0.065) PP04àPP01 0.037 (0.065)   PP01àPP02 0.034 (0.060) PP02àPP02 0.889 (0.097) PP03àPP02 0.033 (0.059) PP04àPP02 0.044 (0.073)   PP01àPP03 0.048 (0.081) PP02àPP03 0.048 (0.082) PP03àPP03 0.857 (0.115) PP04àPP03 0.048 (0.081)   PP01àPP04 0.026 (0.047) PP02àPP04 0.026 (0.046) PP03àPP04 0.026 (0.047) PP04àPP04 0.923 (0.074)     PP01àPP01 0.853 (0.103) PP02àPP01 0.037 (0.065) PP03àPP01 0.037 (0.065) PP04àPP01 0.073 (0.086) PP Neutral   PP01àPP02 0.033 (0.059) PP02àPP02 0.833 (0.098) PP03àPP02 0.033 (0.059) PP04àPP02 0.100 (0.090)     PP01àPP03 0.048 (0.081) PP02àPP03 0.048 (0.081) PP03àPP03 0.809 (0.114) PP04àPP03 0.095 (0.104)     PP01àPP04 0.026 (0.047) PP02àPP04 0.026 (0.046) PP03àPP04 0.026 (0.047) PP04àPP04 0.923 (0.074) TL Outlier   TL01àTL01 0.879 (0.091) TL02àTL01 0.030 (0.054) TL03àTL01 0.061 (0.073) TL04àTL01 0.030 (0.054)   TL01àTL02 0.022 (0.041) TL02àTL02 0.889 (0.077) TL03àTL02 0.067 (0.066) TL04àTL02 0.022 (0.041)   TL01àTL03 0.028 (0.050) TL02àTL03 0.056 (0.067) TL03àTL03 0.860 (0.091) TL04àTL03 0.057 (0.070)   TL01àTL04 0.042 (0.072) TL02àTL04 0.042 (0.072) TL03àTL04 0.049 (0.083) TL04àTL04 0.868 (0.111) TL Neutral   TL01àTL01 0.850 (0.095) TL02àTL01 0.030 (0.054) TL03àTL01 0.090 (0.085) TL04àTL01 0.030 (0.054)   TL01àTL02 0.022 (0.041) TL02àTL02 0.867 (0.080) TL03àTL02 0.089 (0.072) TL04àTL02 0.022 (0.041)   TL01àTL03 0.028 (0.050) TL02àTL03 0.056 (0.067) TL03àTL03 0.889 (0.085) TL04àTL03 0.028 (0.050)   TL01àTL04 0.042 (0.072) TL02àTL04 0.042 (0.072) TL03àTL04 0.123 (0.106) TL04àTL04 0.794 (0.107)   93 Table 4.7. Mixed-effects models explaining chronic long-term stress as measured by hair corticosterone in American pikas sampled in North Cascades National Park, WA, USA. Fixed effects included cranial diameter (cranial), sex, maximum ambient summer temperature (Amb_max) and individuals estimates of inbreeding (Fis). All variables loaded negatively with the excpetion of Fis which loaded positively in each model.   Variable Parameters AICc ∆AICc Weight Evidence Ratio Marginal R2 Fis, Cranial, Amb_max 6 263.35 - 0.278 1.00 0.376 Fis, Cranial 5 263.63 0.28 0.242 1.15 0.269 Fis, Cranial, Sex, Amb_max 7 264.13 0.78 0.188 1.48 0.404 Cranial, Sex, Amb_max 6 264.59 1.24 0.150 1.86 0.368 Fis, Cranial, sex 6 264.80 1.45 0.135 2.06 0.288 Fis, Amb_max 5 270.49 7.14 0.008 35.52 0.199 Fis 4 271.11 7.77 0.006 48.42 0.048 Fis, Sex, Amb_max 6 271.39 8.04 0.005 55.70 0.221 Fis, Sex 5 272.25 8.90 0.003 85.63 0.064  94   Figure 4.1. Map of Pyramid Peak (PP) and Thornton Lakes (TL) transects in the North Cascades National Park, Washington, USA. Topographic lines represent 100 m change in elevation. Highway 20 is indicated by the grey line.  95                        Figure 4.2. Summary plots for genotype-environment associations in American pika. LFMM analysis for the Pyramid Peak (PP) transect (panel A) and Thornton Lakes (TL) transect (panel B); note that p-values are shown on an inverted log10 axis. Panels C and D show the regression coefficient (Beta) and Bayes Factor (BF) in decibans for the BAYPASS analysis for PP and TL transects, respectively. Analyses conducted with the first environmental principal component (circles) and second (squares) are shown for each analysis. Significant outliers identified in both transects are shown in solid blue. Locus ID is shown for all annotated outliers, blue leader lines indicate significance in both transects (See Table 4.4).   96  Figure 4.3. Admixture proportions for American pika along Pyramid Peak (panels A and C) and Thornton Lakes (panels B and D) transects using only neutral (panels A and B) or outlier loci (panels C and D) with inset cross-entropy (CE) estimations for each value of K. Plots show admixture proportions with the lowest CE in each case.| PP01 | PP02 | PP03 | PP0400.250.50.751| TL01 | TL02 | TL03 | TL0400.250.50.751| TL01 | TL02 | TL03 | TL04Neutral0.50.7511 6CEKA) B)| PP01 | PP02 | PP03 | PP04Outliers0.7511.251.51 6CEK0.50.7511.251 6CEKC) D)0.511.51 6CEKPP01 PP02 PP03 PP04PP01 PP02 PP03 PP04TL TL0 TL03 TL04TL TL TL03 TL04 97 Table S4.1. Climate variables from iButton data obtained along the elevational transects established in the North Cascades National Park, WA, USA. See Table 4.1 and Figure 4.1 for site locations and Figure S4.1 for variable definitions. All temperatures are in degrees Celsius.   Site Amb_10 Amb_25 Amb_MAT Amb_MST Amb_MWT Tal_MAT Tal_MST Tal_MWT PP1T 0 29 10.03 19.76 2.76 7.57 16 0.79 PP2T 3 22 8.82 18.62 1.65 6.12 13.74 0.25 PP3T 6 3 6.14 15.39 -0.14 4.48 13.24 -1.18 PP4T 5 0 4.58 13.8 -0.84 3.86 11.83 0.02 TL1T 0 16 10.58 19.25 4.19 9.96 18.13 3.93 TL2T 3 22 8.95 17.93 2.64 7.33 15.16 1.57 TL3T 4 1 5.83 14.06 0.32 3.69 10.46 -0.45 TL4T 8 0 4.21 12.46 -1.14 3.67 11.17 -0.14    98 Table S4.2. ClimateWNA data for each site along the elevational transects established in the North Cascades National Park, WA. See Table 4.1 and Figure 4.1 for site locations and Figure S4.1 for variable definitions. All temperatures are in degrees Celsius, precipitation is in mm, and radiation is in MJ m-2 d-1.  Site MAT  MWMT  MCMT  TD  MAP  MSP  AHM  SHM  Tmax_wt  Tmax_sp  Tmax_sm  Tmax_at  Tmin_wt  Tmin_sp  PP1T 8.2 18.1 -0.9 19.0 1749 265 10.4 68.4 3.0 12.9 23.7 13.1 -3.0 1.9 PP2T 6.7 16.4 -1.9 18.3 2125 324 7.9 50.5 2.0 11.0 21.6 11.5 -4.0 0.5 PP3T 4.2 13.6 -3.9 17.5 2382 369 6.0 36.8 0.2 8.0 18.4 8.9 -6.5 -2.2 PP4T 3.2 12.4 -4.6 17.1 2621 417 5.0 29.8 -0.3 7.0 17.1 8.0 -7.4 -3.2 TL1T 8.7 17.5 0.8 16.7 2003 343 9.3 51.0 4.7 12.8 22.3 13.3 -1.7 2.7 TL2T 7.8 16.6 0.1 16.5 2099 366 8.5 45.3 4.2 11.9 21.4 12.5 -2.4 1.6 TL3T 4.4 13.5 -3.3 16.8 2602 464 5.5 29.2 1.1 7.6 17.7 9.3 -5.9 -1.9 TL4T 2.8 12.1 -4.8 16.9 2819 501 4.6 24.2 -0.3 5.5 15.8 7.7 -7.5 -3.4   Table S4.2. continued   Site Tmin_sm  Tmin_at  Tave_wt  Tave_sp  Tave_sm  Tave_at  PPT_wt  PPT_sp  PPT_sm  PPT_at  Rad_wt  Rad_sp  Rad_sm  Rad_at  PP1T 10.3 3.7 0.0 7.4 17.0 8.4 708 336 125 580 3.8 14.1 20.0 6.8 PP2T 8.7 2.7 -1.0 5.7 15.2 7.1 857 410 152 706 3.3 14.4 20.5 6.7 PP3T 6.2 0.4 -3.1 2.9 12.3 4.7 957 466 173 787 3.2 14.3 20.4 6.4 PP4T 5.2 -0.5 -3.9 1.9 11.1 3.8 1046 521 196 859 4.2 15.2 21.4 7.6 TL1T 10.4 4.8 1.5 7.7 16.4 9.1 796 407 162 639 3.8 15.0 21.0 7.8 TL2T 9.4 4.0 0.9 6.7 15.4 8.2 832 427 172 667 5.7 16.5 22.3 10.7 TL3T 6.4 1.0 -2.4 2.8 12.0 5.1 1012 530 222 838 7.1 17.2 23.4 11.7 TL4T 5.2 -0.3 -3.9 1.0 10.5 3.7 1087 575 240 918 7.4 17.3 23.6 11.6     99 Table S4.3. Gene Ontology terms associated with significant outlier loci Blast hits.   Gene   GO terms Coiled-coil domain containing 77 (CCDC77)   Centrosome   Membrane 2-acylglycerol O-acyltransferase 3-like   2-acylglycerol O-acyltransferase activity   Diacylglycerol O-acyltransferase activity   Transferase activity   Transferase activity, transferring acyl groups   Transferase activity, transferring acyl groups other than amino-acyl      groups  Cytochrome P450 2C18-like   Monooxygenase activity   Iron ion binding   Arachidonic acid epoxygenase activity   Steroid hydroxylase activity   Oxidoreductase activity MAD1 mitotic arrest deficient-like 1   Colocalizes with kinetochore   Condensed chromosome kinetochore   Spindle pole   Nucleus   Colocalizes with nuclear pore 2-oxoglutarate and iron-dependent oxygenase domain containing 2   Iron ion binding   Oxidoreductase activity   Oxidoreductase activity, acting on paired donors, with incorporation or      reduction of molecular oxygen   L-ascorbic acid binding   Metal ion binding Tubulin alpha-1A chain   Nucleotide binding   GTPase activity   Structural molecule activity   Structural constituent of cytoskeleton   Protein binding Acyl-coenzyme A thioesterase 1-like   Very long-chain fatty acid metabolic process   Long-chain fatty acid metabolic process   Acyl-CoA metabolic process   Long-chain fatty-acyl-CoA biosynthetic process TNF receptor superfamily member 11a (TNFRSF11A)   Receptor activity   Transmembrane signaling receptor activity   Tumor necrosis factor-activated receptor activity   Protein binding   Cytokine binding Interleukin 20 (IL20)   Cytokine activity   Interleukin-20 receptor binding   Interleukin-22 receptor binding DEAH (Asp-Glu-Ala-His) box polypeptide 34 (DHX34)   RNA binding   ATP-dependent RNA helicase activity   Helicase activity   ATP binding   Hydrolase activity Monophosphate deaminase 3 (AMPD3), transcript variant X6, mRNA   AMP deaminase activity   Hydrolase activity   Deaminase activity   Metal ion binding  100  Figure S4.1. Climate variables and loadings for the first (PC1) and second principal component (PC2) of the environmental PCA for the Pyramid Peak and Thornton Lakes transects. The first eight variables were directly measured using temperature sensors, and all other variables were estimated using the ClimateWNA model. Bars represent the relative magnitude (negative to the left and positive to the right of the dashed line) of loading for each variable. All temperatures are in degrees Celsius, precipitation is in mm, and radiation is in MJ m-2 d-1.     101           Figure S4.2. Admixture plot for all the sites sampled in NOCA (see Table 4.1 and Figure 4.1 for site locations). Inset cross-entropy (CE) estimations for each value of K show in the inset graph.| PP01 | PP02 |PP03 | PP04 | TL01 | TL02 | TL03 | TL040.50.650.81 10CEKPP01 P 02 PP03 PP04 TL TL TL03 TL04 102 Chapter 5: Conclusion  5.1 Research findings This thesis describes three independent approaches to assess biological responses in the American pika to rapid environmental change. The results of these chapters taken together provide consistent evidence for several biological responses occurring in pikas. Here, I will summarize how these studies add to our understanding of dispersal patterns, thermal stress, and metabolic patterns occurring in the American pika. These patterns were apparent in multiple chapters of this thesis and were delineated by using independent genetic and physiological evidence. Where appropriate, I will outline potential ramifications of each pattern for the continued viability of the species and explore management implications for the American pika and other climate-sensitive species.     5.1.1 Evidence for restricted dispersal  The limited dispersal capabilities of the American pika are well documented. Population genetic analyses have demonstrated restricted gene flow (Henry et al. 2012; Robson et al. 2015; Castillo et al. 2014), while demographic analyses have shown minimal dispersal, typically in the form of juveniles settling in the next available territory (Smith & Ivins 1983). These dispersal characteristics of American pikas have important implications to the metapopulation dynamics of the species (Peacock & Smith 1997a; Moilanen et al. 1998) and for the movement of potentially beneficial climate adaptations (Hoffmann & Willi 2008; Sgrò et al. 2011).   Human modifications to the landscape have had genetic and demographic impacts in numerous species and are a leading cause of biodiversity loss (Fischer & Lindenmayer 2007).  103 However, American pikas typically inhabit high alpine and relatively undisturbed habitat [although see Peacock & Smith (1997)]. In Chapter 2, I investigated the genetic connectivity of pika populations in a human-modified landscape; this approach allowed me to assess the impacts of road construction, which restricted gene flow in this species. Given the ubiquity of roads on the landscape, this finding has potential implications for range shifts associated with climate change (Opdam & Wascher 2004). Specifically, in this study, I found evidence that a major roadway severed gene flow between southern and northern pika populations. If the southern sites become uninhabitable due to climate stress, then northern migration becomes untenable if this roadway is not permeable to dispersal. On a broader scale, such restrictions to dispersal could inhibit the poleward range shift in the American pika as the climate warms (Root et al. 2003; Parmesan & Yohe 2003).   In Chapter 4, I used genome-wide markers to investigate climate-driven responses in pikas distributed along two elevational gradients. One of the strengths of using genome-wide markers is the increased power to detect demographic patterns relevant to conservation (Primmer 2009). By applying a Bayesian framework to detect directional migration in a large genomic data, I assessed the fine-scale movement of pika along elevational gradients. This analysis supported previous evidence showing resource-dependent dispersal in pika where territory availability is the primary driver of movement (Peacock 1997). Importantly, since high elevation populations typically occur at relatively higher densities (Smith 1974b), elevational dispersal patterns in pikas are likely to have a consistent downslope bias. The elevational pattern of movement in pikas has important consequences to the species’ response to climate change. It is likely that low elevation pika populations will lose viability faster than their high elevation counterparts, if dispersal in pika continues to be density  104 dependent then this process will likely reinforce the downward bias of gene flow. Consequently, thermal adaptations that may evolve in lower warmer sites may be impeded from upslope movement due to the downslope direction of gene flow. Since this pattern is likely present in other alpine species, this study represents an important example of elevational gene flow in a thermally-sensitive species.   5.1.2 Evidence for cold stress American pikas typically exhibit body temperatures only a few degrees below their lethal threshold (MacArthur & Wang 1974); consequently, the thermal sensitivity of this species has received a lot of attention (Smith 1974b; Beever et al. 2011; Ray et al. 2012; Stewart et al. 2015). We resolved consistent patterns that add to the growing body of evidence for cold stress in this species (Beever et al. 2010, 2011; Ray et al. 2012; Jeffress et al. 2013; Schwalm et al. 2016). In Chapter 3, I showed a decrease in individual stress levels associated with higher ambient temperatures in early spring. Additionally, Chapter 4 highlighted several metabolic genes related to cold acclimation that appear to be under selective pressure. Antagonistic pleiotropy is when one gene impacts two or more traits that are in opposition to one another. In the case of American pikas, it seems their thermal sensitivity may be setting up an antagonistic pleiotropic relationship. For instance, some small mammals exhibit limits to seasonal metabolic plasticity (Boratyński et al. 2016). If such limits are present in pikas, then an increase in metabolic rate to endure the cold winter could leave individuals exposed to thermal stress in the hot summer since the lower limit of their summer metabolic rate would be set by the species’ metabolic plasticity. The action of natural selection on these metabolic genes could therefore leave pika disproportionately exposed to thermal stress in  105 the summer. This is a concerning result for two reasons; first, climate change predictions include an increase in seasonal extremes and second, snowpack is decreasing across the United States (Karl et al. 2009). The result of climate change could therefore leave pikas without a protective snowpack in the winter, which would necessitate further cold adaptations at the cost of increased summer heat stress. Continued monitoring of winter snowpack and thermal stress in this species remain critical steps to predict the continued viability of populations.   5.1.3 Metabolic implications The American pika has a disproportionately high metabolic rate compared to other lagomorphs and even many rodent species (Lovegrove 2003). It is likely that occupying the niche of an alpine small mammal that does not hibernate has placed unique metabolic demands on the American pika. In Chapter 3, I showed a strong relationship between body size and levels of corticosterone, a glucocorticoid stress hormone. Corticosterone is a metabolic byproduct and therefore linked to metabolic rates; among mammalian species the relationship between body size and glucocorticoids is ascribed to the species mass specific metabolic rate (Haase et al. 2016). Therefore, high levels of corticosterone in small pikas were likely due to elevated metabolic rates to compensate for heat loss due to the elevated surface area to mass ratio of smaller individuals. Interestingly, in Chapter 4, I showed evidence for several hypoxia-related genes being under natural selection from an approximate 10% reduction in available oxygen at the high-elevation sites. However, these small changes in oxygen concentrations would be exacerbated by the species’ high metabolic rate, especially in smaller individuals who exhibit elevated metabolic rates. This result adds  106 another dimension to the problem of living in a high alpine environment as a small mammal. Continued research is warranted to further delineate the role of hypoxia in elevational adaptation in the American pika. For example, the hypoxia inducible factor system, which was shown to be under directional selection in Chapter 4, matches metabolic demands of organisms to available oxygen (Pugh & Ratcliffe 2003; Schofield & Ratcliffe 2005). Delineating elevational patterns of the genes involved in the hypoxia inducible factor system could further highlight the role hypoxia has played in the evolution of American pika.    5.2 Limitations  5.2.1 Sampling limitations The initial experimental design for Chapter 4 used non-invasive hair snares previously applied to genetic assessments of the American pika (Henry & Russello 2011). It was predicted that this technique could provide ample high quality genetic material suitable for genomic analysis via RAD sequencing (Baird et al. 2008). Unfortunately, I found that the non-invasive hair sampling technique delivered two orders of magnitude less genetic material than expected and suffered from cross-species contamination (Russello et al. 2015). While this technique still provided sufficient genetic material for a high resolution genomic-based population assessment, the quality and quantity of genomic information was insufficient for a genome-wide scan for climate adaptations. Overcoming this limitation necessitated live-trapping pika to obtain tissue samples, which limited the sample sizes used in Chapters 3 and 4. Additionally, the difficulty of trapping these elusive mammals and the stress induced from trapping made direct validation of the hair corticosterone measurements used in Chapter 3 infeasible.   107 5.2.2 Genetic limitations In Chapter 2, I concluded the presence of a major highway restricted gene flow sufficiently for discrete genetic units of pika to develop. However, I also provided an alternative hypothesis that natural barriers to gene flow in the form of a seasonal stream and topography could be the driving factors of this genetic divergence. Coalescence theory provides a possible method to differential between these two competing hypotheses. By using approximate Bayesian computation, researchers have been able to determine the age of divergence between two populations (Beaumont et al. 2002; Cornuet et al. 2008). This method could be applied to these pika populations; if the age of divergence corresponds to the construction of the highway then the evidence would support anthropogenic barriers to gene flow. I attempted this calculation using the microsatellite dataset; however, the number of parameters needed to optimize the calculation inhibited the fine-scale assessment of divergence. Future work could incorporate genome-wide sampling to assess the fine scale genetic divergence associated with barriers to gene flow in the American pika and in other species suspected of being impacted by fragmentation.      While the American pika reference genome is an invaluable resource, it exists in a relatively incomplete form. Contigs are assembled into approximately 18,760 scaffolds with an average length of ~18,370 base pairs and are only sequenced at 1.93X coverage ( Additionally, large portions of repetitive elements and low complex regions are masked in this genome. During the reference assemble of sequencing data in Chapter 4, 76.1% (SD = 0.44%) of reads aligned to this genome potentially because of the genomes’ incomplete state. Moreover, the relatively short scaffolding lengths inhibited some powerful analytical approaches. Sliding window  108 analyses scan the genome for evidence of natural selection in the form of runs of decreased diversity from selective sweeps (Stephan 2016). Interestingly, I saw some evidence of this genomic pattern where sequential SNPs along several scaffolds showed outlier patterns. A more completely annotated reference genome would have further facilitated highlighting the functional role of outliers in Chapter 4.   5.2.3 Geographic limitations The environmental pressures limiting American pika occupancy are known to vary over the broad geographic range of the species. For instance, Jeffress et al. (2013) found the relative importance of specific climate stressors changed while sampling across a bioclimatic gradient encompassing eight U.S. national parks. This study showed that heat stress was most influential on pika occupancy in the driest parks while in other parks it appeared that high elevation, cold temperatures, and high precipitation acted to limit the upper elevational distribution of pika populations. Therefore, it is likely that the selective pressures being applied to pika populations and potential resulting adaptations will also vary geographically. Logistical considerations limited the geographic extent of our sampling for Chapters 3 and 4 to a relatively narrow segment of the American pika’s northern range. Given this northern distribution it is unsurprising that we primarily resolved evidence for cold stress. A broader geographic sampling of genomic patterns in the American pika is currently being conducted that aims to resolve the genomic response to climate stress in other regions of the pika’s distribution.       109 5.3 Management implications In Chapter 2, I documented restricted gene flow associated with road development and landscape modifications that potentially caused the fragmentation of pika populations around the Highland Valley Copper mine. In this chapter, I documented extensive colonization of artificial habitat that was generated by mining activities. A study contrasting pikas on artificial sites and natural talus in the Highland Valley Copper mine found no difference in body condition or survival (Blair, in preparation). This indicates that artificial talus could be used to provide habitat corridors which would promote gene flow and mitigate the barrier to movement between the northern and southern sites. Such corridors have been successful in promoting movement in other small mammals (Debinski & Holt 2000), but have rarely been applied to American pika (e.g. I-90 Snoqualmie Pass East Project; Ernest, unpublished).    In Chapter 4, I documented a natural downslope bias to movement in the American pika that may act as a potential barrier to the upslope movement of individuals as the climate warms. It is likely that this biased movement will be prevalent in other alpine species where dispersal is density dependent. One potential management response is assisted gene flow from low populations suspected to contain thermal adaptations to higher populations suspected to be suffering reduced viability due to thermal stress. Aitken & Whitlock (2013) concluded that assisted gene flow between populations has the potential to address maladaptation due to climate change when it is used to introduce genotypes that are preadapted to the new local climate or increase the frequency of these genotypes in the existing population. However, one concern that needs to be carefully considered is the risk of outbreeding depression. In the case of American pika, the risk of outbreeding depression is unlikely across the restricted geographic range of elevational gradients. We resolved  110 significant genetic differentiation between transects in Chapter 4. However, the consistent downward movement of individuals within transects likely keeps any genomic incompatibilities to a minimum. If assisted gene flow is implemented in this system, it may be important to use individuals from the closest low elevation population to reduce the potential of outbreeding depression when they are introduced to higher elevation populations.    5.4 Significance and future directions During this thesis, I evaluated a new sample collection method for the genomic assessment of elusive mammals, applied a novel technique of assessing long-term chronic stress in natural populations using hair samples, and provided one of the first genomic assessments of climate adaptation in a thermally-sensitive mammal. The work in this thesis will provide other researchers new and exciting genomic and physiological tools applicable to elusive mammals and other climate-sensitive species. Additionally, I have provided unique insights into several dispersal and metabolic patterns as well as providing further evidence for cold stress in the American pika. This knowledge adds to our understanding of biotic and climate stressors this species faces and will inform conservation interventions. However, more information is needed to more fully assess the continued viability of pika and other climate-sensitive species.  The geographic scope of the genotype-environment associations in Chapter 4 could be extended by incorporating lineage-wide genomic data. The Cascade lineage of pika extends from central British Columbia through central Oregon (Galbreath et al. 2009). A detailed genomic assessment of patterns along the latitudinal gradient of pika populations in the Cascade lineage presents the opportunity to scan for broad-scale genomic patterns  111 associated with climate variation. Such an analysis could further highlight potential adaptations in the American pika to changing climate conditions. Furthermore, evidence of lineage-wide adaptation to climate patterns could further inform population viability analyses and niche modeling exercises aimed at forecasting future range contractions both in pika and other thermally sensitive species.       112 References Aitken SN, Whitlock MC (2013) Assisted gene flow to facilitate local adaptation to climate change. Annual Review of Ecology, Evolution, and Systematics, 44, 367–388. Alexander DH, Novembre J, Lange K (2009) Fast model-based estimation of ancestry in unrelated individuals. Genome Research, 19, 1655–1664. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ (1990) Basic local alignment search tool. Journal of Molecular Biology, 215, 403–410. Baguette M, Mennechez G, Petit S, Schtickzelle N (2003) Effect of habitat fragmentation on dispersal in the butterfly Proclossiana eunomia. Comptes Rendus Biologies, 326, 200–209. Baird NA, Etter PD, Atwood TS et al. (2008) Rapid SNP discovery and genetic mapping using sequenced RAD markers. PloS one, 3, 1–7. Balding DJ, Nichols RA (1995) A method for quantifying differentiation between populations at multi-allelic loci and its implications for investigating identity and paternity. Genetica, 96, 3–12. Balkenhol N, Waits LP (2009) Molecular road ecology: exploring the potential of genetics for investigating transportation impacts on wildlife. Molecular Ecology, 18, 4151–4164. Ballantyne AP, Alden CB, Miller JB, Tans PP, White JWC (2012) Increase in observed net carbon dioxide uptake by land and oceans during the past 50 years. Nature, 488, 70–2. Bates D, Mächler M, Bolker B, Walker S (2015) Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67, 1–48. Beaumont MA, Balding DJ (2004) Identifying adaptive genetic divergence among populations from genome scans. Molecular Ecology, 13, 969–980. Beaumont MA, Nichols RA (1996) Evaluating loci for use in genetic analysis of population structure. Proceedings of the Royal Society of London Series B-Biological Sciences, 263, 1619–1626. Beaumont MA, Zhang W, Balding DJ (2002) Approximate Bayesian Computation in Population Genetics. Genetics, 162, 2025–2035. Beever EA, Brussard PF, Berger J (2003) Patterns of apparent extirpation among isolated populations of pika (Ochtona princeps) in the Great Basin. Journal of Mammalogy, 84,  113 37–54. Beever EA, Dobrowski SZ, Long J, Mynsberge AR, Piekielek NB (2013) Understanding relationships among abundance, extirpation, and climate at ecoregional scales. Ecology, 94, 1563–71. Beever EA, Hall LE, Varner J et al. (2017) Behavioral flexibility as a mechanism for coping with climate change. Frontiers in Ecology and the Environment, 15, 299–308. Beever EA, O’Leary J, Mengelt C et al. (2016) Improving conservation outcomes with a new paradigm for understanding species’ fundamental and realized adaptive capacity. Conservation Letters, 9, 131–137. Beever EA, Ray C, Mote PW, Wilkening JL (2010) Testing alternative models of climate-mediated extirpations. Ecological Applications, 20, 164–178. Beever EA, Ray C, Wilkening JL, Brussard PF, Mote PW (2011) Contemporary climate change alters the pace and drivers of extinction. Global Change Biology, 17, 2054–2070. 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. Western North American Naturalist, 68, 8–14. Belkhir K, Borsa P, Chikhi L, Raufaste N, Bonhomme F (2004) GENETIX 4.05, logiciel sous Windows TM pour la génétique des populations. Laboratoire génome populations, interactions, CNRS UMR 5000, Universite de Montpellier II, Montpellier (France). Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B, 57, 289–300. Blas J, Bortolotti GR, Tella JL, Baos R, Marchant TA (2007) Stress response during development predicts fitness in a wild, long lived vertebrate. Proceedings of the National Academy of Sciences of the United States of America, 104, 8880–8884. Bloodgood MA, Gizikoff KG, Jones CE (1998) Biosolids amendment of waste rock and tailings as reclamation treatment at Highland Valley Copper. In: 22nd British Columbia Mine Reclamation Symposium, pp. 15–25. Penticton, BC, Canada. Bonier F, Martin PR, Moore IT, Wingfield JC (2009) Do baseline glucocorticoids predict  114 fitness? Trends in Ecology and Evolution, 24, 634–642. Boratyński JS, Jefimow M, Wojciechowski MS (2016) Phenotypic flexibility of energetics in acclimated Siberian hamsters has a narrower scope in winter than in summer. Journal of Comparative Physiology B, 186, 387–402. Boutin S, Lane JE (2014) Climate change and mammals: evolutionary versus plastic responses. Evolutionary Applications, 7, 29–41. Bradshaw WE, Holzapfel CM (2001) Genetic shift in photoperiodic response correlated with global warming. Proceedings of the National Academy of Sciences of the United States of America, 98, 14509–11. Briggs D (1997) Fundamentals of the physical environment (P Smithson, K Atkinson, Eds,). Routledge, London and New York, New York. British Columbia Ministry of Transportation and Highways (2009) British Columbia Traffic Data Program. Traffic data for: Logan Lake - 28-950NS-N. Brown CJ, O’Connor MI, Poloczanska ES et al. (2016) Ecological and methodological drivers of species’ distribution and phenology responses to climate change. Global Change Biology, 22, 1548–1560. Buchmann CM, Schurr FM, Nathan R, Jeltsch F (2013) Habitat loss and fragmentation affecting mammal and bird communities—the role of interspecific competition and individual space use. Ecological Informatics, 14, 90–98. Busch DS, Hayward LS (2009) Stress in a conservation context: a discussion of glucocorticoid actions and how levels change with conservation-relevant variables. Biological Conservation, 142, 2844–2853. Cabezas S, Blas J, Marchant TA, Moreno S (2007) Physiological stress levels predict survival probabilities in wild rabbits. Hormones and Behavior, 51, 313–320. Cahill AE, Aiello-Lammens ME, Fisher-Reid MC et al. (2012) How does climate change cause extinction? Proceedings of the Royal Society B: Biological Sciences, 280, 1–9. Castillo JA, Epps CW, Davis AR, Cushman SA (2014) Landscape effects on gene flow for a climate-sensitive montane species, the American pika. Molecular Ecology, 23, 843–856. Castillo JA, Epps CW, Jeffress MR et al. (2016) Replicated landscape genetic and network analyses reveal wide variation in functional connectivity for American pikas. Ecological  115 Applications, 26, 1660–1676. Catchen JM, Amores A, Hohenlohe P, Cresko W, Postlethwait JH (2011) Stacks: building and genotyping loci de novo from short-read sequences. Genes, genomes, genetics, 1, 171–182. Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA (2013) Stacks: an analysis tool set for population genomics. Molecular Ecology, 22, 3124–3140. Cavicchi S, Guerra D, Torre V La, Raymond BH (1995) Chromosomal analysis of heat-shock tolerance in Drosophila melanogaster evolving at different temperatures in the laboratory. Evolution, 49, 676–684. Chen I, Hill JK, Ohlemüller R, Roy DB, Thomas CD (2011) Rapid range shifts of species of climate warming. Science, 333, 1024–1026. Comin A, Zufferli V, Peric T et al. (2012) Hair cortisol levels determined at different body sites in the New Zealand White rabbit. World Rabbit Science, 20, 149–154. Coop G, Witonsky D, Di Rienzo A, Pritchard JK (2010) Using environmental correlations to identify loci underlying local adaptation. Genetics, 185, 1411–23. Cornuet J-MM, Santos F, Beaumont MA et al. (2008) Inferring population history with DIY ABC: A user-friendly approach to approximate Bayesian computation. Bioinformatics, 24, 2713–2719. Daly C, Gibson WP, Taylor GH, Johnson GL, Pasteris P (2002) A knowledge-based approach to the statistical mapping of climate. Climate Research, 22, 99–113. Danecek P, Auton A, Abecasis G et al. (2011) The variant call format and VCFtools. Bioinformatics, 27, 2156–2158. Davis MB, Shaw RG, Etterson JR (2005) Evolutionary responses to changing climate. Ecology, 86, 1704–1714. Debinski DM, Holt RD (2000) A Survey and Overview of Habitat Fragmentation Experiments. Conservation Biology, 14, 342–355. Devlin B, Roeder K (1999) Genomic control for association studies. Biometrics, 55, 997–1004. Duforet-Frebourg N, Blum MGB (2014) Nonstationary patterns of isolation-by-distance: inferring measures of local genetic differentiation with Baysian kriging. Evolution, 68,  116 1110–1123. Duke K (1951) The external genitalia of the pika. Journal of Mammalogy, 32, 169–173. Durka W (1999) Genetic diversity in peripheral and subcentral populations of Corrigiola litoralis L. (Illecebraceae). Heredity, 83, 476–484. van Dyck H, Baguette M (2005) Dispersal behaviour in fragmented landscapes: Routine or special movements? Basic and Applied Ecology, 6, 535–545. Earl DA, VonHoldt BM (2011) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genetics Resources, 4, 359–361. Eckert CG, Samis KE, Lougheed SC (2008) Genetic variation across species’ geographical ranges: the central–marginal hypothesis and beyond. Molecular Ecology, 17, 1170–1188. Ellegren H (2014) Genome sequencing and population genomics in non-model organisms. Trends in Ecology & Evolution, 29, 51–63. Estes-Zumpf WA, Rachlow JL, Waits LP, Warheit KI (2010) Dispersal, gene flow, and population genetic structure in the pygmy rabbit (Brachylagus idahoensis). Journal of Mammalogy, 91, 208–219. Etter PD, Bassham S, Hohenlohe PA, Johnson EA, Cresko WA (2011) Molecular Methods for Evolutionary Genetics. In: Molecular Methods for Evolutionary Genetics Methods in Molecular Biology. (eds Orgogozo V, Rockman M V.), pp. 157–178. Humana Press, Totowa, NJ. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 14, 2611–20. Excoffier L, Lischer HEL (2010) Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources, 10, 564–7. Fischer J, Lindenmayer DB (2007) Landscape modification and habitat fragmentation: a synthesis. Global Ecology and Biogeography, 16, 265–280. Foll M, Gaggiotti O (2008) A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics, 180, 977–93. François O, Martins H, Caye K, Schoville SD (2016) Controlling false discoveries in genome  117 scans for selection. Molecular Ecology, 25, 454–469. Frankham R (1995) Effective population size/adult population size ratios in wildlife: a review. Genetical Research, 66, 95–107. Frankham R (1996) Relationship of genetic variation to population size in wildlife. Conservation Biology, 10, 1500–1508. Franks SJ, Kane NC, O’Hara NB, Tittes S, Rest JS (2016) Rapid genome-wide evolution in Brassica rapa populations following drought revealed by sequencing of ancestral and descendant gene pools. Molecular Ecology, 25, 3622–31. Franks SJ, Sim S, Weis AE (2007) Rapid evolution of flowering time by an annual plant in response to a climate fluctuation. Proceedings of the National Academy of Sciences, 104, 1278–1282. Freberg MR, Gizikoff KG (1999) Development and utilization of an end land use plan for highland valley copper. In: Proceendings of the 23rd Annual British Columbia Mine Reclamation Symposium, pp. 46–56. Kamploops, BC, Canada. Frichot E, François O (2015) LEA: An R package for landscape and ecological association studies. Methods in Ecology and Evolution, 6, 925–929. Frichot E, Schoville SD, Bouchard G, François O (2013) Testing for associations between loci and environmental gradients using latent factor mixed models. Molecular Biology and Evolution, 30, 1687–1699. Funk WC, McKay JK, Hohenlohe PA, Allendorf FW (2012) Harnessing genomics for delineating conservation units. Trends in Ecology & Evolution, 27, 489–496. Gaines MS, Diffendorfer JE, Tamarin RH, Whitttam TS (1997) The effects of habitat fragmentation on the genetic structure of small mammal populations. Journal of Heredity, 88, 294–304. Galbreath KE, Hafner DJ, Zamudio KR (2009) When cold is better: climate-driven elevation shifts yield complex patterns of diversification and demographi in an alpine specialist (American pika, Ochotona princeps). Evolution, 63, 2848–2863. Gardiner SM (2012) Ethics and Global Climate Change. Nature Education Knowledge, 3, 1–7. Gautier M (2015) Genome-wide scan for adaptive divergence and association with  118 population-specific covariates. Genetics, 201, 1555–1579. Gienapp P, Lof M (2013) Predicting demographically sustainable rates of adaptation: can great tit breeding time keep pace with climate change? Philosophical transactions of the Royal Society of London. Series B, Biological sciences, 368, 1–10. Gienapp P, Teplitsky C, Alho JS, Mills JA, Merilä J (2008) Climate change and evolution: disentangling environmental and genetic responses. Molecular ecology, 17, 167–78. Gilk WR, Richardson S, Spiegelhalter DJ (1996) Markov Chain Monte Carlo in practice. Chapman and Hall, London/New York. Gittleman JL, Thompson D (1988) Energy allocation in mammalian reproduction. Integrative and Comparative Biology, 28, 863–875. Gow R, Thomson S, Rieder M, Van Uum S, Koren G (2010) An assessment of cortisol analysis in hair and its clinical applications. Forensic Science International, 196, 32–37. Grayson DK (2005) A brief history of Great Basin pikas. Journal of Biogeography, 32, 2103–2111. Günther T, Coop G (2013) Robust identification of local adaptation from allele frequencies. Genetics, 195, 205–20. Guo B, Lu D, Liao WB, Merilä J (2016) Genomewide scan for adaptive differentiation along altitudinal gradient in the Andrew’s toad Bufo andrewsi. Molecular Ecology, 25, 3884–3900. Guo SW, Thompson EA (1992) Performing the exact test of Hardy-Weinberg proportion for multiple alleles. Biometrics, 48, 361–372. Haase CG, Long AK, Gillooly JF (2016) Energetics of stress : linking plasma cortisol levels to metabolic rate in mammals. Biology letters, 12, 20150867. Hafner DJ (1993) North American pika (Ochotona princeps) as a late Quaternary biogeographic indicator species. Quaternary Research, 39, 373–380. Hafner DJ, Sullivan RM (1995) Historical and ecological biogeography of neararctic pikas (Lagomorpha: Ochotonidae). Journal of Mammalogy, 76, 302–321. Helyar SJ, Hemmer-Hansen J, Bekkevold D et al. (2011) Application of SNPs for population genetics of nonmodel organisms: new opportunities and challenges. Molecular ecology resources, 11 Suppl 1, 123–36. Hendry AP, Farrugia TJ, Kinnison MT (2008) Human influences on rates of phenotypic  119 change in wild animal populations. Molecular Ecology, 17, 20–29. Henry P, Russello MA (2011) Obtaining high-quality DNA from elusive small mammals using low-tech hair snares. European Journal of Wildlife Research, 57, 429–435. Henry P, Russello MA (2013) Adaptive divergence along environmental gradients in a climate-change-sensitive mammal. Ecology and Evolution, 3, 3906–17. Henry P, Sim Z, Russello MA (2012) Genetic evidence for restricted dispersal along continuous altitudinal gradients in a climate change-sensitive mammal: the American pika. PloS one, 7, e39077. Hoffman RS, Smith AT (2005) Order Lagomorpha. In: Mammal Species of the World: A Taxonomic and Geographic Reference (eds Wilson DE, Reeder DM), pp. 185–193. The Johns Hopkins University Press, Baltimore. Hoffmann A, Hallas RJ, Dean JA, Schiffer M (2003) Low potential for climatic stress adaptation in a rainforest Drosophila species. Science, 301, 100–102. Hoffmann AA, Sgrò CM (2011) Climate change and evolutionary adaptation. Nature, 470, 479–485. Hoffmann AA, Willi Y (2008) Detecting genetic responses to environmental change. Nature Reviews Genetics, 9, 421–432. Hohenlohe PA, Bassham S, Etter PD et al. (2010) Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags. PLoS Genetics, 6, 1–23. Howie R (2007) Habitat locations for the American pika (Ochotona princeps) in the Highmont area. Aspen Park Consulting. Kamloops, BC, Canada, 1–13. Hutchinson MF (1989) A new objective method for spatial interpolation of meteorological variables from irregular networks applied to the estimation of monthly mean solar radiation, temperature, precipitation and windrun. Division of Water Resources. CSIRO Australia. Technical Memorandum 89/5, 104 p. 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. Ecological Applications, 23, 864–878. Jeffreys H (1961) The theory of probability. Oxford University Press, London/New York/Oxford. Jensen JL, Bohonak AJ, Kelley ST (2005) Isolation by distance, web service. BMC genetics,  120 6, 13. Jones MR, Forester BR, Teufel AI et al. (2013) Integrating landscape genomics and spatially explicit approaches to detect loci under selection in clinal populations. Evolution, 67, 3455–3468. Jones PD, Osborn TJ, Briffa KR (2001) The evolution of climate over the last millennium. Science, 292, 662–667. Joost S, Bonin A, Bruford MW et al. (2007) A spatial analysis method (SAM) to detect candidate loci for selection: towards a landscape genomics approach to adaptation. Molecular Ecology, 16, 3955–3969. Kalinowski ST (2005) Hp-Rare 1.0: A computer program for performing rarefaction on measures of allelic richness. Molecular Ecology Notes, 5, 187–189. Karl TR, Melillo JM, Peterson TC (2009) Global climate change impacts in the United States. Cambridge University Press, Cambridge. Keller LF, Waller DM (2002) Inbreeding effects in wild populations. Trends in Ecology and Evolution, 17, 230–241. Keyghobadi N (2007) The genetic implications of habitat fragmentation for animals. Canadian Journal of Zoology, 85, 1049–1064. Koren L, Mokady O, Geffen E (2008) Social status and cortisol levels in singing rock hyraxes. Hormones and Behavior, 54, 212–216. Koren L, Mokady O, Karaskov T et al. (2002) A novel method using hair for detirmining hormal levels in wildlife. Animal Behaviour, 63, 403–406. Krear HR (1965) An ecological and ethological study of the pika (Ochotona princeps saxatilis Bangs) in the Front Range of Colorado. PhD Dissertation. University of Colorado, 329 p. Kreuzer MP, Huntly NJ (2003) Habitat-specific demography: evidence for source-sink population structure in a mammal, the pika. Oecologia, 134, 343–349. Lamb CT, Robson KM, Russello MA (2013) Development and application of a molecular sexing protocol in the climate change-sensitive American pika. Conservation Genetics Resources, 6, 17–19. Lanier HC, Olson LE (2009) Inferring divergence times within pikas (Ochotona spp.) using mtDNA and relaxed molecular dating techniques. Molecular Phylogenetics and  121 Evolution, 53, 1–12. Leberg PL (2002) Estimating allelic richness: Effects of sample size. Molecular Ecology, 11, 2445–2449. Lefcheck JS (2016) piecewiseSEM: Piecewise structural equation modelling in R for ecology, evolution, and systematics. Methods in Ecology and Evolution, 7, 573–579. Lemay MA, Henry P, Lamb CT, Robson KM, Russello MA (2013) Novel genomic resources for a climate change sensitive mammal: characterization of the American pika transcriptome. BMC genomics, 14, 1–11. Lemay MA, Russello MA (2015) Genetic evidence for ecological divergence in kokanee salmon. Molecular Ecology, 24, 798–811. Lesica P, Allendorf FW (1995) When are peripheral valuable populations for conservation? Conservation Biology, 9, 753–760. Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 25, 1754–1760. Li H, Handsaker B, Wysoker A et al. (2009a) The sequence alignment/map format and SAMtools. Bioinformatics, 25, 2078–2079. Li H, Ren Y, Guo S et al. (2009b) The protein level of hypoxia-inducible factor-1α is increased in the plateau pika (Ochotona curzoniae) inhabiting high altitudes. Journal of Experimental Zoology Part A: Ecological Genetics and Physiology, 311A, 134–141. Lotterhos KE, Whitlock MC (2015) The relative power of genome scans to detect local adaptation depends on sampling design and statistical method. Molecular Ecology, 24, 1031–1046. Lovegrove BG (2003) The influence of climate on the basal metabolic rate of small mammals: a slow-fast metabolic continuum. Journal of Comparative Physiology B, 173, 87–112. Lovejoy TE, Hannah L (2005) Climate change and biodiversity. Yale University Press, New Haven CT. Luikart G, England PR, Tallmon D, Jordan S, Taberlet P (2003) The power and promise of population genomics: from genotyping to genome typing. Genetics, 4, 981–94. Luikart G, Ryman N, Tallmon DA, Schwartz MK, Allendorf FW (2010) Estimation of census and effective population sizes: the increasing usefulness of DNA-based  122 approaches. Conservation Genetics, 11, 355–373. MacArthur RA, Wang LCH (1973) Physiology of Thermoregulation in the pika. Canadian Journal of Zoology, 51, 11–16. MacArthur A, Wang CH (1974) Behavioral thermoregulation in the pika Ochotona princeps: a field study using radiotelemetry. Canadian Journal of Zoology, 52, 353–358. Macbeth BJ, Cattet MRL, Stenhouse GB, Gibeau ML, Janz DM (2010) Hair cortisol concentration as a noninvasive measure of long-term stress in free-ranging grizzly bears (Ursus arctos): considerations with implications for other wildlife. Canadian Journal of Zoology, 88, 935–949. Martínez-Cruz B, Godoy JA, Negro JJ (2007) Population fragmentation leads to spatial and temporal genetic structure in the endangered Spanish imperial eagle. Molecular ecology, 16, 477–86. Mastromonaco GF, Gunn K, Edwards DB (2014) Validation and use of hair cortisol as a measure of chronic stress in eastern chipmunks (Tamias striatus). Conservation Physiology, 2, 1–12. Mayr E (1963) Animal species and evolution. Harvard University Press, Cambridge, MA. McEwen BS, Wingfield JC (2003) The concept of allostasis in biology and biomedicine. Hormones and Behavior, 43, 2–15. Meirmans PG, Van Tienderen PH (2004) GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms. Molecular Ecology Notes, 4, 792–794. Meredith SJ (2002) The impact of habitat spatial structure on pika (Ochotona princeps) dispersal dynamics. MSc Thesis. University of Nevada, Reno, NV, 62p. Merilä J (2012) Evolution in response to climate change: in pursuit of the missing evidence. BioEssays, 34, 811–818. Merilä J, Hendry AP (2014) Climate change, adaptation, and phenotypic plasticity: the problem and the evidence. Evolutionary Applications, 7, 1–14. Meyer J, Novak M, Hamel A, Rosenberg K (2014) Extraction and analysis of cortisol from human and monkey hair. Journal of Visualized Experiments, 83, 1–6. Millspaugh JJ, Washburn BE (2004) Use of fecal glucocorticoid metabolite measures in conservation biology research: considerations for application and interpretation.  123 General and Comparative Endocrinology, 138, 189–199. De Mita S, Thuillet AC, Gay L et al. (2013) Detecting selection along environmental gradients: analysis of eight methods and their effectiveness for outbreeding and selfing populations. Molecular Ecology, 22, 1383–1399. Moilanen A, Smith AT, Hanski I (1998) Long-term dynamics in a metapopulation of the American pika. American Naturalist, 152, 530–542. Monclús R, Rödel HG, Palme R, Von Holst D, De Miguel J (2006) Non-invasive measurement of the physiological stress response of wild rabbits to the odour of a predator. Chemoecology, 16, 25–29. Moyer-Horner L, Mathewson PD, Jones GM, Kearney MR, Porter WP (2015) Modeling behavioral thermoregulation in a climate change sentinel. Ecology and Evolution, 1–13. Muir AP, Biek R, Thomas R, Mable BK (2014) Local adaptation with high gene flow: temperature parameters drive adaptation to altitude in the common frog (Rana temporaria). Molecular Ecology, 23, 561–574. Nakagawa S, Schielzeth H (2013) A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution, 4, 133–142. Narum SR, Buerkle CA, Davey JW, Miller MR, Hohenlohe PA (2013) Genotyping-by-sequencing in ecological and conservation genomics. Molecular Ecology, 22, 2841–2847. Nei M (1978) Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics, 89, 583–590. Nei M, Maruyama T, Chakraborty R (1975) The bottleneck effect and genetic variability in population. Evolution, 29, 1–10. Nieuwenhuis R, Te Grotenhuis M, Pelzer B (2012) influence.ME: tools for detecting influential data in mixed effects models. R Journal, 4, 38–47. Noël S, Ouellet M, Galois P, Lapointe F-J (2006) Impact of urban fragmentation on the genetic structure of the eastern red-backed salamander. Conservation Genetics, 8, 599–606. Nozawa Y (2011) Adaptive regulation of membrane lipids and fluidity during thermal  124 acclimation in Tetrahymena. Proceedings of the Japan Academy, Series B, 87, 450–462. O’Grady JJ, Brook BW, Reed DH et al. (2006) Realistic levels of inbreeding depression strongly affect extinction risk in wild populations. Biological Conservation, 133, 42–51. Oleksyk TK, Smith MW, O’Brien SJ (2010) Genome-wide scans for footprints of natural selection. Philosophical transactions of the Royal Society of London. Series B, Biological sciences, 365, 185–205. van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes, 4, 535–538. Opdam P, Wascher D (2004) Climate change meets habitat fragmentation: linking landscape and biogeographical scale levels in research and conservation. Biological Conservation, 117, 285–297. Osborn SD, Applebee DB (2011) Evalutation of the petition from the Center for Biological Diversity to list the American pika (Ochotona princeps) as threatened. California Department of Fish and Game Memorandum, 32p. Ozgul A, Childs DZ, Oli MK et al. (2010) Coupled dynamics of body mass and population growth in response to environmental change. Nature, 466, 482–5. Pacifici M, Visconti P, Butchart SHM et al. (2017) Species’ traits influenced their response to recent climate change. Nature Climate Change, 7, 205–208. Parmesan C (2006) Ecological and evolutionary responses to recent climate change. Annual Review of Ecology, Evolution, and Systematics, 37, 637–669. Parmesan C, Ryrholm N, Stefanescu C et al. (1999) Poleward shifts in geographical ranges of butterfly species associated with regional warming. Nature, 399, 579–583. Parmesan C, Yohe G (2003) A globally coherent fingerprint of climate change impacts across natural systems. Nature, 421, 37–42. Patterson N, Price AL, Reich D (2006) Population structure and eigenanalysis. PLoS Genetics, 2, 2074–2093. Peacock MM (1997) Determining natal dispersal patterns in a population of North American pikas (Ochotona princeps) using direct mark-resight and indirect genetic methods.  125 Behavioral Ecology, 8, 340–350. Peacock AJ (1998) ABC of oxygen: oxygen at high altitude. BMJ, 317, 1063–1066. Peacock MM, Kirchoff VS, Merideth SJ (2002) Identification and characterization of nine polymorphic microsatellite loci in the North American pika, Ochotona princeps. Molecular Ecology Notes, 2, 360–362. Peacock MM, Smith AT (1997a) The effect of habitat fragmentation on dispersal patterns, mating behavior, and genetic variation in a pika (Ochotona princeps) metapopulation. Oecologia, 112, 524–533. Peacock MM, Smith AT (1997b) Nonrandom mating in pikas Ochotona princeps: evidence for inbreeding between individuals of intermediate relatedness. Molecular Ecology, 6, 801–811. Peakall R, Smouse PE (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes, 6, 288–295. Pearson RG, Dawson TP (2003) Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful? Global Ecology and Biogeography, 12, 361–371. Peric T, Comin A, Corazzin M et al. (2017) Relocation and hair cortisol concentrations in New Zealand white rabbits. Journal of Applied Animal Welfare Science, 20, 1–8. Piry S, Alapetite A, Cornuet JM et al. (2004) GENECLASS2: A software for genetic assignment and first-generation migrant detection. Journal of Heredity, 95, 536–539. Pop M, Salzberg SL (2008) The impact of next-generation sequencing technology on genetics. Trends in Genetics, 24, 133–141. Primmer CR (2009) From conservation genetics to conservation genomics. Annals of the New York Academy of Sciences, 1162, 357–368. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics, 945–959. Pugh CW, Ratcliffe PJ (2003) Regulation of angiogenesis by hypoxia: role of the HIF system. Nature Medicine, 9, 677–684. Van der Putten WH, Macel M, Visser ME (2010) Predicting species distribution and abundance responses to climate change: why it is essential to include biotic interactions across trophic levels. Philosophical Transactions of the Royal Society B: Biological  126 Sciences, 365, 2025–2034. Queller DC, Goodnight KF (1989) Estimating relatedness using genetic markers. Evolution, 43, 258–275. R Core Team (2015) R: A language and environment for statstical computing. R Foundation for Statstical Computing, Vienna, Austria. Ramirez J-M, Folkow LP, Blix AS (2007) Hypoxia tolerance in mammals and birds: from the wilderness to the clinic. Annual Review of Physiology, 69, 113–143. Rannala B, Mountain JL (1997) Detecting immigration by using multilocus genotypes. Proceedings of the National Academy of Sciences, 94, 9197–9201. Ray C, Beever EA, Loarie SA (2012) Retreat of the American pika: up the mountain or into the void? In: Conserving wildlife populations in a changing climate. (eds Brodie JF, Post E, Doak D), pp. 245–270. University of Chicago Press, Chicago. Rayman M, Rousset F (1995) GENEPOP (version 1.2): population genetics software for exact test and ecumenicism. Journal of Heredity, 86, 248–249. Reeder DM, Kramer KM (2005) Stress in free-ranging mammals: integrating physiology, ecology, and natural history. Journal of Mammalogy, 86, 225–235. Rellstab C, Gugerli F, Eckert AJ, Hancock AM, Holderegger R (2015) A practical guide to environmental association analysis in landscape genomics. Molecular Ecology, 24, 4348–4370. Reusch TBH, Wood TE (2007) Molecular ecology of global change. Molecular Ecology, 16, 3973–3992. Rice WR (1989) Analyzing tables of statistical tests. Evolution, 43, 223–225. 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. Journal of Mammalogy, 97, 464–472. Rodhouse TJ, Hovland M, Jeffress MR (2017) Variation in subsurface thermal characteristics of microrefuges used by range core and peripheral populations of the American pika (Ochotona princeps). Ecology and Evolution, 7, 1514–1526. Root TL, Price JT, Hall KR, Schneider SH (2003) Fingerprints of global warming on wild animals and plants. Nature, 421, 57–60. Rousset F (2008) genepop’007: a complete re-implementation of the genepop software for  127 Windows and Linux. Molecular Ecology Resources, 8, 103–106. Ruiz-Jaen MC, Aide TM (2005) Restoration success: how is it being measured? Restoration Ecology, 13, 569–577. Russell E, Koren G, Rieder M, Van Uum S (2012) Hair cortisol as a biological marker of chronic stress: current status, future directions and unanswered questions. Psychoneuroendocrinology, 37, 589–601. Russello MA, Waterhouse MD, Etter PD, Johnson EA (2015) From promise to practice: pairing non-invasive sampling with genomics in conservation. PeerJ, 3, e1106. Sapolsky RM, Romero LM, Munck AU (2000) How do glucocorticoids influence stress responses? Integrating permissive, supressive, stimulatory and preparative actions. Endocrine Reviews, 21, 55–89. Saunders DA, Hobbs RJ, Margules CR (1991) Biological consequences of ecosystem fragmentation: a review. Conservation Biology, 5, 18–32. Sauvajot RM, Buechner M, Kamradt DA, Schonewald CM (1998) Patterns of human disturbance and response by small mammals and birds in chaparral near urban development. Urban Ecosystems, 2, 279–297. Savolainen O, Lascoux M, Merilä J (2013) Ecological genomics of local adaptation. Nature Reviews Genetics, 14, 807–820. Schilthuizen M, Kellermann V (2014) Contemporary climate change and terrestrial invertebrates: evolutionary versus plastic changes. Evolutionary Applications, 7, 56–67. Schofield CJ, Ratcliffe PJ (2005) Signalling hypoxia by HIF hydroxylases. Biochemical and Biophysical Research Communications, 338, 617–626. Schwalm D, Epps CW, Rodhouse TJ et al. (2016) Habitat availability and gene flow influence diverging local population trajectories under scenarios of climate change: a place-based approach. Global Change Biology, 22, 1572–1584. Sgrò CM, Lowe AJ, Hoffmann AA (2011) Building evolutionary resilience for conserving biodiversity under climate change. Evolutionary Applications, 4, 326–337. Sheriff MJ, Dantzer B, Delehanty B, Palme R, Boonstra R (2011) Measuring stress in wildlife: techniques for quantifying glucocorticoids. Oecologia, 166, 869–887. Sheriff MJ, Krebs CJ, Boonstra R (2010) Assessing stress in animal populations: do fecal and plasma glucocorticoids tell the same story? General and Comparative Endocrinology,  128 166, 614–619. Sheriff MJ, Wheeler H, Donker SA et al. (2012) Mountain-top and valley-bottom experiences: the stress axis as an integrator of environmental variability in arctic ground squirrel populations. Journal of Zoology, 287, 65–75. Simpson WG (2009) American pikas inhabit low-elevation sites outside the species’ previously described bioclimatic envelope. Western North American Naturalist, 69, 243–250. Smith AT (1974a) The distribution and dispersal of pikas: consquences of insular population structure. Ecology, 55, 1112–1119. Smith AT (1974b) The distribution and dispersal of pikas: influences of behavior and climate. Ecology, 55, 1368–1376. Smith AT (1980) Temporal changes in insular populations of the pika (Ochotona princeps). Ecology, 61, 8–13. Smith JM, Haigh J (1974) The hitch-hiking effect of a favourable gene. Genetical Research, 23, 23–35. Smith AT, Ivins BL (1983) Colonization in a pika population: dispersal vs philopatry. Behavioral Ecology and Sociobiology, 13, 37–47. Smith O, Wang J (2014) When can noninvasive samples provide sufficient information in conservation genetics studies? Molecular Ecology Resources, 14, 1011–1023. Smith AT, Weston ML (1990) Ochotona princeps. Mammalian Species, 352, 1–8. Solomon S, Plattner G-K, Knutti R, Friedlingstein P (2009) Irreversible climate change due to carbon dioxide emissions. Proceedings of the National Academy of Sciences of the United States of America, 106, 1704–1709. Somero GN (2010) The physiology of climate change: how potentials for acclimatization and genetic adaptation will determine “winners” and “losers”. The Journal of Experimental Biology, 213, 912–20. Stephan W (2016) Signatures of positive selection: from selective sweeps at individual loci to subtle allele frequency changes in polygenic adaptation. Molecular Ecology, 25, 79–88. Stewart JAE, Perrine JD, Nichols LB et al. (2015) Revisiting the past to foretell the future: summer temperature and habitat area predict pika extirpations in California. Journal of  129 Biogeography, 42, 880–890. Stocker TF, Qin D, Plattner GK et al. (2013) Summary for Policymakers. In: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge. Sundqvist MK, Sanders NJ, Wardle DA (2013) Community and ecosystem responses to elevational gradients: processes, mechanisms, and insights for global change. Annual Review of Ecology, Evolution, and Systematics, 44, 261–280. Taylor AC, Walker FM, Goldingay RL et al. (2011) Degree of landscape fragmentation influences genetic isolation among populations of a gliding mammal. PloS one, 6, e26651. Teck Resources Limited (2012) Generations: Teck Sustainability Report. Vancouver, BC, Canada. Telonis-Scott M, Guthridge KM, Hoffmann AA (2006) A new set of laboratory-selected Drosophila melanogaster lines for the analysis of desiccation resistance: response to selection, physiology and correlated responses. The Journal of Experimental Biology, 209, 1837–1847. Tempel DJ, Gutierrez RJ (2004) Factors related to fecal corticosterone levels in California Spotted Owls: implications for assessing chronic stress. Conservation Biology, 18, 538–547. Tigano A, Friesen VL (2016) Genomics of local adaptation with gene flow. Molecular Ecology, 25, 2144–2164. Touma C, Palme R (2005) Measuring fecal glucocorticoid metabolites in mammals and birds: the importance of validation. Annals of the New York Academy of Sciences, 1046, 54–74. United States Fish and Wildlife Service (USFWS) (2010) Endangered and threatened wildlife and plants; 12-month finding on a petition to list the American pika as threatened or endangered. CFR 17, 1-34. Urban MC, De Meester L, Vellend M, Stoks R, Vanoverbeke J (2012) A crucial step toward realism: responses to climate change from an evolving metacommunity perspective. Evolutionary Applications, 5, 154–167. VanDerWal J, Murphy HT, Kutt AS et al. (2012) Focus on poleward shifts in species’  130 distribution underestimates the fingerprint of climate change. Nature Climate Change, 3, 239–243. Varner J, Dearing MD (2014) The importance of biologically relevant microclimates in habitat suitability assessments. PloS one, 9, e104648. De Villemereuil P, Frichot É, Bazin É, François O, Gaggiotti OE (2014) Genome scan methods against more complex models: when and how much should we trust them? Molecular Ecology, 23, 2006–2019. Visser ME (2008) Keeping up with a warming world; assessing the rate of adaptation to climate change. Proceedings of the Royal Society B: Biological Sciences, 275, 649–659. Vitalis R, Dawson K, Boursot P (2001) Interpretation of variation across marker loci as evidence of selection. Genetics, 158, 1811–1823. Vitti JJ, Grossman SR, Sabeti PC (2013) Detecting natural selection in genomic data. Annual Review of Genetics, 47, 97–120. Wade GN, Schneider JE (1992) Metabolic fuels and reproduction in female mammals. Neuroscience and Biobehavioral Reviews, 16, 235–272. Walther G, Post E, Convey P et al. (2002) Ecological responses to recent climate change. Nature, 416, 389–395. Wang T, Hamann A, Spittlehouse DL, Murdock TQ (2012) ClimateWNA—High-resolution spatial climate data for western North America. Journal of Applied Meteorology and Climatology, 51, 16–29. Wang J-M, Zhang Y-M, Wang D-H (2006) Seasonal thermogenesis and body mass regulation in plateau pikas (Ochotona curzoniae). Oecologia, 149, 373–82. Waterhouse MD, Sjodin B, Ray C et al. (2017) Individual-based analysis of hair corticosterone reveals factors influencing chronic stress in the American pika. Ecology and Evolution, 7, 4099–4108. Weir BS, Cockerham CC (1984) Estimating F-statistics for the analysis of population structure. Evolution, 38, 1358–1370. Whitlock MC (2005) Combining probability from independent tests: the weighted Z-method is superior to Fisher’s approach. Journal of Evolutionary Biology, 18, 1368–1373. Wilcox BA, Murphy DD (1985) Conservation strategy: the effects of fragmentation on  131 extinction. American Naturalist, 125, 879–887. Wilkening JL, Ray C (2016) Characterizing predictors of survival in the American pika (Ochotona princeps). Journal of Mammalogy, 97, 1366–1375. Wilkening JL, Ray C, Beever EA, Brussard PF (2011) Modeling contemporary range retraction in Great Basin pikas (Ochotona princeps) using data on microclimate and microhabitat. Quaternary International, 235, 77–88. Wilkening JL, Ray C, Sweazea KL (2013) Stress hormone concentration in Rocky Mountain populations of the American pika (Ochotona princeps). Conservation Physiology, 1, 1–13. Wilson GA, Rannala B (2003) Bayesian inference of recent migration rates using multilocus genotypes. Genetics, 163, 1177–1191. Wingfield JC, Maney DL, Breuner CW et al. (1998) Ecological bases of hormone-behavior interactions: the “Emergency Life History Stage.” American Zoologist, 38, 191–206. Wingfield JC, Romero LM (2011) Adrenocortical responses to stress and their modulation in free-living vertebrates. In: Comprehensive Physiology, pp. 211–236. John Wiley & Sons, Inc., Hoboken, NJ, USA. Yandell M, Ence D (2012) A beginner’s guide to eukaryotic genome annotation. Nature Reviews Genetics, 13, 329–342. Yandow LH, Chalfoun AD, Doak DF (2015) Climate tolerances and habitat requirements jointly shape the elevational distribution of the American pika (Ochotona princeps), with implications for climate change effects. PLoS ONE, 10, 1–21. Yang HZ, Lan J, Yan JM, Xue JW, Dail WH (1998) A preliminary study of steroid reproductive hormones in human hair. Journal of Steroid Biochemistry and Molecular Biology, 67, 447–450. Yang J, Wang ZL, Zhao XQ et al. (2008) Natural selection and adaptive evolution of leptin in the Ochotona family driven by the cold environmental stress. PloS one, 3, e1472. Zaphiropoulos PG (1997) Exon skipping and circular RNA formation in transcripts of the human cytochrome P-450 2C18 gene in epidermis and of the rat androgen binding protein gene in testis. Molecular and Cellular Biology, 17, 2985–2993.   


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items