Stress-associated metabolites vary with both season and habitat across populations of a climate sentinel species

ABSTRACT

Relating physiological stress to habitat quality could refine conservation efforts. Habitat quality, which is often inferred from patch occupancy or demographic rates, might be measured in a more timely and nuanced way using metrics of physiological stress. To understand whether stress-associated hormones vary with metrics of habitat quality, we measured fecal glucocorticoid metabolite (FGM) levels in the American pika (Ochotona princeps), a small mammal with well-defined habitat (talus), which can vary in quality depending on the presence of subsurface ice features. In spring and fall 2018, we collected feces noninvasively from pika territories in taluses “with” or “without” subsurface ice to capture seasonal variation in FGM between habitat types. We used linear mixed effects models to explore the interactions among season, habitat metrics (including subsurface ice status), and subsurface temperature as predictors of FGM. We found support for interacting effects on FGM levels, which covaried with season, elevation, putative ice presence, graminoid to forb ratio, graminoid cover, and measures of acute subsurface heat exposure. However, only one subsurface temperature metric differed according to putative presence of subsurface ice. Our results contribute to the growing evidence that FGMs might be developed as a tool to assess habitat quality.

KEYWORDS:

• Ochotona princeps saxatilis
• habitat quality
• stress-associated hormone
• glucocorticoid metabolite
• rock glacier
• American pika

Introduction

Climate change is predicted to influence range dynamics for many species over the coming decades, with general trends showing upslope and poleward shifts in species distributions (Chen et al. 2011). However, developing methods for evaluating these predictions in the short term, rather than waiting for range adjustments, would speed the refinement of conservation efforts. Wildlife management has long relied on monitoring techniques that measure population size and distribution to understand population health. However, these techniques can only detect population decline after it has already happened. For example, in a review of 950 papers on relationships of vertebrate and invertebrate animals to habitat fragmentation (Ellis, McWhorter, and Maron 2012), the majority of studies focused on distribution-related responses such as presence/absence or abundance. In contrast, only thirty-six studies (<4 percent) measured a physiological response, like stress, to habitat factors (Ellis, McWhorter, and Maron 2012).

Wildlife populations live in a patchwork of habitats across landscapes and constantly experience changes in localized environments. How animals cope with these changing conditions can affect their fitness (Sheriff et al. 2011). Animals in low-quality habitats may not have access to the energy required for allostasis (physiological stability in the presence of actual or perceived stressors), resulting in allostatic overload (a cumulative effect of chronic stress) that reduces individual fitness (McEwen and Wingfield 2003). Thus, stress in individuals can be influenced by habitat quality. By evaluating patterns of the individual stress response among populations, we can gain insights into drivers of population-level health (Homyack 2010; Ellis, McWhorter, and Maron 2012; Madliger and Love 2014; Lea et al. 2018).

When an unpredictable event (e.g., predator–prey interaction, extreme weather, food shortage) challenges an animal’s physiological and behavioral stability (homeostasis), a cascade of physiological processes is activated to regain equilibrium (McEwen and Wingfield 2003; Reeder and Kramer 2005). One part of this process is the release of glucocorticoid hormones, which aid in modulating the animal’s response to the current stressor (Sapolsky, Romero, and Munck 2000; MacDougall-Shackleton et al. 2019). A rapid glucocorticoid response is often beneficial, because it stimulates a rapid physical and/or behavioral response (like running from a predator, seeking shelter from inclement weather, or moving to find better resources) and thus speeds the subsequent process of replenishing energy and restoring homeostasis (McEwen and Wingfield 2003). However, if the stressor remains persistent, prolonged elevation of glucocorticoid levels can lead to decreased body mass, reduced reproduction, and increased mortality (Romero, Dickens, and Cyr 2009; Dantzer et al. 2014). Metrics of physiological stress therefore have the potential to help conservationists understand habitat-mediated effects on populations (Homyack 2010). For example, in red-legged kittiwakes (Rissa brevirostris), juveniles who were fed poor-quality diets grew more slowly and exhibited higher levels of plasma corticosterone hormones (Kitaysky et al. 2006). Likewise, fecal corticosterone metabolite levels were higher and body mass was lower in wild wood mice (Apodemus sylvaticus) living in poor quality pine forest habitats that provided fewer resources, compared to those in higher quality holm oak wood habitats (Navarro-Castilla and Barja 2019).

Understanding the habitat features that result in higher average stress at the population level could help managers prioritize conservation efforts. For example, studies on the influence of traffic noise on stress hormone levels in spotted owls (Strix occidentalis caurina and Strix occidentalis occidentalis) helped inform conservation strategies for the species and highlighted both potential costs and benefits of roads on population dynamics (USDA Forest Service 2019). Here, we investigate whether habitat quality can influence stress in American pikas (Ochotona princeps). Specifically, we ask whether individual fecal glucocorticoid metabolites (FGMs) vary with metrics of habitat quality that are understood to limit the species’ range. We focus on this species because pikas have been long studied as a model for biogeographic response to climate. They also have strong habitat associations and a recent history of population declines attributed to climate change.

The American pika, hereafter “pika,” depends on the relatively cool and mesic habitats created by fractured rock habitats, such as talus, boulder fields, and lava beds (Hafner 1994; Henry, Henry, and Russello 2012; A. T. Smith, Nagy, and Millar 2016; Figure 1). Such habitats are typically produced at higher elevations by freeze–thaw processes, and pikas are generally associated with permafrost and periglacial features (Hafner 1994; C. I. Millar and Westfall 2010). Range retractions (Beever et al. 2016; Jeffress, van Gunst, and Millar 2017) and local extirpations (Nichols, Klingler, and Peacock 2016; Jeffress, van Gunst, and Millar 2017; Stewart, Wright, and Heckman 2017) have been documented for pikas, with climate implicated as a driver of most losses (Hafner 1994; Beever, Brussard, and Berger 2003; Grayson 2005; Beever et al. 2010, 2011; Wilkening et al. 2011; Stewart et al. 2015; Wilkening, Ray, and Varner 2015; Stewart, Wright, and Heckman 2017; Rodhouse et al. 2018). Several range projections suggest dramatic losses in response to climate change during this century (Galbreath, Hafner, and Zamudio 2009; Calkins et al. 2012; Stewart et al. 2015). However, climate vulnerability predictions for pikas in eight national parks suggest a more complex future in which some populations might persist into the next century, whereas others are likely to be extirpated (Schwalm et al. 2016). Likewise, other analyses have shown that “ecoregion” is the most important predictor of pika range associations with climate (A. B. Smith et al. 2019), supporting predictions that pika response to climate change may be regionally variable and idiosyncratic (Jeffress et al. 2013). Furthermore, subsurface microclimates can extend the range and facilitate local persistence of this species in otherwise hostile environments (Varner and Dearing 2016; Ray, Beever, and Rodhouse 2016; Camp et al. 2020), which compounds regional variability between population responses and complicates vulnerability assessments based on surface climate projections.

Figure 1. An American pika (Ochotona princeps) in the Rocky Mountains. Photograph courtesy of Kristi Odom.

