Like mother, like daughter: heritability of female Richardson's ground squirrel Urocitellus richardsonii cortisol stress responses.

Activation of the hypothalamic-pituitary-adrenal (HPA) axis liberates glucocorticoids, which provides an acute indication of an individual's response to stressors. The heritability of the stress response in wild mammals, however, remains poorly documented. We quantified the cortisol stress response of female Richardson's ground squirrels (RGSs) to handling and physical restraint, testing for: (1) the effects of individual age, time of day, and sample latency; (2) repeatability within individuals; (3) narrow-sense heritability; and (4) differences among individuals owing to potential genetic and/or environmental effects. We detected a positive linear relationship between baseline plasma cortisol (BL-cortisol) concentration and stress-induced plasma cortisol (SI-cortisol) concentration that defined each individual's cortisol stress response. BL-cortisol, SI-cortisol, and stress response did not differ according to the time the sample was taken, or by subject age. Cortisol stress response was highly repeatable within individuals, had a mother-offspring heritability of h 2 = 0.40 ± 0.24 (mean ± SE), full-sibling heritability of h FS 2   =   0.37   ±   0.71 , and half-sibling heritability of h HS 2 =   0.75   ±   1.41 . Stress responses of sibling groups, immediate-family groups, and squirrels within a given area did not differ, whereas those of individuals from more distantly related matrilines did. Our results highlight the natural variability in HPA axis reactivity among individuals by quantifying both BL- and SI-cortisol levels, demonstrate partial heritability of the stress response that is not attributable to environmental variation, and suggest that at least part of an individual's stress response can be accounted for by differences in matrilineal history.

The acute stress response promotes behavioral and physiological changes that facilitate the maintenance of homeostasis in the face of challenging conditions, and is manifested in part via the activation of the hypothalamic–pituitary–adrenal (HPA) axis triggering the release of glucocorticoids such as corticosterone and cortisol into circulation (Barton 2002; Moore and Jessop 2003; Romero 2004; Palme et al. 2005; Herman et al. 2012). These “stress hormones” occur at baseline levels that vary seasonally with changes in breeding status and mating (Coe and Levine 1995; Kitaysky et al. 1999; Boonstra et al. 2001; Nunes et al. 2006; Chauke et al. 2011). They also vary over the course of the day in diurnal species, with peak concentrations upon waking, and declining levels throughout the day (Weitzman et al. 1971; Coe and Levine 1995; Tsang et al. 2014). Against that background, circulating glucocorticoid concentrations vary over the course of minutes (Sheriff et al. 2011; Hare et al. 2014; Jenkins et al. 2014) in response to acute stressors associated with the presence of a predator (Eilam et al. 1999; Cockrem and Silverin 2002; Chauke et al. 2011), conspecifics (Boonstra et al. 2001; Nunes et al. 2006; Yue et al. 2006), or novel and stressful conditions such as handling (Boonstra et al. 2001; Blas et al. 2007; Atwell et al. 2012; Bennett et al. 2012; Hare et al. 2014).

The magnitude of baseline glucocorticoid concentration varies relative to resource availability and body condition, behavioral traits, and the age of the individual. Reduced food availability and thus diminished body condition is associated with decreased baseline glucocorticoids and fitness according to the Cort-Fitness hypothesis [see Bonier et al. (2009) for a review; Kitaysky et al. 1999; Bókony et al. 2009]. Pregnant or laying mothers experiencing low resource availability have heightened baseline glucocorticoid levels (Love et al. 2005, 2009; Monclús et al. 2011), particularly under the threat of predation (Sheriff et al. 2010a). Baseline glucocorticoid levels may reflect prevailing environmental conditions (e.g., resource availability and predation risk) or internal processes (behavioral syndrome and age-related allostatic loading) rather than representing an adaptive response per se (Elliott et al. 2014). Higher baseline glucocorticoid levels are also associated with certain behavioral syndromes (Atwell et al. 2012; Montiglio et al. 2012; Clary et al. 2014; Dosmann et al. 2014) that enhance reproductive success but decrease survival, such as boldness (Smith and Blumstein 2008), and the age and reproductive state of the individual (Reeder and Kramer 2005; Hämäläinen et al. 2015; Homberger et al. 2015), with older and more reproductively inactive individuals exhibiting exaggerated stress responses and longer recovery times to return to baseline levels (Reeder and Kramer 2005; Elliott et al. 2014; Homberger et al. 2015). The exaggerated stress response and decrease in recovery latency of older individuals is likely not adaptive, but rather the result of increased allostatic load owing to cumulative wear and tear throughout an organism’s life (McEwen and Wingfield 2003; Romero et al. 2009). Further, increased glucocorticoids during embryogenesis may drastically affect offspring development, senescence, and oxidative stress (Haussmann et al. 2012).