A growing body of work has demonstrated a close relationship between rock-ice features (RIFs) and pika distribution (Hafner 1994; C. I. Millar and Westfall 2010; C. I. Millar, Westfall, and Delany 2014; Wilkening, Ray, and Varner 2015), suggesting that active RIFs may indicate better habitat quality for pikas. This association is supported by evidence that pikas cannot tolerate temperatures above 25°C without access to cooler temperatures (MacArthur and Wang 1974; A. T. Smith 1974; Hafner 1993, 1994; Henry, Henry, and Russello 2012; A. T. Smith, Nagy, and Millar 2016). Biogeographic patterns also support the association of pikas and RIFs: a survey of pika sign (current or old) at 420 sites in the Southwestern United States found that the majority of current occurrences were on RIFs (C. I. Millar and Westfall 2010). As recently as 1994, all pika populations in the southern Rocky Mountains were estimated to be within 5 km of current permafrost (Hafner 1994). RIFs are also associated with subsurface microclimates that are both more moderate than ambient temperatures (cooler in the summer and warmer in the winter) and less variable year-round (C. I. Millar and Westfall 2010; C. I. Millar, Westfall, and Delany 2013; Wilkening, Ray, and Varner 2015). In the first study to measure pika FGM levels as a function of habitat type, pikas inhabiting rockslides associated with RIFs exhibited lower FGM levels than those in an adjacent area without RIFs (Wilkening, Ray, and Varner 2015). Therefore, we consider sites with subsurface ice to be of higher habitat quality than sites lacking RIFs.

FGMs are metabolized stress-associated hormones present in feces that can be used as a measure of glucocorticoids in the body (Keay et al. 2006; Homyack 2010; Dantzer et al. 2014). FGMs represent a time-averaged metric of circulating glucocorticoids (Sheriff et al. 2011; Dantzer et al. 2014), which is dependent on gut passage time and varies by species (Palme et al. 1996). For pikas, the gut passage time over which FGMs are averaged is about 12.5 hours (Wilkening, Ray, and Sweazea 2013). FGMs are therefore used as a metric of relatively “chronic” stress (Keay et al. 2006; Dantzer et al. 2014), as opposed to metrics of “acute” stress like corticosterone in blood, which can vary within minutes of exposure to a stressor (Romero and Reed 2005; Sheriff et al. 2011; Ellis, McWhorter, and Maron 2012). Previous studies have confirmed that FGM measurements provide a biologically relevant measure of physiological stress in pikas (Wilkening, Ray, and Sweazea 2013; Wilkening, Ray, and Varner 2015; Wilkening and Ray 2016). Furthermore, pikas that exhibited higher FGM levels were less likely to survive the following year in the wild (Wilkening and Ray 2016), pointing to the utility of FGMs as an indicator of individual fitness. After appropriate validation, Wilkening, Ray, and Varner (2015) found that FGM levels varied between two sites, one with and one without subsurface ice. As such, we expect that FGM levels will vary with habitat quality across our study sites. The use of FGMs as a potential indicator of habitat quality could therefore provide timely information about environmental changes that may be missed by studies of range retraction and/or population declines that have already occurred.

Here, we focus on FGM levels to address questions about habitat quality for pikas occupying RIF and non-RIF habitats in Boulder and Larimer Counties, Colorado. Specifically, we investigated whether pikas living in close association with RIFs (presumably higher quality habitat) appear to experience lower levels of FGM compared to those living in habitats without RIFs (lower quality habitat). Expanding on Wilkening, Ray, and Varner (2015), we sampled pika scat from two types of RIFs: “active” rock glaciers thought to harbor subsurface ice recently and “fossil” rock glaciers considered long devoid of subsurface ice (as classified by Janke 2005, 2007). To understand other influences of habitat quality on FGMs, we also measured fine-scale habitat differences between active and fossil rock glaciers, including potential metrics of thermal stress (microclimate) and available resources (talus depth, clast size, and land cover metrics related to forage). To account for possible seasonal variation in pika FGMs (Whipple et al. 2021), scat samples were collected in the alpine spring and fall. We then modeled FGM as a function of habitat-driven differences (rock glacier status or fine-scale habitat metrics) and season after accounting for effects of individual sex, or as a function of temperature-driven differences between habitat types. We predicted that pikas living within recently active rock glaciers would exhibit lower FGM levels than pikas living within fossil rock glaciers. Developing FGM measurements as an assay of habitat quality would help pinpoint important habitat features on the landscape that are otherwise difficult to identify, which could ultimately help refine conservation efforts and inform predictions of future pika distributions.

Materials and methods

Study system

The American pika is a small (~150 g) mammal that lives in rock piles (talus slopes, boulder fields, and lava beds) in Western North America (A. T. Smith and Weston 1990). Pikas defend relatively small territories (14–34 m in diameter; A. T. Smith and Weston 1990) and typically do not relocate after establishing a territory. Furthermore, pikas do not migrate or hibernate, so each individual is exposed year-round to the quality of habitat contained within its territory. Each territory owner usually builds a winter food cache or “haypile” consisting largely of graminoids (such as grasses and sedges) and forbs (herbaceous flowering plants). Haypiles are amassed during the growing season and consumed over winter (Dearing 1997). Because pikas are central-place foragers, haypiles are typically located near the centroid of a pika’s territory (Huntly, Smith, and Ivins 1986; Roach, Huntly, and Inouye 2001). We determined territory locations by the presence of fresh haypiles in fall or by the presence of fresh defecations at old haypile sites in spring. Although an individual sometimes constructs multiple distinct caches within its territory, each haypile sampled in this study was separated by at least 50 m from other haypiles, which exceeds the typical distance between territory centers (A. T. Smith and Weston 1990). Thus, our samples likely derive from distinct individuals.

Fecal pellets were collected in 2018 from pika territories located in rock glaciers within eight sites along the Front Range of Colorado: four in Rocky Mountain National Park (ROMO; 40°19′N, 105°39′W) and four adjacent to the Niwot Ridge Long-Term Ecological Research site (NWT; 40°3′N, 105°36′W). To identify the status of the rock glaciers at each site, we used geographic data layers created by Janke (2005, 2007). Janke (2005, 2007) used a combination of photographs, field investigations, and previous studies of terminus (leading-edge) motion to classify rock glaciers along the Front Range of Colorado as “active,” “inactive,” or “fossil.” Active rock glaciers were classified as moving, exhibiting steep frontal (leading-edge) slopes and pronounced ridges and furrows due to the presence of internal ice (Janke 2005, 2007). Both inactive and fossil rock glaciers were classified by the absence of movement, indicated by gentler frontal slopes. However, inactive rock glaciers may still contain ice, whereas fossil rock glaciers are thought to contain no ice (Whalley and Martin 1992; Barsch 1996). Fossil rock glaciers were further indicated by increasing cover of vegetation within the boundaries of a formerly active rock glacier (Janke 2005, 2007). We chose to study only active and fossil rock glaciers (omitting those that were inactive) because these represented the most extreme differences in recent conditions. Active and fossil rock glaciers were paired in the same drainage, when possible, to help control for spatially varying environments and conditions. Within each geographic area (ROMO or NWT), we selected two active and two fossil rock glaciers (Figure 2). Collectively, these study sites occurred at a range of elevations and aspects (Table 1).

Figure 2. Study site locations within Rocky Mountain National Park (Larimer County) and adjacent to Niwot Ridge (Boulder County), Colorado, USA. Blue triangles = active rock glaciers, red triangles = fossil rock glaciers. Black boundary delimits the southeastern portion of Rocky Mountain National Park.

Table 1. Paired active (A) and fossil (F) rock glaciers studied in Rocky Mountain National Park (ROMO) and adjacent to the Niwot Ridge long-term ecological research site (NWT), Boulder County, Colorado, USA.

Site name	Abbreviation	Location	Rock glacier status	Elevation (m)	Predominant aspect	Total number of territories sampled
Tyndall Gorge	A1	ROMO	Active	3,534	East	10
Lake Haiyaha	F1	ROMO	Fossil	3,122	East	15
Sky Pond	A2	ROMO	Active	3,342	East	10
Sky Pond	F2	ROMO	Fossil	3,321	South	13
Lake Isabelle	A3	NWT	Active	3,560	Southeast	10
Lake Isabelle	F3	NWT	Fossil	3,356	South	10
Green Lake 5	A4	NWT	Active	3,620	North	12
Green Lake 1	F4	NWT	Fossil	3,438	South	20

At each of the eight sites, we collected fresh (perched, sensu Stewart and Wright 2012) fecal pellets from at least ten unique territories, sampling each territory in both spring and fall. Fecal pellet sampling stations (hereafter, “stations”) were established by searching systematically for haypiles in the talus, starting at the terminus/toe (lower margin) and working toward the head (upper margin) of the talus slope. When a conspicuous haypile was found, we searched within a 5-m radius for a scat station where the territory owner appeared to deposit pellets regularly. The station location was recorded with a handheld Global Positioning System, and the search for new haypiles would continue after moving 50 m laterally or 50 m up the fall line (steepest slope) of near the edge of the talus.

Sample collection

FGM levels were assayed from fecal pellets collected noninvasively from free-ranging individuals in natural habitats. Pikas, like all lagomorphs, create two types of feces. Cecal feces are the soft portion of waste that has passed through the gut only once. After re-ingestion and further digestion, cecal material is re-excreted as firm fecal pellets (A. T. Smith and Weston 1990; Nichols 2010). Fecal pellets change color and texture as they age, and fresh pellets (less than one week since deposition) are greenish in color and fibrous inside (Nichols 2010). We collected only fresh (greenish, fibrous) fecal pellets from adult pikas (>2.75 mm in diameter) because FGM levels can vary between adults and juveniles (Wilkening, Ray, and Sweazea 2013). We targeted green pellets that were “perched” precariously (freshly) on rocks rather than pellets of less determinate age that had collected in depressions or on soil (Stewart and Wright 2012). Pellets start to lose their green color after one week (Nichols 2010) and remain perched (glued to the rock by urine) only until wind, rain, or a passing animal dislodges them. These points give us confidence that the pellets we collected were less than a week old (Nichols 2010). We note that samples could be a collection of multiple defecation events.

Postdeposition effects of climate should be considered when collecting feces from the field (Sheriff et al. 2011; Dantzer et al. 2014). In a study where Snowshoe hare (Lepus americanus) fecal pellets were experimentally exposed to environmental conditions for six days (all combinations of warm or cool and wet or dry), FGM levels varied under all treatments when compared to the control; however, FGM levels varied the least when exposed to cool and dry conditions (Lafferty et al. 2019). Similarly, a study that exposed American pika pellets to ambient conditions in pika latrine sites in several different ecoregions across the species range found that FGM levels varied inversely with precipitation and directly with temperature during the month prior to sample collection (Wilkening, Ray, and Varner 2016). However, these effects were only evident when FGM levels were compared across different ecoregions (e.g., between the Rocky Mountains and the Oregon Cascades or lava bed habitats) and not when samples were compared within the southern Rocky Mountain ecoregion (Wilkening, Ray, and Varner 2016). The southern Rocky Mountain ecoregion, where our study sites are located, has a relatively cool and dry climate, and because our sites are located within 23 km of each other in the same ecoregion, we do not expect postdepositional climate to have an effect on FGM levels in our study.

At each site (n = 8), the coordinates of multiple collection stations were established in the spring and revisited in the fall. We collected pellets in 2018 from 12 June through 16 July (“spring” in the alpine) and 9 September through 29 October (“fall”). After the collection of fresh pellets in the spring, any remaining pellets were cleared from each station to facilitate identification of fresh pellets in the fall. Samples were stored in paper coin envelopes labeled with sampling date and coordinates. To avoid contamination between samples, we scooped up pellets using the envelope flap or pushed pellets into the envelope using a stick or rock. We collected ten to fifteen pellets (0.1 g dry mass total) per sample to facilitate FGM analyses (Wilkening, Ray, and Sweazea 2013) and sex determination assays (described below). Samples were stored at −20°C within 24 hours of sampling, where they remained until analysis.

Sex determination

Because FGM levels can vary between sexes (Wilkening, Ray, and Sweazea 2013), we attempted to use a portion of each sample for pika sex determination through genetic analysis. When samples contained surplus pellets following the FGM assay, these surplus pellets were used to provide data on pika sex by season at most sampling stations.

Genetic samples from ROMO and NWT were analyzed in two different laboratories using slightly different assays for sex identification; however, both protocols relied on established, validated methods and included appropriate controls (Castillo et al. 2014; Lamb, Robson, and Russello 2014). Detailed polymerase chain reaction (PCR) recipes and conditions can be found in Supplement S1. In both analyses, sex was determined based on consensus across all PCR replicates to account for allelic dropout and PCR errors resulting from low-quality and low-quantity DNA template. Samples were identified as male, female, or unknown according to the following criteria: male = SRY fragments visible in two of three replicates or microsatellite and SRY fragments present in the same replicate; female = microsatellite fragments were present in at least one of the three replicates and SRY fragments were never present; unknown = only SRY fragments were present in one replicate, target genes never amplified, or samples were too small to subset for genetic analyses.

Habitat surveys

In 2019, we recorded pika-associated habitat metrics within a 12-m radius of each station to understand local habitat characteristics and fine-scale land cover types. Surveys were conducted between 11 July and 2 September, to coincide with the timing of highest vegetation cover in the alpine. We modeled our survey protocols after methods developed in previous studies (Rodhouse et al. 2010; Erb, Ray, and Guralnick 2011; Jeffress et al. 2013; Ray, Beever, and Rodhouse 2016). Recorded habitat metrics were described in Whipple et al. (2021) and included talus depth (visual estimate of the depth, in increments of 0.5 m, of the deepest crevice observed within 12 m of scat station), rock size (longest axis of largest rock; LR₁), tenth-largest rock (LR₁₀; within the 12-m radius) and Daubenmire’s (1959) estimates of percentage land cover for six cover types (rock, bare ground, graminoid, forb, shrub, and tree) within the 12-m radius. Specifically, the forb category included all herbaceous plants; graminoids included all grasses, rushes, and sedges; and shrub included all woody plants. We then used these data to calculate additional metrics of habitat, including mean rock size (0.5[LR₁ + LR₁₀]), rock size variation ([LR₁ − LR₁₀]/[0.5(LR₁ + LR₁₀)]), and several vegetation metrics that could affect diet estimated from percentage cover types: graminoid-to-forb ratio, forb cover, grass cover, and total available forage (sum of graminoid, forb, and shrub cover).

Temperature data collection

At each of the eight sites (Table 1), one HOBO temperature-sensing data logger (hereafter, “sensor”) was placed under the surface of the talus to record subsurface temperatures. All eight sensors (models H08 pro, H08, or H02; Onset Computer Corporation, Bourne, MA) recorded temperature (±0.2°C) every 1 or 2 hours for one year, starting in September 2018. Sensors were placed at similar depths (~0.75 m) near the center of each rock glacier. Sensor location and elevation were recorded with a handheld Global Positioning System unit and photos were taken to aid in sensor retrieval. A generalized additive model was used to estimate mean, maximum, and minimum daily temperatures from the interval data recorded by each sensor. Generalized additive models were fitted using the mgcv package (Wood 2011) in R 3.5.1 (R Core Team 2018).

Hormone analysis