Elevated glucocorticoid levels and stress in mothers can also have a lasting effect on offspring through pre- and post-natal programming [see Love et al. (2013) for a review; Sheriff et al. 2010a] and thus offspring fitness. Maternal glucocorticoids influence offspring sex (Cameron et al. 2008; Ryan et al. 2011), growth rate (Dantzer et al. 2013), neuromotor development (Schneider and Coe 1993; DiPietro 2012; Giesbrecht et al. 2015), and stress physiology. For example, high maternal glucocorticoids at parturition can alter the offspring stress response (i.e., the change in plasma cortisol concentration from pre- to post-stressor levels; Sheriff et al. 2010b; Homberger et al. 2015) by weakening internal negative feedback mechanisms that normally limit the offspring stress response (e.g., Kanitz et al. 2003; Kapoor et al. 2006). Thus, offspring of highly stressed mothers may be more physiologically responsive to external stressors (Kanitz et al. 2003; Kapoor et al. 2006; Love et al. 2013), although increased post-natal care can attenuate the stress response (Bókony et al. 2009; McGhee and Bell 2014). These developmental differences may serve to adapt offspring to stressful environmental conditions such as high population density (Dantzer et al. 2013), high predation risk (Boonstra et al. 1998; Sheriff et al. 2010a; Vitousek et al. 2014), or unpredictable resource availability (Homberger et al. 2015). The offspring stress response is thus plastic in response to the mother’s environmental and internal stressors relative to offspring baseline glucocorticoids. Therefore, the magnitude of the stress response may experience increased selection compared to baseline glucocorticoid concentrations alone (Bókony et al. 2009; Rensel and Schoech 2011; Jenkins et al. 2014; Homberger et al. 2015) and thus may better predict the resulting behavior and fitness of wild individuals than baseline glucocorticoids (Vitousek et al. 2014; Homberger et al. 2015).

Despite the importance of maternal glucocorticoid levels in determining offspring stress response, few studies have examined whether the stress response is repeatable over time, the extent to which it is heritable, and its susceptibility to environmental effects, such as differences in resource availability and predation risk. While baseline glucocorticoid measures are likely most repeatable when using integrative measures obtained from urine, hair, or feces (Bosson et al. 2009), the stress response is most repeatable when considering point-based methods (Rensel and Schoech 2011; Cook et al. 2012; Grace and Anderson 2014; Small and Schoech 2015). However, some studies have reported low repeatability for fecal glucocorticoids (an integrative measure) in wild Sciurid populations (Smith et al. 2012; Dantzer et al. 2016). Baugh et al. (2014) also reported that point-based measures of baseline and stress-induced glucocorticoids were not repeatable in great tits Parus major, suggesting that local environmental stressors are the most important determinant of these levels. While both integrative and point-based measures may or may not be repeatable, both are affected by the methods employed in the study, with more reliable measures obtained within a small temporal window or with captive or lab-reared individuals (Smith et al. 2012; Grace and Anderson 2014; Boulton et al. 2015). Intra-individual repeatability and inter-individual differences indicate that the stress response is somewhat fixed and thus is susceptible to natural selection (Smith et al. 2012), contingent on the trait’s heritability (Jenkins et al. 2014; Homberger et al. 2015).

Heritability in the narrow sense (h2) is the proportion of phenotypic variance attributable to additive genetic variance, independent of environmental factors, and is measured as h2 on a continuum from 0 (not heritable) to 1 (entirely heritable). Baseline and stress-induced plasma glucocorticoid concentrations are at least partly heritable in some mammals, fish, and birds (h2 = 0.15–0.62; Pottinger and Carrick 1999; Tanck et al. 2001; Bartels et al. 2003; Jenkins et al. 2014); however, few studies have explored the heritability of the stress response in wild animals. Fewer still have quantified the narrow-sense heritability of the glucocorticoid stress response in mammals. These findings highlight the need for studies of the stress response among free-living animals, thereby encompassing effects of environmental variation that are not present in a controlled lab setting (Jenkins et al. 2014).

Due to the acute nature of the stress response (Sheriff et al. 2011), its potential heritability, and its correlation with biological fitness, we used 2 point-based measures (baseline and stress-induced concentrations) of total plasma cortisol to assess the stress response in Urocitellus richardsonii. Richardson’s ground squirrels (RGSs) are a moderately social, diurnal species (Michener 1979b, 1981, 1983) that have been the subject of several studies documenting the repercussions of variation in both fecal and plasma glucocorticoids (e.g., Ryan et al. 2011, 2014; Clary et al. 2014; Hare et al. 2014). While total circulating cortisol (both free and bound circulating cortisol) includes biologically inactive cortisol and not just free circulating cortisol (Breuner et al. 2013), we have previously characterized the RGS stress response in terms of total cortisol (Ryan et al. 2011; Hare et al. 2014) and this appears to provide a good measure of underlying stress physiology in this species. Individual RGSs experience high predation pressure owing to a broad suite of terrestrial and avian predators (Michener 1979a), which imposes selection pressure on the stress response. Additionally, matrilineal relatedness data are easily collected on populations through live trapping around the natal burrows at juvenile emergence within the first weeks of a newborn’s life (Michener 1985). These attributes, in concert with a recently validated glucocorticoid assay (Hare et al. 2014) and variation in resource availability and predation pressure across the study area, render this species ideal for studying the stress response and its intergenerational transmission.

We characterized the impacts of age, diel variation, sample retrieval latency, and individual repeatability on the plasma cortisol stress response of female RGSs. We also assessed the heritability of the stress response by comparing littermates raised in similar environments to estimate narrow-sense heritability coefficients. Finally, we applied linear models contrasting sibling, immediate-family, distant/matrilineal-family, and individuals residing in the same versus different spatial areas within a colony to elucidate the factors contributing to differences in the cortisol stress response among these groups. We predict that the RGS stress response will be repeatable within individuals and heritable among familial groups if there is a genetic basis for variation in the stress response.

Materials and Methods

Marking and tracking squirrels