To measure FGM level (concentration) in fecal pellet samples, we used a commercially available Corticosterone Enzyme Immunoassay kit (cat. No. K014-H5; Arbor Assay Design, Inc., Ann Arbor, MI) that has been previously validated for O. princeps (Wilkening, Ray, and Sweazea 2013). We followed the manufacturer’s protocols to extract FGMs from dry scat as described in Whipple et al. (2021). For each sample, a random subset of pellets was dried using a DNA SpeedVac for 30 minutes and then weighed to 0.1 ± 0.02 g. Next, 200-proof ethanol was added in proportion to sample mass and pellets were smashed using a small metal spatula in a vial. These smashed pellets in ethanol were shaken for 30 minutes and then centrifuged to extract the supernatant, which was transferred to a new vial and concentrated using a DNA SpeedVac. Concentrated supernatant was stored at −20°C until all samples were prepared for the corticosterone assay. Then, a mix of assay buffer and ethanol was used to reconstitute the supernatant to its original volume (<5 percent ethanol content by volume) before loading onto enzyme immunoassay plates according to the manufacturer’s instructions. The commercial kit antibody is a sheep polyclonal antibody specific for corticosterone with a cross-reactivity of 100 percent. The cross-reactivity of the antibody with other steroids is very low—0.24 percent, 0.03 percent, and <0.08 percent for progesterone, testosterone, and estradiol, respectively—so we do not expect reproductive events to have affected our results. We assayed samples, standards (5,000, 2,500, 1,250, 625, 312.50, 156.25, and 78.125 pg/mL), and controls (maximum binding and nonspecific binding) in triplicate according to the kit instructions. Final fecal FGM concentrations were calculated from a standard curve and reported as pg FGM/g dried feces. Interassay coefficients of variation (%CV = CV × 100) were 3.6 percent and 13.9 percent for the nonspecific and maximum binding controls, respectively, and 8.9 percent for pooled fecal extracts. Intra-assay coefficients of variation were between 5.0 percent and 11.6 percent for nine plates (n = two controls, seven standards, two pooled extracts, and 21 samples per plate).

Linear mixed-effects models

FGM concentration was log-transformed and used as the response variable in linear mixed-effects models to characterize variation in FGM levels among and within rock glacier habitats. We evaluated two sets of candidate models, each based on up to five fixed effects (model sets A and B, Table 2), plus random effects of territory nested within site to account for repeated measures and location effects using our full data set (n = 184).

Table 2. Fixed effects and associated hypotheses used to construct three sets of linear mixed-effects models of glucocorticoid metabolite concentrations measured from fecal sampling in individual territories of the American pika.

Model set	Fixed effect(s)	Predictor category	Hypothesized mechanism
A	Rock glacier status (active or fossil) elevation (m), location	Broad-scale habitat	Spatial trends in the environment that are generally shared among several territories within a site
B	Graminoid-to-forb ratio, % forb cover, % graminoid cover, graminoid + forb + shrub cover, rock size (LR1, LR10), mean rock size, rock size variation, deepest crevice	Fine-scale habitat	Spatial trends in the environment that often vary among several territories within a site
C	No. of days below −5°C, no. of days below −10°C	Subsurface temperature	Acute exposure to cold temperature
C	No. of days above 10°C, no. of days above 15°C	Subsurface temperature	Acute exposure to hot temperature
C	Average early-winter temperature (October–December), average late-winter temperature (April–June), early-winter temperature variation (October–December), late-winter temperature variation (April–June)	Subsurface temperature	Chronic exposure to cold temperature
A, B, C	Sex	Inter-individual	Expected sex-associated variation in FGM levels
A, B, C	Season	Temporal	Temporal trends shared among sites

FGM = fecal glucocorticoid metabolite.

In the first set of models (model set A), we explored qualitative effects of broad-scale habitat metrics and season on FGMs, because stress-associated hormones like FGMs are known to vary seasonally in many animals (Baker, Gobush, and Vynne 2013), including pikas (Whipple et al. 2021). FGM level was modeled using fixed and interacting effects of season and rock glacier status (active or fossil) or using fixed and interacting effects of season and elevation in meters (mean centered and divided by standard deviation). We expected FGM levels to be higher in pikas living in fossil rock glacier habitats. Alternatively, elevation was used as a proxy for the many factors that vary with altitude, including rock glacier status. To account for effects of latitude and land management, we also included a fixed effect of location (ROMO or NWT; Figure 2) and its interaction with season. Finally, we included a fixed effect of sex along with all broad-scale habitat metrics, to account for the expected variation in FGM levels between sexes (Wilkening, Ray, and Sweazea 2013).

In the second set of candidate models (model set B), we modeled FGM level as a linear response to fine-scale habitat differences. To account for potential differences in fine-scale habitat differences between active and fossil rock glaciers, FGM level was modeled in response to either (1) a single fixed effect of one of our fine-scale habitat metrics plus fixed effects of season and/or sex or (2) fixed and interacting effects of one fine-scale habitat metric and rock glacier status plus fixed effects of season and/or sex. Fine-scale habitat metrics represented differences among territories, such as talus characteristics (talus depth, largest rock size, tenth-largest rock size, mean rock size, rock size variation) and vegetation metrics (graminoid-to-forb ratio, forb cover, graminoid cover, and the total available forage). Talus depth (coded as a categorical variable in 0.5-m increments) was explored because talus may act as an insulator from surface conditions (Wilkening et al. 2011; Henry, Henry, and Russello 2012; C. I. Millar, Westfall, and Delany 2014, 2016; Varner and Dearing 2016; Wilkening, Ray, and Varner 2015), and deeper talus may provide a better refuge during extreme heat (Moyer-Horner et al. 2015). Rock characteristics were explored because we lack prior knowledge of how structural complexity of the habitat might affect habitat quality. Vegetation metrics were explored because vegetation could affect pika diet and because these metrics have been predictive of pika occupancy in other studies (e.g., Wilkening et al. 2011; Wilkening, Cole, and Beever 2019).

We also conducted an exploratory analysis to investigate the effects of subsurface temperature on FGM levels. In this analysis, we evaluated three sets of candidate models (broad- and fine-scale habitat metrics along with subsurface temperature metrics) using a subset of our data (n = 170). Again, each model set was based on up to five fixed effects (Table 2) plus random effects of territory nested within rock glacier. Covariates remained the same for broad- and fine-scale habitat model sets as above. For the set of subsurface temperature models (model set C), we modeled FGM level as a linear response to each of eight potential fixed effects, each a metric of thermal stress in pikas, based on previous studies (Beever et al. 2011; Wilkening et al. 2011; Varner et al. 2015; Wilkening, Ray, and Varner 2015). These metrics were intended to represent acute exposure to cold temperatures (number of days below −5°C or −10°C, 23 October 2018 to 19 July 2019), acute exposure to hot temperatures (number of days above 10°C or 15°C, 23 October 2018 to 19 July 2019), and chronic exposure to cold temperatures (mean and variation of temperatures during the two months in early and late winter, 23 October 2018 to 23 December 2018 and 1 April 2019 to 1 June 2019, respectively; Table 2). Our metric of chronic cold exposure represents a period of time when snowpack is variable in our study area, which could expose subsurface habitats to acutely cold temperatures. Once a thick blanket of snow covers the talus, subsurface temperatures stabilize near 0°C, providing a stable refuge for pikas during subfreezing surface temperatures. Temperature variation (%CV) was calculated from the standard deviation of subsurface temperatures recorded during the specified period divided by the mean of those same temperature measurements. To account for potential differences in subsurface temperatures between active and fossil rock glaciers, FGM level was modeled in response to either (1) a single fixed effect of one of our temperature metrics plus fixed effects of season and/or sex or (2) fixed and interacting effects of one temperature metric and rock glacier status plus fixed effects of season and/or sex. The sample size was reduced in this model set because the single temperature logger in one of the active rock glaciers malfunctioned and failed to record subsurface temperature. We also used Student’s t-test to determine whether the means of each temperature metric differed between active and fossil rock glaciers. Because we performed eight t-tests (for eight temperature metrics), we used a Bonferroni correction to reduce the chance of obtaining false-positive results, adjusting our alpha value to .00625 (0.05/8) to determine whether means were significantly different.

Collinearity among predictor variables was calculated using Pearson’s rho (continuous variables), Spearman’s rank correlation coefficient, or Kendall’s tau, and pairs with a high correlation (>0.7) were not included in the same model. Models were fitted using the lme4 package (Bates et al. 2015) in R 3.5.1 (R Core Team 2018). We ranked models using Akaike’s information criterion adjusted for a small sample size (AICc) (Burnham and Anderson 2002). Relative support for each model was calculated using ΔAICc for each model i within a set of candidate models (set A, B, or C), defined as AICc_i − AICc_min, where AICc_min was the lowest AICc in the model set. Following convention, we considered models with ΔAICc > 2 to have lower support than the “top” model (ΔAICc = 0). We evaluated the relative likelihood of each model using Akaike weights, w[i] = exp(0.5*dAICc[i])/sum[j](dAICc[i = j]) (Burnham and Anderson 2002). Additionally, we evaluated pseudo-R² to estimate the variance explained by the fixed effects only (marginal R²) and the variance explained by both the fixed effects and random effects (conditional R²) compared to the null model with only random effects using the r.squaredGLMM function in R package MuMIn (Barton 2022). In each model set, we included a null model of FGM as a function of random effects (territory nested within site) and no fixed effects, as well as models based on the fixed and interacting effects described above and listed in Table 2. Model set A included fixed and interacting effects of season and broad-scale habitat metrics (rock glacier status, elevation, location) plus a fixed effect of sex. Model set B included fixed and interacting effects of rock glacier status and fine-scale habitat metrics (graminoid-to-forb ratio, percentage forb cover, percentage graminoid cover, graminoid + forb + shrub cover, rock size [LR₁, LR₁₀], mean rock size, rock size variation, deepest crevice) plus fixed effects of season and/or sex. Model set C included fixed and interacting effects of rock glacier status and a measure of thermal stress in pikas (acute exposure to cold temperature, acute exposure to hot temperature, and chronic exposure to cold temperature) plus fixed effects of season and/or sex.

Results

During the spring and fall of 2018, 190 fecal samples were collected from 100 sampling stations in eight rock glacier sites along the Front Range of Colorado. However, our model sample size was reduced to 184 fecal samples from 95 sampling stations due to persistent snow hindering our ability to conduct habitat surveys when we returned to stations in 2019. Of the 184 useful samples, 77 were collected from the four active rock glaciers (39 in the spring, 38 in the fall) and 107 were collected from the four fossil rock glaciers (55 in the spring, 52 in the fall). The elevation of sampled territories ranged from 3,122 m to 3,703 m (Figure 3).

Figure 3. Elevation range of sampled territories grouped by site (n = 8). Abbreviated site names distinguish active (“A”) from fossil (“F”) rock glacier sites (Table 1).

Fecal pellet mass was sufficient to allow genetic analysis for 168 of 184 samples to determine the sex of the putative territory owner. Of the samples collected in the spring (n = 94), 36 were deposited by males, 18 were deposited by females, and 40 were deposited by pikas of unknown sex due to autosomal DNA amplification failure (n = 29) or insufficient sample mass (n = 11). Of the samples collected in the fall (n = 90), 50 were deposited by males, 15 were deposited by females, and 25 were deposited by pikas of unknown sex due to DNA amplification failure (n = 20) or insufficient sample mass (n = 5). In model comparisons, we categorized sex either as a three-class variable (n(sex) = 86 male, 33 female, 65 unknown; Figure 4) or a four-class variable (n(sex) = 86 male, 33 female, 49 PCR failure, 16 insufficient sample mass). In both cases, there was no significant difference in log(FGM) levels among the sex categories (analysis of variance, F_(2,181) = 0.118, p = .889; Figure 4).

Figure 4. Pika fecal glucocorticoid metabolite (FGM) concentration, log-scaled by sex and season. (a) In the spring, eighteen samples were classified as female (pink circles) and thirty-six were classified as male (blue circles); the remaining forty samples were classified as unknown due to DNA amplification failure (black circles) or insufficient size (black triangles). (b) In the fall, fifteen samples were classified as female, fifty as male, and the remaining twenty-five as unknown.

Because we determined sex from samples collected at a number of territories both in the spring and fall, we were also able to estimate turnover in the sex of putative territory owners. We found that of the eighteen territories occupied by females in the spring, nine were occupied by males in the fall and six were occupied by pikas of unknown sex. Of the thirty-six territories occupied by males in the spring, six were occupied by females in the fall and fourteen were occupied by pikas of unknown sex.

When we compared models in sets A and B using the full data set (n = 184), the best-supported model included a positive interaction of elevation and season (Figure 5). Two other models were also well supported (ΔAICc < 2), one containing a positive interaction of graminoid-to-forb ratio and rock glacier status plus a negative effect of season, followed by a model containing a positive interaction of percentage graminoid cover and rock glacier status, plus a negative effect of season (Figure 6). A model with a negative interaction of rock glacier status and season (Figure 7) and a model with similar terms to those of the top model with the addition of a sex covariate were moderately supported (∆AICc = 2.65 and 3.69, respectively; Table 3). Location (ROMO or NWT) garnered less support (∆AICc = 4.74, Table 3). A model containing a negative effect of season alone had very weak support (ΔAICc > 10). Including sex as a covariate did not improve model rank, regardles of whether “sex” was categorized using three or four categories, as described above. Akaike weights indicated that the interaction between standardized elevation and season was the strongest predictor of FGM levels (w[i] = 0.30), followed by the interaction of graminoid to forb ratio and rock glacier status plus season (w[i] = 0.26). Additionally, one of the highest conditional R² values was associated with model 1A, indicating that its fixed and random effects estimated the most variance in FGM compared to other models. Generally, marginal and conditional R² values were similar across models and the addition of our random effects always increased model fit.

Figure 5. Fecal glucocorticoid metabolite (FGM) concentration by elevation and season. Each point represents a unique fecal sample colored by the season in which it was collected (green = spring, orange = fall). Lines show mean FGM concentration by season with 95 percent confidence intervals.

Figure 6. Fecal glucocorticoid metabolite (FGM) concentration by percentage graminoid cover (left) or graminoid-to-forb ratio (right) colored by rock glacier status (blue = active, red = fossil). Lines show mean FGM concentration with 95 percent confidence intervals.

Figure 7. Fecal glucocorticoid metabolite (FGM) concentration for all sampled territories grouped by rock glacier status and season. When available, fecal samples were collected from each territory in the spring and in the fall of 2018 from locations along the Front Range of Colorado. These box-and-whisker plots depict median (heavy horizontal line), interquartile range (box), range (whiskers), and outliers (dots).

Table 3. Relative support for models of fecal glucocorticoid metabolite (FGM) concentration (log-FGM) in pika territories sampled within Rocky Mountain National Park and adjacent to the Niwot Ridge long-term ecological research site (n = 184).

Model rank (no., set)	Model fixed effects	Effect sign	AICc	ΔAICc	Akaike weights	k	Marginal R^2	Conditional R^2
	Null model		171.29	54.90	0	4
1, A	Standardized elevation	(−)	116.44	0	0.30	7	0.23	0.57
	Season (fall)	(−)
	Standardized elevation: season (fall)	(+)