Research was conducted on a free-living population of U. richardsonii at the Assiniboine Park Zoo (49°52′N, 97°14′W) in Winnipeg, Manitoba, Canada, where squirrels have been marked individually since 2003. National or Tomahawk live traps (Tomahawk Live Trap Co., Hazelhurst, WI, USA) baited with No Name® peanut butter (Loblaw Companies Ltd., Brampton, ON, Canada) were used to capture ground squirrels. All squirrels were ear-tagged (National Band and Tag Co., Monel #1, Newport, KY, USA) at first capture for permanent identification and marked for rapid visual identification using a unique pattern of hair dye applied to their dorsal pelage (Clairol Hydrience, No. 52 Pearl Black, Stamford, CT, USA). All juveniles were live trapped and marked as they emerged from their natal burrows (mid-May through mid-June in Manitoba) allowing accurate records of matrilines to be maintained from year to year (Michener 1985).

Blood sampling

Blood samples were drawn following the protocol outlined in Hare et al. (2014). Sampling occurred between 18 and 28 June 2013 and 17 June through 29 July 2014, always subsequent to weaning of juveniles by their dam. Approximately 16% of the colony was subject to blood sampling. We targeted individuals who had not been trapped at least 1 h before blood sampling to preclude sampling individuals who had premature activation of the stress axis. BL-cortisol samples were taken from the medial saphenous vein of either the right or left hind leg using a 25 gauge × 5/8” heparinized needle and 1 mL syringe. We refer to the baseline sample latency (BSL) time as the time interval from the closure of the live trap with the squirrel inside to the end of the BL-cortisol blood sample draw (mean ± SE: 4.5 min ± 12 s). The BSL of all BL-cortisol samples was between 2 and 7 min post-capture, and sampling was done in a natural field setting to minimize premature HPA axis activation in response to handling (Hare et al. 2014). Squirrels were then transported to the on-site veterinary hospital—an unnatural and presumably stressful situation—in cloth-covered live traps where “stress-induced” (SI-cortisol) samples were drawn between 9 and 29 min post-capture (i.e., stress-induced sample latency time, SSL, 15 min ± 1 s), or 6–26 min following the BL-cortisol sampling (SSL–BSL, mean ± SE: 11 min ± 29 s); the time interval from the closure of the live trap with the squirrel inside to the end of the stress-induced blood sample draw. After sampling, direct pressure was applied to the sampled vein for roughly 1 min to ensure bleeding had ceased. Squirrels were then returned to and released at their point of capture, and not targeted for blood sampling again for at least 24 h if they had no visual signs of hematoma, or at least 48 h if they had developed a hematoma. Trapping effort was relatively uniformly distributed throughout the field season and avoided areas we knew to have been previously disturbed by natural predators and/or humans on the same day of blood sampling, thereby minimizing any confounding effect on BL-cortisol.

Whole blood was transferred to a 2-mL microcentrifuge tube and centrifuged for 3 min at 13,000 × g. The resulting plasma was pipetted and stored in 2-mL polycarbonate vials marked with a unique identifier referencing the date and time of collection, as well as the donor squirrel’s ear tag number. Samples were stored in liquid nitrogen and subsequently transferred to a −80 °C freezer until cortisol extractions were performed. All methods employed were approved under protocol F12-014 of the University of Manitoba, Fort Garry Campus, Protocol Management and Review Committee, by the Research Review Panel of the Assiniboine Park Zoo, and by Manitoba Conservation under Wildlife Scientific Permit WB14952.

Quantification of cortisol

The radioimmunoassay (RIA; Sheriff et al. 2011) method validated by Hare et al. (2014) for U. richardsonii was used to measure plasma concentrations of total cortisol. Upon thawing the frozen samples, plasma was diluted 1/100 by adding 10 µL of plasma to 990 µL of ice-cold 95% ethanol and vortexing for 5 s. Of this 1/100 solution, 100 µL was then centrifuged at 13,000 × g for 5 min at 4 °C. After decanting the resulting supernatant into a new tube, the pellet was resuspended in 500 µL of 95% ethanol through vortexing. The resuspended pellet was centrifuged as described and the resulting supernatant was combined with the first and placed in a sample concentrator to evaporate the ethanol (Savant SpeedVac; Thermo Scientific, Waltham, MA, USA). The extract was then stored at −20 °C until the day of the assay, which typically occurred within 1 week of sample extraction.

Samples were reconstituted before measurement by RIA by vortexing in 250 µL of buffer consisting of 0.1 M phosphate, 0.9% NaCl [w/v], and 0.5% bovine serum albumin [w/v]. To perform the cortisol RIA, 100 µL of cortisol-specific antibody (1:9,000 dilution; Fitzgerald Industries, North Acton, MA, USA) was combined in an assay tube with 100 µL of tritiated cortisol (5,000 ± 250 disintegrations per min per 100 µL; GE Healthcare, Piscataway, NJ, USA) and either 100 µL of sample or a known amount of unlabeled cortisol (Steraloids, Newport, RI, USA). Tubes were then incubated at room temperature for 1 h and cooled at 4 °C overnight until the reaction was terminated by adding 100 μL of dextran-coated (0.5% w/v) charcoal (5% w/v). Tubes were then placed on ice for 15 min and centrifuged at 2,500 × g for 30 min at 4 °C. The supernatant was then decanted into a 7-mL scintillation vial, with 4 mL of Ultima Gold scintillation fluid (Perkin Elmer, Waltham, MA, USA). The radioactivity of tritiated cortisol in each tube was measured using a liquid scintillation counter (TriCarb; Perkin Elmer, Waltham, MA, USA). Known concentrations were used to interpolate sample concentrations using GraphPad Prism software (GraphPad Software Inc., La Jolla, CA, USA). Our RIA’s extraction efficiency, minimum detectable levels, inter-assay variation, and intra-assay variation were all comparable to those values previously reported for RGS plasma (Ryan et al. 2011; Hare et al. 2014).

Sample sizes

CPR collected and quantified the total BL- and SI-cortisol in plasma samples from 19 different female squirrels (N = 19 BL- and SI-cortisol samples) in the summer/fall of 2013, whereas 34 different females (N = 46 BL- and SI-cortisol samples, N = 11 SI-cortisol samples only) were obtained for total cortisol by KBN in the summer/fall of 2014. Overall, N = 65 BL-cortisol and N = 76 SI-cortisol samples were acquired, with 61 of the individual squirrels having paired BL- and SI-cortisol measurements, as well as some repeated samples within an individual (N = 17 repeated BL- and SI-cortisol). BL- and/or SI-cortisol for 19 breeding females (1 age-4, 3 age-3, 3 age-2, and 12 yearlings) in June 2013 (n2013 = 19), 27 breeding females (1 age-4, 1 age-3, 8 age-2, and 17 yearlings) in June/July 2014 and 11 juveniles in June/July 2014 (n2014 = 38) were obtained, with repeat samples between years (2013 and 2014) for 4 individuals. Thus, we sampled a total of N = 53 different females in our study. However, 15 individuals had 2 or more samples taken of both BL- and SI-cortisol, 4 of these 15 being re-sampled in both 2013 and 2014 and the remainder (N = 11) being sampled twice or more in 2014. Additionally, 2 individuals were repeat-sampled in 2014 alone for only SI-cortisol, making the total N = 17 for repeated SI-cortisol samples. There were N = 51 stress response measures (i.e., with paired BL- and SI-cortisol measures) from 45 different individuals, and upon averaging repeated samples and removing outliers 45 measures remained. Stress response from 20 mothers and 29 daughters (1.45 offspring/mother on average) were compared to calculate narrow-sense heritability, h2 using a parent–offspring regression. Full- and half-sibling heritability estimates of the stress response were also calculated among siblings for N = 15 daughters from 7 dams. Stress responses of individuals within 7 sibling groups (N = 15), 19 immediate-family groups (N = 44), 12 matrilineal-family groups (N = 41), and individuals from 5 different spatial locations (N = 45) were contrasted to explore the relative contributions of genetic/epigenetic, as well as environmental factors on the stress response.

Statistical analysis

A measure of proportional change in BL- and SI-cortisol concentrations was used as the “stress response” for further analyses including the calculation of repeatability and heritability. The proportional difference was used because of the log-linear and bivariate normal relationship observed when plotting BL- versus SI-cortisol concentrations of samples (Bókony et al. 2009). While many studies simply use SI-cortisol (e.g., Blas et al. 2007; Rensel and Schoech 2011; Vitousek et al. 2014; Homberger et al. 2015), we wanted to preserve the rise-over-run relationship of SI-cortisol by BL-cortisol to examine what drives these differences in the ratio of SI- to BL-cortisol, or in effect, what places individuals above or below the line of best fit between SI- and BL-cortisol measures. The individual’s “stress response” is thus given as: where relatively large values above the line of best fit (Figure 1) represent higher SI-cortisol levels than predicted by the BL-cortisol and relatively small values below the line of best fit represent lower-than-predicted SI-cortisol relative to the BL-cortisol level. The stress response variable satisfied all assumptions of parametric tests, and was used for further analyses of heritability and differences among genetic/environmental groups.

Figure 1

Linear relationship between log-BL-cortisol and log-SI-cortisol concentrations (in ng/mL; r = 0.24, N= 65) in RGS U. richardsonii at the Assiniboine Park Zoo (Winnipeg, Manitoba, Canada).

A datum was deemed an outlier and removed from analyses if it was less than the lower quartile −1.5 * IQR (interquartile range) or greater than the upper quartile +1.5 * IQR (Rousseeuw and Croux 1993). One outlier for stress response (2.040 > 1.671) and 3 outliers for log SSL (1800, 1900, 3145 > 1578.5 s) were thus removed (∼4% of data) and all other data were employed in subsequent analyses. All assumptions of parametric tests were satisfied upon log10 transformation of BL-cortisol, SI-cortisol, and SSL, whereas BSL and stress response did not require transformation (Shapiro–Wilks W ≤ 0.97, P ≥ 0.066; Bartlett’s F ≤ 1.04, P ≥ 0.40). All tests were 2-tailed and results were considered significant where P < 0.05. Measures of central tendency and dispersion are presented as mean ± SE and back-transformed when log-transformation was applied. The effects of sample latency time (BSL and SSL), age, sex, and hour of sampling on the 3 plasma cortisol measures (BL-cortisol, SI-cortisol, and stress response) were examined using analysis of variance (ANOVA) and simple linear regression models.

To assess the repeatability of the stress response, we looked at whether the differences between first and second samples of BL-cortisol, SI-cortisol, and stress response within the same individuals differed from 0 using Student’s paired t-tests. Additionally, we performed a simple linear regression with a fixed intercept at (0,0) of the first and second samples of stress response. The latter test was used as we assumed that the stress response must pass through or near the origin if stress response is ever equal or near to 0. Repeated measurements of cortisol concentrations for inter-year (2013 and 2014) and intra-year samples (2014 only) were analyzed together due to our limited sample size (BL-cortisol and stress response: N = 4 and 11, respectively; SI-cortisol: N = 4 and 13, respectively). Because initial and repeated samples quantifying the stress response of individuals within 2014 varied in latency between 0 and 12 days apart, we tested for any relationship between latency and the absolute difference between the initial and subsequent stress response of 12 individual squirrels using linear regression. Differences between the variances of first and second samples for BL- and SI-cortisol were also examined using the Pitman’s (1939) paired variances test. Repeated measures of first and second samples were averaged for both BL- and SI-cortisol, and these measures were used for all subsequent analyses of stress response involving heritability (h2) and differences among genetic/environmental groups.

Two separate methods were employed to explore the narrow-sense heritability of the stress response, the first using a parent–offspring regression method, the second using an ANOVA model comparing sibling groups. While parent–offspring regressions are commonly used and more powerful, the sibling-ANOVA method avoids problems of maternal environment being correlated with the offspring environment (Becker 1984; Brodie and Garland 1993). BL- and SI-cortisol samples were analyzed using 1-way ANOVAs to calculate maternal family full-sibling (h2FS) and half-sibling (h2HS) heritability coefficients (Becker 1984, pp. 47–54): where is the among-group variance, within is the within-group variance calculated from ANOVA table mean squares (MS) estimates among and within related-sibling groups (MSamong and MSwithin, respectively), and k1 is the effective number of offspring for every dam/mother when there are unequal numbers of progeny sampled among dams, calculated as k =  = 2.13 for our study, where D is the number of dams or mothers, n. is the total number of offspring, and k is the number of offspring from the i-th dam (Becker 1984, p. 51). Full-sibling heritability assumes offspring are full siblings fathered and mothered by the same 2 individuals (∼50% shared genetics) and includes additive genetic variance as well as the variance from non-genetic drivers, such as maternal effects (Brodie and Garland 1993). Half-sibling heritability is similar but assumes offspring are fathered by different males and thus are half-siblings sharing ∼25% of their genetics with fellow littermates (Becker 1984). An underestimate of full-sibling heritability is thus expected when multiple males sire offspring (thus, some offspring are half-siblings; Brodie and Garland 1993), as is the case in RGSs, where multiple insemination results in multiple paternity within litters (Hare et al. 2004). We employed both full-sibling and half-sibling heritabilities to obtain a range containing the sibling h2 of our study population, as U. richardsonii are likely to be somewhere in the middle of both h2FS and h2HS. Mother–offspring heritability (h) was calculated for the stress response using 2 times the slope of a linear regression of mother–offspring stress responses (Becker 1984, pp. 103–106), weighted by the number of offspring sampled for each dam. This method of assessing heritability assumes no changes in the stress response of an individual with age (Falconer 1981), which is unlikely to be the case given the previously discussed age-related variation in the activity of the HPA axis (Reeder and Kramer 2005; Hämäläinen et al. 2015). A Student’s paired-sample t-test was used to compare mother–daughter cortisol measures in terms of differences in magnitude.

Similarities in the stress response, allowing contributions of genetic, environmental, and other influences, were examined using 1-way ANOVAs for sibling groups (e.g., sisters), immediate-family groups of mothers and their respective offspring (i.e., within 1 generation of one another, e.g., mother and her offspring), as well as for broader matrilineal families (where individuals share matrilineal ancestry within 2–3 generations, e.g., mother and offspring, grandmother, great-grandmother, aunt, etc.). To examine the relative contribution of spatial area and environment on stress response, we examined the stress response of individuals occupying different areas within the colony (areas separated by physical barriers such as asphalt or gravel paths or physical distances exceeding 100 m, over which breeding females and their offspring seldom interact) using a 1-way ANOVA. RGSs are predicted to be affected by their environmental variance if their above-ground environments are affected by different predation pressures, vegetation quality, and conspecific density/competition (Michener 1979a; Davis 1984; Michener and Koeppl 1985), all of which apply to our study site. If spatial area and hence, the environment, have an effect on the stress response, differences among the stress responses of individuals from different areas are expected. However, differences among individuals from different spatial areas may also be explained by shared genes, as female RGSs often inherit part of their mother’s territory after 9–10 weeks (Michener 1981, 1985), and thus the heritability of the stress response confounds our measure of environmental variance when considered alone. To determine if 1-way ANOVAs were significant relative to random models, we produced 1,000 randomized bootstraps with replacement in R (Wickham, 2009; R v.3.2.3, R Core Team 2015) and applied the same ANOVA model to each bootstrap to determine the significance of the observed model relative to random models based on the resulting distribution of P-values (“permutation ANOVA” or “percentile bootstrap method” in Haddon 2001). If the observed P-values were outside of the 95% percentiles of the 1,000 random ANOVA models (i.e., less than the 2.5th percentile or greater than the 97.5th percentile), then the grouping factor (immediate-family groups, matrilineal-family groups, or groups by spatial area) was considered significantly different from random models. Post hoc pairwise comparisons among groups were analyzed using Tukey Honest Significant Difference (HSD) tests.

Results

Cortisol and sample latency

SI-cortisol was higher than BL-cortisol for most samples (Student’s paired t64 = 7.72, P < 0.0001; back-transformed mean BL-cortisol: 49.9 ± 1.1 ng/mL, mean SI-cortisol: 86.3 ± 1.1 ng/mL), with 52/65 paired samples having higher SI- than BL-cortisol. There was a positive relationship between BL- and SI-cortisol concentrations for all samples (Figure 1; simple linear regression F1,63 = 20.3, P < 0.0001, r2 = 0.24, β0 = 1.27, β1 = 0.41), which was used to compute the normally distributed stress response (mean stress response: 1.17 ± 0.023). This relationship remained strong and consistent even when considering averaged samples with outliers removed (F1,44 = 16.6, P = 0.0002, r2 = 0.27, β0 = 1.27, β1 = 0.41).

Mean BSL ± SE was 271 ± 12 s (N = 65), whereas the SSL was 918 ± 1 s (N = 54). BSL did not affect BL-cortisol (F1,63 = 0.48, P = 0.49) nor SI-cortisol (F1,63 = 1.15, P = 0.29), and SSL did not affect SI-cortisol (F1,49 = 2.72, P = 0.11) meaning that earlier or later sample collections within the 0–7 min for BL-cortisol and 9–29 min for SI-cortisol were similar in terms of their cortisol values and did not change significantly within these time intervals. The stress response measure was also unaffected by sample latency (BSL: F1,59 = 0.40, P = 0.53; SSL: F1,49 = 1.14, P = 0.29), or by differences in transport time between BSL and SSL samples (SSL–BSL: F1,59 = 0.38, P = 0.54).

Age and time of day effects

Age of the squirrel sampled (age groups: 0 = juveniles, 1 = yearlings, 2+ = adults) did not exhibit any difference in terms of BL-cortisol, SI-cortisol, or stress response (Table 1). BL- and SI-cortisol concentrations decreased throughout the day (Table 1). However, the hour of sample collection (i.e., time of day) did not affect the stress response measure (Table 1), and, therefore, we did not include time of day as a covariate in subsequent models involving stress response as the dependent variable.

Table 1

Differences in BL-cortisol, SI-cortisol, and the stress response (i.e., SI-/BL-cortisol) as a function of age, time of day, or sex in RGS U. richardsonii subjected to an acute capture and handling stress

Explanatory variable	[log] BL-cortisol	[log] SI-cortisol	Stress response
Age	ANOVA	ANOVA	ANOVA
(0 = juveniles, 1 = yearlings, 2 + = adults)	F2,38 = 1.34, N = 45, P = 0.27	F2,43 = 1.13, N = 46, P = 0.33	F2,38 = 0.19, N = 41, P = 0.83
Time of day	Linear regression	Linear regression	Linear regression
(hours past 00h00 CDT)	F1,60 = 7.00, N = 62, P = 0.010*, β0 = 2.39, β1 = −0.052, r2 = 0.10	F1,60 = 6.51, N = 62, P = 0.013*, β0 = 2.51, β1 = −0.041, r2 = 0.10	F1,56 = 0.20, N = 58, P = 0.66
Sexa	No data	Wilcoxon Rank Sums	No data
(male or female)		W = 455, N = 55, P < 0.0001*, females: 91.46 ±1.10 ng/mL, nfemales = 42, males: 59.29 ± 1.07 ng/mL, nmales = 13

Note: Test results (where data exist) are presented, with those yielding significant differences including additional information about effect sizes (β0, β1, and r2 or mean ± SE).

Sex is presented as supplementary material for this study, as it was not directly examined. P values set in bold font with an * denote statistical significance at P ≤ 0.05.

Repeatability

Neither BL-cortisol (Student’s paired t14 = 0.54, P = 0.60), SI-cortisol (t16 = 1.70, P = 0.11), nor stress response (t14 = 0.85, P = 0.41) differed significantly between first and second sampling for repeatedly sampled individuals. A fixed-intercept linear regression of repeated measures revealed a positive relationship between first and second measures of stress response (F1,14 = 361.02, P < 0.0001, r2 = 0.967). Further, the number of days between the first and second sampling did not affect stress response repeatability within years (F1,10 = 0.333, N = 12, P = 0.58, r2 = 0.032; 4 ± 1 days between first and second sampling). Intra-individual variance in BL-cortisol was greater in the first (variance = 1.22 ng/mL) relative to the second sampling (variance = 1.05 ng/mL; Pitman’s test t13 = 2.84, P = 0.014, N = 15), but first and second samples did not differ in terms of variability in SI-cortisol or stress response (t15 = 0.01, P = 0.99; t13 = 1.75, P = 0.10; respectively).

Heritability

Heritability was quantified to compare stress responses among related individuals by calculating narrow-sense heritability estimates among mothers and daughters as well as among siblings. General mother–offspring heritability of the stress response was estimated as: h2 = 0.40 ± 0.24 (slope = 0.20, intercept = 0.97, N = 20). Full-sibling heritability and half-sibling heritability were calculated as h2FS = 0.37 ± 0.71 and h2HS = 0.75 ± 1.41, respectively (ANOVA F6,8 = 1.49, N = 15, P = 0.29). These estimates of heritability assume that mothers and offspring shared a common environment (Becker 1984), and the constrained regression assumes no difference in the magnitude between BL- and SI-cortisol for mothers versus offspring, which was upheld for mother–offspring pairs tested (paired t-test t19 = 0.92, N = 20, P = 0.37).

Group differences

Close-familial comparisons (sibling, immediate family) using the percentile bootstrap method simulating ANOVA models to evaluate the observed model yielded no significant differences with regards to the female stress response (Table 2). The degree of relatedness for siblings and immediate-family members may vary between 0.25 and 0.5 on average depending on the extent of multiple mating in female RGSs (Hare et al. 2004). Differences among stress responses were detected only for matrilineal groups (Table 2, Figure 2). There were no differences in the stress responses of individuals from the same spatial location (k = 5 locations; Table 2).

Table 2

RGS U. richardsonii stress response differences among sibling, immediate-family, and matrilineal-family groups, as modeled by permutation ANOVA models employing 1,000 simulated ANOVA models from the data selected randomly with replacement for comparison to the observed ANOVA models for each

Comparison	Coefficient of relatedness (r)	Observed F statistic	df	Observed P-value	Bootstrapped P-value, 95% CI
Siblings	0.25–0.5	1.49	6, 8	0.29	(0.0135, 0.9710)
Immediate family	0.25–0.5	1.49	18, 25	0.17	(0.0273, 0.9796)
Matrilineal family	0.0625–0.5	2.98	11, 29	0.009*	(0.0219, 0.9683)
Spatial location	N/A	1.55	4, 40	0.21	(0.0244, 0.9724)

Note: An asterisk denotes observed model significance at α = lower bootstrapped confidence interval. Bold font indicates a comparison for which a statistically significant difference was detected, while a P-value with an * specifically denotes significance at P ≤ 0.05.

Figure 2

Mean (± SE) stress response of all 12 RGS U. richardsonii matrilineal-family groups with N= 41 females. Sample sizes for each matriline are displayed in parentheses and the overall mean stress response of all individuals (1.18) is illustrated as the dotted line. Pair-wise significance using the Tukey HSD test is denoted with an asterisk and the associated P-value.

Discussion

BL- and SI-cortisol and sample latency

The strong positive relationship between BL- and SI-cortisol within each subject allowed the use of a log-proportional linear stress response measure, as has been employed where individuals are challenged with an acute stressor such as physical handling or a novel environment (Boonstra et al. 2001; Atwell et al. 2012; Hare et al. 2014). This finding is consistent with the correlation between baseline and stress-induced measures of corticosterone among birds, as revealed in a meta-analysis performed by Bókony et al. (2009). Additionally, neither BSL nor SSL affected BL-cortisol, SI-cortisol, or the stress response. This confirms that our samples were drawn within the timeframe suggested by Hare et al. (2014) for RGS plasma glucocorticoids to be detectable in response to an acute stressor (ca. 7 min post-capture) and that HPA activity increased plasma cortisol with no apparent decline since the physiological challenge.

While age may have a pronounced effect on an individual’s stress response (e.g., Elliott et al. 2014), we found no effect of age in the present study. This may be attributable to the limited number of adults sampled (age-3, N = 4; age-4, N = 2), or to correlated stress responses within matrilineal families. The time of day at which the samples were taken affected BL- and SI-cortisol, but not stress response. This is consistent with previous reports that circadian rhythms influence circulating cortisol (Coe and Levine 1995; Boonstra et al. 2001; Nunes et al. 2006; Chauke et al. 2011). Diurnal mammals exhibit a steady decline in cortisol throughout the day and afternoon until midnight, at which point the trend begins to reverse (Tsang et al. 2014). Our data highlight this decline in daytime BL- and SI-cortisol, but this temporal trend is eliminated when both measures are considered together as the stress response.

A positive correlation between the first and second measures of the stress response was observed when subjected to a 0-constrained regression (r2 = 96.8%). Further, the differences between the first and second sampling for BL-, SI-cortisol, and stress response were not significantly different from 0. High repeatability of the stress response indicates that our measures for each individual are relatively consistent for the first and second samplings. This conclusion is supported by the fact that the number of days between first and second samples did not affect the magnitude of the difference in stress responses between repeated measures in 2014. Both baseline and stress-induced glucocorticoids are repeatable when using integrative and point-based measures, respectively (Bosson et al. 2009; Rensel and Schoech 2011; Cook et al. 2012; Grace and Anderson 2014; Small and Schoech 2015); however, glucocorticoid concentrations may be influenced by the environment (Smith et al. 2012; Baugh et al. 2014; Boulton et al. 2015) through environmental variation in diet (Boonstra et al. 2001; Dantzer et al. 2011), or predation events (Boonstra et al. 1998; Scheuerlein et al. 2001; Clinchy et al. 2011; Clinchy et al. 2013).

While not directly comparable to our fixed-intercept regression model or to our integrated stress response measure, similar studies in other species examining SI-plasma glucocorticoids and related measures have demonstrated high repeatability. In sciurids, comparable fecal glucocorticoid repeatability has been observed in yellow-bellied marmots Marmota flaviventris, Eurasian red squirrels Sciurus vulgaris (R = 0.52), and Columbian ground squirrels U.columbianus (R = 0.57), with a general trend toward wild-caught individuals exhibiting higher and less consistent BL-cortisol levels than their captive counterparts (Bosson et al. 2009; Dantzer et al. 2016; Smith et al. 2012). Repeatability studies in Sciuridae have generally employed integrative measures and emphasized BL-cortisol. In other non-mammalian species, repeatability is also high for stress response-related measures (multiple linear regression r2 = 0.50 in Lepomis macrochirus, Cook et al. 2012; R = 0.50–0.67 in Aphelocoma coerulescens, Rensel and Schoech 2011; Small and Schoech 2015). Often there is high repeatability for stress-induced plasma glucocorticoids (R = 0.32), but little or no repeatability for baseline-circulating levels (Rensel and Schoech 2011; Grace and Anderson 2014; Small and Schoech 2015). Other studies have found weak support for the repeatability of SI-cortisol (∼5% repeatable; e.g., Boulton et al. 2015), or no support for the repeatability of either measure (e.g., Baugh et al. 2014). Measured concentrations of both BL- and SI-cortisol were similar, however, within first and second samples of stress response obtained from our wild-caught females, and thus both appear to be somewhat consistent within individuals over our study. Additionally, there was no evidence of habituation to the acute stressor within individuals, and thus repeated samples were averaged in further analyses. The minimal latency between first and second sampling (mean ± SE: 4 ± 1 days) likely contributes to the repeatability of the concentrations measured in our study (Boulton et al. 2015).

Heritability

The heritability of physiological traits is likely governed by genetic and environmental factors, and interactions thereof, that can have a synergistic effect on offspring physiological development (Gillespie et al. 2009; Monclús et al. 2011). Our results suggest that the stress response is at least partially heritable in RGSs (mother–daughter h2 = 0.40 ± 0.24) and there is some evidence of additive genetic variance contributing to variation in the stress response. However, females in this species typically inherit a portion of their mother’s territories (Michener 1981) and thus likely experience similar environmental conditions to their mothers, which may lead to phenotypic variation being falsely attributed to genetic variation. Our estimates for full-sibling and half-sibling heritability were h2FS = 0.37 ± 0.71 and h2HS = 0.75 ± 1.41, respectively. Our mother–daughter stress response heritability is more similar to full-sibling heritability estimates, despite multiple paternity within litters (Hare et al. 2004). However, the full- and half-sibling heritability estimates presented here must be interpreted cautiously, owing to their extensive variance, which may reflect the behaviorally plastic nature of the stress response in accord with variation in local environmental conditions (Rensel and Schoech 2011; Jenkins et al. 2014; Homberger et al. 2015), or other factors considered above.

Many species exhibit similar narrow-sense heritabilities to those observed in our study. Artificial selection experiments on zebra finches Taeniopygia guttata reported heritabilities of h2 = 0.08 and h2 = 0.24 when selected for low and high stress-induced corticosterone levels (Evans et al. 2006), respectively, and similarly in Japanese quail Coturnix japonica, h2 = 0.14 and 0.30, respectively (Odeh et al. 2003). Outbred CD-1 mice selected for low- and high-stress reactivity exhibit similar narrow-sense heritability to that reported in our study (h2 = 0.40, Touma et al. 2008). For lab-raised rainbow trout Oncorhynchus mykiss and common carp Cyprinus carpio, SI-cortisol levels as a result of stress are strongly heritable (h2 = 0.41 and h2 = 0.60, respectively; Pottinger and Carrick 1999, Tanck et al. 2001), whereas European sea bass Dicentrarchus labrax appear to have negligible heritability of SI-cortisol when subjected to confinement (Volckaert et al. 2012). A meta-analysis of human twin studies (N = 399) by Bartels et al. (2003) reported a high narrow-sense heritability estimate of BL-cortisol (h2 = 0.62). In wild barn swallows, baseline and stress-induced corticosterone levels show some level of heritability (h2 = 0.152 and h2 = 0.343, respectively; Jenkins et al. 2014). Wild barn owls Tyto alba also exhibit similar physiological regulation of stress-induced corticosterone as their mothers (Almasi et al. 2010). The extent of heritability thus depends on the species, experimental setting (e.g., artificial selection vs. wild line, lab-reared vs. field experiment), and the measure in question (baseline vs. stress-induced glucocorticoids). Yet in most cases, both baseline and stress-induced glucocorticoids appear to be heritable, as observed for the stress response measure in the present study.

Groups of RGSs defined by direct descent or by area inhabited did not differ in terms of stress response from one another. However, the stress responses of distantly related individuals in the same larger matrilineal-family group differed from those of individuals in other matrilineal groups. Thus, while mothers, daughters, and sisters did not differ in their stress responses when viewed on a larger scale, matrilines composed of great-grandmothers, grandmothers, mothers, daughters, and sisters do differ from other matrilines in terms of how they respond to an acute stressor. No differences between immediate-family groups may be attributable to the reduced statistical power of tests examining only direct familial associations due to the limited sample size within each immediate-family group. For example, the average sample size per each matrilineal-family group was 3.42, whereas immediate-family groups averaged 2.32 individuals sampled per group, and sibling groups had 2.14 individuals per group. Therefore, it is possible that the latter 2 comparisons lacked statistical power relative to the comparison among matrilines. That said, spatial location groups were sampled nearly 2–4 times more often than all other comparisons, with an average of 8.20 individuals sampled per group. Therefore, if a spurious relationship were to be detected, it would be most likely to occur among groups of individuals by spatial location and not among matrilines. Additionally, our modest heritability coefficients calculated between mother and daughters do suggest an important genetic component associated with the stress response.

Different matrilineal groups may have experienced differing selection for the stress response, possibly to cope with unique past local environments, and these differences are likely inherited over a broad familial scale. Increased selection for the stress response in some matrilines may have left a signal in the larger matrilines, but not in the smaller immediate-family groups. Such a selective regime would result in squirrels showing similar responses when considering broader ancestry, but reveal no apparent differences when contrasting immediate-family groups or siblings, particularly in light of recent reductions in selection for the stress response or increased intra-familial variation induced by multiple paternity (Hare et al. 2004). Juvenile RGSs inherit part of their mother’s territory after 9–10 weeks, form kin clusters with their mothers, and aid in maintaining and defending this territory aggressively when non-kin or males intrude (Michener 1981, 1985). However, spatial area did not significantly predict stress response and thus genetic variation may better explain the differences in stress response among matrilineal-family groups, rather than environmental variation as observed in other studies (Boonstra et al. 2001; Dantzer et al. 2011; Baugh et al. 2014).

As previously discussed, the adaptive nature of the stress response is evident in terms of maternal influences shaping the stress response early in the offspring’s life (Sheriff et al. 2010a; Homberger et al. 2015), the high repeatability of point-based measures of glucocorticoids in other species, and increased plasticity in response to environmental perturbations relative to baseline glucocorticoids (Rensel and Schoech 2011; Jenkins et al. 2014; Homberger et al. 2015). Thus, selection may shape the stress response more so than BL-cortisol in RGSs. Studies of both BL-cortisol and the stress response have found that the individual fitness is better explained by the latter (Vitousek et al. 2014; Homberger et al. 2015).

Overall, our results demonstrate the variable nature of the stress response and highlight the apparent heritability of stress responses among related individuals. Thus, while distantly related squirrels may share a similar stress response, variation in that response occurs even among closely related individuals overlapping in their use of space, and thus that experience similar selective regimes. This suggests that shared genetic variation produces similar stress responses on a broader scale of matrilineal family, while immediate families do not differ significantly from one another, highlighting the inherent variability in proximate measures of HPA axis activity. Our findings have important methodological implications for studies measuring stress hormones, as well as applications in personality and behavioral research, evolutionary ecology, and population ecology. Future studies should examine the proximate causes of correlated stress responses among related individuals. Direct manipulation of circulating cortisol would also prove useful in teasing apart the contributions of genetic versus environmental factors shaping the mammalian stress response.