2, B	Graminoid-to-forb ratio	(−)	116.74	0.31	0.26	8	0.35	0.50
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	Graminoid-to-forb ratio: rock glacier status (fossil)	(+)
3, B	Percentage cover graminoid	(−)	118.28	1.84	0.12	8	0.35	0.50
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	Percentage cover graminoid: rock glacier status (fossil)	(+)
4, A	Rock glacier status (fossil)	(+)	119.08	2.65	0.08	7	0.33	0.52
	Season (fall)	(−)
	Rock glacier status (fossil): season (fall)	(−)
5, A	Sex (female = 0, male = 1, unknown = 2)	(var)	120.12	3.69	0.05	9	0.24	0.57
	Standardized elevation	(−)
	Season (fall)	(−)
	Standardized elevation: season (fall)	(+)
6, A	Rock glacier status (fossil)	(+)	120.67	4.23	0.04	6	0.32	0.50
	Season (fall)	(−)
7, B	Graminoid-to-forb ratio	(−)	120.70	4.23	0.04	10	0.35	0.50
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	Sex (female =0, male =1, unknown = 2)	(+)
	Graminoid-to-forb ratio: rock glacier status (fossil)	(+)
8, A	Management area (NWT)	(+)	121.18	4.74	0.03	8	0.33	0.52
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	Rock glacier status (fossil): season (fall)	(−)
9, B	Tenth-largest rock	(+)	121.96	5.53	0.02	8	0.33	0.50
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	Tenth-largest rock: rock glacier status (fossil)	(−)
10, B	Percentage cover graminoid	(−)	122.54	6.10	0.01	10	0.35	0.50
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	Sex (female = 0, male = 1, unknown = 2)	var
	Graminoid: rock glacier status (fossil)	(+)
11, A	Sex (female = 0, male = 1, unknown = 2)	var	122.86	6.43	0.01	9	0.33	0.52
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	Rock glacier status: season (fall)	(−)

Models are ranked (no.) in order of increasing AICc, and Akaike weights provide a measure of relative support for each model. Only models with ΔAICc <7 and the null model are shown. Top models (in bold) have equivalent support (ΔAICc < 2). Pseudo-R² estimates the variance explained by the fixed effects only (marginal R²) and the variance explained by both the fixed effects and random effects (conditional R²) compared to the null model with only random effects. The effect sign associated with each respective covariate is reported as positive (+), negative (−), or variable (var). For each model, sample size (n) was 184 and number of fitted parameters (k) is reported. Models in set C were based on a lower sample size and were not included in this comparison.

As an exploratory analysis, we compared models in sets A and B with models in set C. In this exploratory analysis, sample size was reduced from 184 (in results presented above) to 170 due to temperature sensor failure in one of the active rock glaciers. Top model rankings were similar to those in Table 3 for broad-scale habitat metrics (model set A); however, fine-scale habitat metrics (model set B) were outcompeted by several subsurface temperature models (model set C). The best-supported model included a positive interaction of elevation and season, followed by a model including a positive interaction of the number of days where subsurface temperature reached above 10°C and rock glacier status and a negative interaction of the number of days where subsurface temperature reached above 15°C and rock glacier status (Table 4). Both fixed effects of the number of days that subsurface temperature reached above 10°C or above 15°C were positively associated with FGM levels. However, only the mean of one temperature metric was significantly different between rock glacier status (Student’s t-tests; Figure 8). Specifically, in late winter (April–June), average temperature was warmer in fossil rock glaciers than in active rock glaciers (t = −6.73, p = .00623). Akaike weights indicated that the interaction between standardized elevation and season was the strongest predictor of FGM levels (model rank 1A; w[i] = 0.14). Additionally, the highest conditional R² value was associated with model 1A, indicating that its fixed and random effects estimated the most variance in FGM compared to all other models. Marginal and conditional R² values were similar across all other models and the addition of our random effects always increased model fit.

Figure 8. Subsurface temperature metrics grouped by rock glacier status (blue = active, red = fossil). Metrics were estimated from one subsurface temperature logger placed in the center of each rock glacier site. Standard box-and-whisker plots depict median (heavy horizontal line), interquartile range (box), range (whiskers), and outliers (dots). P-values are reported for each Student’s t-test comparing means of subsurface temperature metrics between rock glacier status. Using Bonferroni correction, significance is determined by an alpha level of .00625.

Table 4. Relative support for models of fecal glucocorticoid metabolite (FGM) concentration (log-FGM) from an exploratory comparison of all model sets (Table 2) using data from all samples with corresponding subsurface temperature data (n =170).

Model rank	Model fixed effects	Effect sign	AICc	ΔAICc	Akaike weights	k	Marginal R^2	Conditional R^2
	Null		157.81	54.79	0	4
1, A	Standardized elevation	(−)	103.02	0	0.14	7	0.24	0.57
	Season (fall)	(−)
	Standardized elevation: season (fall)	(+)
2, C	No. of days above 10°C	(+)	103.99	0.96	0.09	8	0.35	0.52
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	No. of days above 10°C: rock glacier status (fossil)	(+)
3, C	No. of days above 15°C	(+)	104.50	1.48	0.07	8	0.35	0.52
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	No. of days above 15°C: rock glacier status (fossil)	(−)
4, A	Rock glacier status (fossil)	(+)	104.54	1.52	0.07	6	0.32	0.52
	Season (fall)	(−)
5, C	Average early-winter temperature (October–December)	(−)	105.07	2.05	0.05	8	0.35	0.51
	Rock glacier status (fossil)	(+)
	Season (fall)	(−)
	Average early-winter temperature (October–December): rock glacier status (fossil)	(+)

Models are ranked (no.) in order of increasing AICc, and Akaike weights provide a measure of relative support for each model. Only models with ΔAICc ≤ 2 and the null model are shown. Pseudo-R² estimates the variance explained by the fixed effects only (marginal R²) and the variance explained by both the fixed effects and random effects (conditional R²) compared to the null model with only random effects. The effect sign associated with each respective covariate is reported as positive (+), negative (−), or variable (var). For each model, sample size (n) was 170 and number of fitted parameters (k) is reported.

Discussion

We found that several broad- and fine-scale habitat metrics could explain variation in the FGM levels of wild American pikas, suggesting an association between habitat quality and individual stress levels. This study corroborates a previous investigation of FGM response to putative habitat quality (Wilkening, Ray, and Varner 2015) and extends those results by providing data from two seasons and from a much larger and well-replicated spatial extent. Further, our results add to a growing body of work that illustrates the relationship between broad-scale habitat quality, season, and glucocorticoid (GC) measurements. For example, capercaillie (Tetrao urogallus), wood mice (Apodemus sylvaticus), and arctic ground squirrels (Urocitellus parryii) have been shown to exhibit higher GC levels in poorer quality habitats (Coppes et al. 2018; Navarro-Castilla and Barja 2019; Sheriff et al. 2017, respectively). Additionally, the seasonal trend in FGM levels found in this study has been mirrored in studies of arctic ground squirrels and snowshoe hares where GC levels were elevated at the end of winter and lowered through the active season (Sheriff et al. 2012; Lavergne et al. 2020, respectively).

Although our noninvasive study was not designed to measure physiological indices (e.g., body weight, immune response, glucose response) that could help shed light on the relationship between fitness and FGM level (Breuner, Delehanty, and Boonstra 2013; Dantzer et al. 2014), a previous study found that pikas exhibiting higher FGM levels were less likely to survive the following year after capture (Wilkening, Ray, and Varner 2016), suggesting a direct link between elevated FGM levels and reduced fitness. We recommend further studies to evaluate the relationship between FGM levels and physiological indices to suggest mechanisms governing pika response to habitat quality.

Broad-scale habitat metrics (model set A)

Variation in FGM level was best explained by the interaction of elevation and season. These results indicate that pikas living at higher elevations have lower average FGM levels and that FGM levels are higher in the spring than in the fall (Figure 5). Pika occupancy surveys have documented population losses at lower elevations in several parts of the species range (Beever et al. 2016; Jeffress, van Gunst, and Millar 2017; Stewart, Wright, and Heckman 2017), and our study suggests that these range retractions might be due to processes that cause higher FGM levels. Furthermore, our results suggest that FGM levels might be useful for predicting populations at risk, given that we detected higher FGM levels in populations at lower elevations. However, it is important to note that elevation is confounded with rock glacier status: active rock glaciers are commonly found at higher elevations and fossil rock glaciers are more likely to occur at lower elevations in our study area (Janke 2007; Figure 3). We used elevation as a proxy for broad spatial trends in the environment that are generally shared among several territories within a site. We propose that abiotic and biotic factors could be driving seasonal trends across elevations, such as the timing of snowmelt or plant green-up. For example, a thick blanket of snow at higher elevations can buffer subsurface talus temperatures from extremely cold above-surface temperatures. In the spring, when cold snaps are common in the mountains, lower (subalpine) elevations may not have adequate snow cover to shield pikas from acute cold stress. Likewise, the timing of plant green-up occurs sooner at lower elevations, which could influence the timing of plant senescence later in the summer, when pikas are most actively harvesting food for the winter (Dearing 1997). Thus, the timing of the period during which pikas can acquire highly nutritious food might be important in determining haypile quality, with implications for pika FGM levels in the following spring.

The interaction of rock glacier status and season was also supported as a predictor of pika FGM levels. This model supports our hypothesis that FGM levels are higher in fossil rock glaciers than in active rock glaciers and refines our understanding by showing that this effect on FGMs is more pronounced in the spring (Figure 7). The coherent difference in FGM levels between active and fossil rock glaciers indicates that habitat quality correlates with a stress metric in our system. Results from this study are consistent with previous research at NWT showing that pikas exhibit lower FGM levels in one area dominated by RIFs compared to habitat not associated with a RIF (Wilkening, Ray, and Varner 2015). This result held in our study when replicated several years later and across several RIF sites, indicating that RIFs may constitute higher quality habitats for pikas.

Further, all of the top models in our study indicate that season is an important driver of FGM levels. Specifically, FGM levels were highest in spring, suggesting that the previous winter might have been stressful for the pikas in our study and that overwinter stress might generally be a limiting factor for pikas. Temperature-based metrics of cold stress have been shown to affect pika FGM levels, occupancy, and persistence in other studies (Beever, Brussard, and Berger 2003; Beever et al. 2010, 2011; Wilkening et al. 2011; Jeffress et al. 2013; Wilkening, Ray, and Varner 2015; Rodhouse, Hovland, and Jeffress 2017). We hypothesize that the reduction in snow cover in North America’s Mountain West (Mote et al. 2005) along with the loss of subsurface ice in our research area (Leopold et al. 2015) might indicate changes in the microhabitats that pikas depend on for thermally buffered refugia. Without these refugia, pikas may be faced with a more stressful and unpredictable thermal environment. Our results are also consistent with previous findings that snow cover and winter cold stress are important determinants of pika distribution and extinction risk (Beever et al. 2010; Erb, Ray, and Guralnick 2011; Guralnick, Erb, and Ray 2012; Schwalm et al. 2016; Rodhouse, Hovland, and Jeffress 2017) and add increasing evidence that pikas experience more stress in the cold season than in the warm season. These results were mirrored in a similar study on a nearby population of pikas at NWT, where FGM levels were highest in the spring when measured every two weeks throughout the growing season (Whipple et al. 2021).

If our methods were replicated and elevated FGM levels in the spring stayed consistant across years and in different regions, it might suggest that localized pika losses are due to sustained stressors related to habitat quality, rather than stochastic events like epizootics or extreme weather events. Elevated FGM levels in the spring may be a product of the cumulative effects of winter cold stress combined with the lack of forage availability at a time when pika haypiles are depleted and decayed (Dearing 1997). However, it is important to note that high FGM levels in the spring could also be attributed to the timing of reproduction. Reproduction likely elevates FGM levels in females (Dantzer et al. 2014) and possibly also in males due to increased reproductive competition as evidenced in arctic ground squirrels (Urocitellus parryii; Delehanty and Boonstra 2011) and yellow-bellied marmots (Marmota flaviventris; J. E. Smith et al. 2012). Parturition usually coincides with snowmelt (J. S. Millar 1972; A. T. Smith 1978), which is when our spring sampling effort occurred.

Fine-scale habitat metrics (model set B)

Our results provide evidence that fine-scale (territory-level) habitat metrics also influence FGM levels in pikas. Other well-supported models include an interaction of rock glacier status and either graminoid-to-forb ratio or percentage graminoid cover. These models suggest that higher graminoid cover is correlated with reduced FGM levels in pikas, at least in active rock glacier habitats (Figure 6). In contrast, studies in the warmer and drier Great Basin have demonstrated that pika occupancy is correlated with higher forb cover and reduced graminoid cover (Wilkening et al. 2011). Although these results may seem contradictory (forb cover being associated with higher occupancy in one region but higher stress in another), we suggest that this difference may be due to regional differences in pika water balance between the hotter and drier Great Basin and the cooler and moister Rocky Mountains. Studies in both the Great Basin (Elliot & Flinders 1984; Beck and Peek 2005) and our study region (Bhattacharyya and Ray 2015; J. Varner, unpublished data) demonstrate that forbs are consistently higher in water content, protein content, and micronutrient content than graminoids in both ecoregions. In the Great Basin, higher graminoid cover occurred at more xeric sites (Wilkening et al. 2011). If pikas get nearly all of their water from the vegetation they eat (Krear 1965; J. A. Smith and Erb 2013), then pikas living in the warmer and drier portions of the Great Basin ecoregion—where graminoid cover is higher and pikas likely lose more water to evaporation—may experience higher water stress. In contrast, the graminoids in our study (including primarily Carex spp., Deschampsia cespitosa, and Juncus spp.) are not necessarily indicative of more xeric sites or reduced water availability in the Rocky Mountain ecoregion, and pikas may be less water-limited on a diet of graminoids at these sites. Pikas in this region tend to eat graminoids immediately, whereas forbs are often stored in haypiles for later consumption during winter (Dearing 1996). As a result, pikas living in territories with higher graminoid cover may have more access to immediate energy needed to carry out daily activities, especially in the spring when energy reserves are low and haypiles are depleted. Finally, it is important to note that higher graminoid cover was an important determinant of FGM levels only at active rock glaciers (Figure 6), where vegetation was concentrated along the talus periphery, as were most of our sampling stations. We propose that because pikas appeared to be nearer to meadow interfaces at active rock glaciers sites, compared to fossil rock glaciers sites where pikas could be found further away from meadow interfaces (more insular territories), the easy access to forage that is apparently adequately nutritious could be driving lower FGM levels at active rock glacier sites.

Inter-individual metrics

The addition of a sex covariate to our top models did not further explain the variation in FGM levels, nor was sex well supported as the sole fixed effect in a model (ΔAICc > 10). Additionally, we found no significant difference in FGM levels among the sex categories.

By determining the sex of the samples collected from the same territory in both seasons, our results suggest turnover in the sex of territory owners, in contrast to previous studies that suggest that pika territories are likely to retain an owner of the same sex through time (Krear 1965; Sharp 1973; Tapper 1973; A. T. Smith and Ivins 1983, 1984). Of the fifty-four territories in which sex could be determined from spring samples, sex determined from fall samples differed in fifteen cases. Turnover rates in site occupancy are highly dependent on local influences but can be as high as 50 percent in some populations (Rodhouse et al. 2018). The changes in sex we observed between samples collected at the same territory could be a result of a change in latrine use or territory intrusion by neighbors and not necessarily a change in the owner of the territory. Additionally, it is possible that multiple pikas contribute pellets at a single latrine and therefore that we could have collected pellets from two different individuals. This is a particularly important consideration, given that our genetic analysis does not distinguish individuals within a sample but rather the presence or absence of the male-specific SRY locus. Thus, if any pellets in a sample came from a male, the whole sample could be definitively considered male, which could also explain why territories more frequently turned from female to male than vice versa. To improve our understanding of pika latrine use and the role of latrines in territorial behaviors, behavioral studies of marked individuals combined with genotyping of marked and unmarked individuals could be conducted at known stations to confirm the sex of pikas who deposit pellets.

Temperature metrics (model set C)

We used our preliminary temperature measurements to explore differences in subsurface habitat between rock glacier status, which we hypothesized could drive habitat quality for pikas. The top models from our exploratory analysis revealed that measures of acute exposure to hot temperature (the number of days in which subsurface temperature reached above 10°C or the number of days in which subsurface temperature reached above 15°C) were associated with higher FGM levels. These results are consistent with previous pika occupancy studies that observed population declines in warmer habitats (Wilkening et al. 2011; Stewart, Wright, and Heckman 2017). However, it is important to note that only one of our metrics differed between rock glacier status (active versus fossil; Figure 8). Specifically, fossil rock glaciers experienced warmer temperatures in late winter compared to active rock glaciers (Figure 8). Active rock glaciers may have been colder because they tended to be higher in elevation than fossil rock glaciers and could have been experiencing colder ambient temperatures during late winter. We acknowledge that our limited sampling may not have fully characterized the subsurface temperature within these rock glacier habitats, especially because we are reporting data from only one sensor per rock glacier and for only nine months of one year (October–July). Due to the timing of our fieldwork, we were unable to collect subsurface temperature data during late summer, which is the warmest time in the alpine and a possible driver of chronic heat stress in pikas. Our temperature sensors were also placed no more than 0.75 m below the rock surface; thus, it is also likely that such shallow temperature readings did not capture the full range of temperatures that a pika could sample in this habitat. Previous geophysical studies on rock glacier A4 from our study detected ice 4 to 5 m below the surface (Leopold et al. 2011), depths that are likely accessible to pikas but often too difficult to outfit with sensors (Ray, Beever, and Rodhouse 2016).

Broader implications

Results from this study showed that FGM levels in populations of a small mammal were correlated not only with season but also with habitat metrics (elevation, vegetation, and rock glacier status). The pattern in FGM levels that we observed across landscapes indicates that stress-associated hormones could be used to evaluate seasonal variation in habitat quality for pikas. In particular, seasonal variation in FGM levels suggests that spring may be more physically demanding for pikas, a trend also found in another small population of pikas on NWT (Whipple et al. 2021). Further, poor quality habitats seemed to amplify the effects of seasonal events, especially in spring. Climate change is expected to cause more extreme events, such as increased wildfires, more extreme temperatures, and shifts in precipitation. For species living in extreme environments that experience regular perturbations, adding new stressors to the mix may overtax the individual to the point that they can no longer cope (Wingfield, Kelley, and Angelier 2011). This could have huge conservation implications for climatically sensitive species that live in extreme environments. Female arctic ground squirrels (Urocitellus parryii) living in habitats with short growing seasons exhibited higher fecal corticosterone metabolite (FCM) levels and produced offspring with higher FCM levels than females living in habitats with longer growing seasons (Sheriff et al. 2017). However, females living in habitats with short growing seasons were larger and in better body condition, probably driven by selection of robust phenotypes that could better cope with the poor habitat conditions (Sheriff et al. 2017). Similarly, alpine chipmunks (Tamias alpinus) living in less climatically suitable habitats, as indicated by ecological niche models, and in areas that experience higher daily maximum temperatures exhibited higher FGM levels, indicating that their environment may be challenging their physiological limits (Hammond, Palme, and Lacey 2018).

Given the effects of season we observed, management actions, such as habitat modifications, that consider season may provide the most benefit for pika conservation. Likewise, management plans for other climate-sensitive species like caribou (Rangifer tarandus granti) and polar bears (Ursus maritimus) have recommended protection of seasonal habitats as a critical component for conservation (Western Arctic Caribou Herd Working Group 2019; U.S. Fish and Wildlife Service 2016, respectively). As conservationists and managers begin to take actions for the 30×30 initiative (preserving 30 percent of the planet’s land and water by 2030), we suggest that alpine habitats containing active rock glaciers be considered for future conservation areas. Not only are active rock glaciers important microhabitats for pikas but they could also similarly provide key climate refugia for other mountain-dwelling animals and plants (C. I. Millar et al. 2015). Further research is needed to characterize current rock glacier status and distribution, as well as recent changes to this distribution, across Western North America. This would also be an invaluable resource for future conservation of mountain ecosystems.

During the summer, the ice within rock glaciers can continue to provide water to surrounding vegetation, creating dense and nutritious plant communities (C. I. Millar and Westfall 2010; C. I. Millar, Westfall, and Delany 2013). In the Sierra Nevada range and Great Basin, this dense vegetation contains a higher diversity of shrubs and herbaceous plants, compared to surrounding alpine and subalpine habitats not associated with a rock glacier (C. I. Millar and Westfall 2010; C. I. Millar, Westfall, and Delany 2013). Additionally, the prolonged water resource can help surrounding plants survive during drought conditions and could provide biotic refugia for arthropods during warmer and drier climates (C. I. Millar and Westfall 2010; C. I. Millar et al. 2015).

Despite their importance to the broader ecosystem, ice features such as rock glaciers and permafrost are declining as a result of warmer temperatures in summer and fall (Leopold et al. 2015). Estimates suggest that if current average temperatures were to increase by 2°C, ~94 percent of permafrost would be lost along the Front Range of Colorado (Janke 2005). However, decline in ice will likely respond slowly to warming temperatures and possibly even slower in rock glacier habitats where the ice is protected by a rocky carapace (Janke 2005; C. I. Millar and Westfall 2008; but see Monk and Ray 2022). For pikas, which some estimates predict could be extirpated from some of our study sites by 2100 (Schwalm et al. 2016), active rock glaciers may contain overlooked microrefugia that could help the species persist and should be considered in future species distribution modeling research.

Acknowledgments

We thank the following individuals for aiding in field and/or lab efforts: Elliot Martin, Nadiv Edelstein, Micaela Seaver, Claire Atkins, Ryan Agabani, Airy Adriana Peralta, Angie Rodriguez, Troy Tanner, and Brad Schrom. We also thank Andrew Martin, Steven Schmidt, and Christopher Lowry for the use of their lab space and equipment. The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the U.S. Fish and Wildlife Service.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are openly available in the Environmental Data Initiative repository at https://doi.org/10.6073/pasta/929e2f61cb0926cf58ebac09d187148f, and https://doi.org/10.6073/pasta/b15009b35926522f75e43e1e9a861d64.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15230430.2022.2146633

Additional information

Funding

This work was supported by the Niwot Ridge LTER program under NSF DEB-1637686, a Boulder County Nature Association research grant, a CU Boulder Ecology and Evolutionary Biology Department research grant, a Warren Wilson College Pugh Grant, and an Indian Peaks Wilderness Alliance Paddon/Gellhorn Wilderness Research grant.