Testing hypotheses about individual variation in plasma corticosterone in free-living salamanders

Across vertebrate taxa, activation of the hypothalamic-pituitary–adrenal/interrenal (HPA/I) axis leads to release of circulating glucocorticoid (GC) hormones such as corticosterone (CORT). In most vertebrates, the HPA/I axis is activated by neural integration of internal and external stimuli that stimulates hypothalamic release of corticotropin-releasing hormone (CRH). CRH acts upon the anterior pituitary to release adrenocorticotropic hormone, which in turn acts on the adrenal/interrenal gland to release GCs. Baseline circulating levels of GCs are necessary for survival and have several homeostatic functions such as maintaining vascular tone and blood glucose (Norris and Carr, 2013). The HPA/I axis is also activated when encountering stressors and induces both behavioral and physiological changes geared towards helping the animal avoid, cope with or counter the original stressor (Sapolsky, 2002; Sapolsky et al., 2000).

Despite having many crucial roles, individual variation in activation of the HPA/I axis is frequently observed. HPA/I activity fluctuates with both daily and life-history events like development, reproduction, migration, etc. (Breuner et al., 1999; Dallman et al., 1993; Kalsbeek et al., 2012; Landys et al., 2006; Thurmond et al., 1986). A review of GC concentrations in free-living vertebrates found that many species of amphibians, reptiles and birds seasonally modulate GC release, with both baseline and stress-induced levels being highest during times of breeding (Romero, 2002). Although far less studied than GCs in terms of patterns of variation, levels of CRH peptide are also variable (Crespi and Denver, 2005; Denver, 1996, 1997; Matsuda et al., 2010; Stengel and Tache, 2015).

There are many hypotheses for why the HPA/I axis, and, in particular, GCs like CORT, is variable (Table 1). The hypotheses can be categorized as follows: the current-condition hypothesis (closely related to the fitness hypothesis), the facilitation hypothesis (also called the adaptation hypothesis) and the trade-off hypothesis (also called the energy re-allocation hypothesis) (reviewed in Patterson et al., 2014). Broadly, the current-condition hypothesis proposes that elevated baseline CORT levels are indicative of poor condition that may translate into reduced performance or reduced fitness of either an individual or a population (Bonier et al., 2009a; Husak and Moore, 2008). This is supported by evidence for negative correlations between baseline CORT and body condition/habitat quality in various birds and reptiles (Johnson, 2007; Müller et al., 2007; Oppliger et al., 1998; Waye and Mason, 2008; Williams et al., 2008). The facilitation hypothesis is based on the idea that GCs such as CORT have a role in glucose metabolism at all times, not just during challenges (Moore and Jessop, 2003). In this way, fluctuations in baseline CORT are necessary to accommodate the changing metabolic demands of various physiological processes such as reproduction and immune function. For example, there is accumulating evidence for positive relationships between GCs and sex steroid hormones, and GCs are elevated during breeding in many species (Bonier et al., 2009a; Moore and Jessop, 2003; Moore et al., 2000; Romero, 2002). Likewise, the HPA/I axis is intricately linked to the immune system and can both affect and be affected by immune activity (Martin, 2009; Turnbull and Rivier, 1999; Webster Marketon and Glaser, 2008). Thus, the facilitation hypothesis predicts that plasma CORT would be positively correlated with energetically costly functions during non-challenging (baseline) times. Conversely, the relationship between CORT and energetically costly functions that are not imperative to survival may shift to negative during times of stress, which is the basis for the trade-off hypothesis, also sometimes called the re-allocation hypothesis (Almasi et al., 2013; Martin et al., 2012; McEwen and Wingfield, 2003; Patterson et al., 2014; Wingfield and Sapolsky, 2003). The trade-off hypothesis predicts that, during times of stress, plasma CORT levels are negatively correlated with functions like reproduction or immune processes that are both energetically costly and not necessary for survival of the threat.

Here, we tested predictions generated by the hypotheses for HPA/I axis variation using data from free-living male and female Allegheny Mountain dusky salamanders (Desmognathus ochrophaeus) (Table 1). We included both sexes because the sexes often differ in patterns of HPA/I activity (Astheimer et al., 1994; Grassman and Hess, 1992a,b; Handa et al., 1994; Rhodes and Rubin, 1999; Tilbrook et al., 2000). These stream-side salamanders have strong seasonal cycles of reproductive activity, energy stores, body condition and plasma CORT (Fitzpatrick, 1973; Petranka, 2010; Ricciardella et al., 2010; Tilley and Tinkle, 1968). Specifically, we asked whether variation in physiological variables could explain patterns of baseline CORT, stress-induced CORT and hypothalamic CRH neurons. We selected physiological variables known to be associated with circulating GCs (Bonier et al., 2009a,b, 2011; Davis et al., 2008; Romero, 2002). We approximated reproductive investment with gonad mass, corrected for body size. We approximated energy stores with abdominal fat body mass, corrected for body size. Finally, we approximated body condition with carcass mass, corrected for body size. To approximate immune function, we examined differential white blood cell counts. The relative amounts of lymphocytes and neutrophils have been linked to circulating CORT in many field and laboratory studies (Davis et al., 2008; Goessling et al., 2015). In general, increases in relative blood neutrophil counts are interpreted as reflecting increased synthesis of these phagocytes in preparation for immune challenges. Likewise, decreases in relative lymphocyte counts are often interpreted as reflecting redistribution to the skin, sites of potential infection and injury (Dhabhar et al., 1996). However, it should be noted that decreases could also reflect reduced cell synthesis or cell death, especially because laboratory studies have linked GCs to lymphocyte apoptosis (Cidlowski et al., 1996).

All animal methods and procedures were approved by Duquesne University's Institutional Animal Care and Use Committee (1209-11). Collecting permits were obtained from the Pennsylvania Fish and Boat Commission. Allegheny Mountain dusky salamanders (Desmognathus ochrophaeus Cope 1859) are stream-side salamanders with a seasonal lifestyle, being active above ground from early spring to late autumn, and retreating underground during the winter. The mating season extends from autumn to spring with an interruption due to winter. Mating pauses in early summer, during which time females lay their eggs and brood their clutches for several months (Fitzpatrick, 1973; Petranka, 2010; Ricciardella et al., 2010; Tilley and Tinkle, 1968).

Adult male and female Allegheny Mountain dusky salamanders were collected along Elk Rock Run in Fayette County, PA, USA (39N°57′05.9″, 79W°34′43.4″). Collections took place across several months to acquire a multi-seasonal data set: spring (mating season) collection occurred from 15 to 30 May 2013; summer (non-mating season) collection occurred from 28 August to 4 September 2013; and autumn (mating season) collection occurred from 30 September to 2 October 2013. Animals were collected by overturning streamside rocks and logs and capturing animals by hand. Sexual maturity was verified by observing elongated premaxillary teeth in males and follicles (visible through the body wall) in females. Initial sample sizes were approximately 10 animals per sex per treatment (baseline versus 30 min after capture) for each sampling season; however, final sample sizes varied as a result of factors such as insufficient plasma volume for measuring CORT, lost tissues, poor blood smear slide quality, poor brain tissue sectioning, etc. We chose a sample size of 10 animals per treatment combination based on a Monte Carlo simulation power analysis (Bolker, 2008) that indicated that samples sizes of 8–12 are necessary to detect stress-induced increases in plasma CORT.

At each seasonal collection, we assessed baseline and stress-induced plasma CORT levels and the number of immunoreactive CRH neurons in both males and females. For baseline CORT, animals were euthanized via decapitation and trunk blood (approximately 3–10 μl) was collected in a heparinized capillary tube within 3 min of capture. For stress-induced CORT, animals were captured, individually placed in Ziploc bags, palpated for 2 min, and contained in the Ziploc for 30 min. This method of capture and handling has been shown previously to increase plasma CORT in Allegheny Mountain dusky salamanders (Woodley et al., 2014). Immediately afterwards, animals were euthanized via decapitation and trunk blood was collected. Animals were assigned to baseline or stress-induced treatments by alternating between treatments in a sequential manner. All blood samples and carcasses were kept on ice until they were processed in the lab, approximately 2–5 h later.

Trunk blood was collected in heparinized capillary tubes and centrifuged to collect plasma; plasma was frozen in heparinized microcentrifuge tubes until being analyzed via radioimmunoassay (RIA). Plasma CORT was assayed at the Endocrine Technology Services Core at the Oregon National Primate Research Center following standard methods. Briefly, a double ether (100% Honeywell Burdick and Jackson, no. 107-4) extraction was performed on up to 3 μl of plasma, after which concentrations of CORT were obtained via RIA using an anti-corticosterone antibody (ab77798 Abcam, 1:20,000 dilution, cross-reactivity with aldosterone is 0.06%). All samples were run in singlicate in a single assay that had an intra-assay coefficient of variation of 11.2%, a recovery of 86.4% and a sensitivity of 1.0 ng ml⁻¹. The assay has been validated for Allegheny Mountain dusky salamanders by demonstrating parallelism from 0.5 to 4 μl at 0.5 μl increments.

At the time of blood collection, a small drop of blood was used to make a blood smear on a clean glass slide. Smears were air-dried and stained using Wright–Giemsa stain (Polysciences, catalog number 24985) as follows: 30 s methanol (100%) dip, 2 min in Wright–Giemsa stain, rinse in phosphate buffer (pH 6.8, Polysciences, catalog number 24984), 4 min in Wright–Giemsa:buffer mixture, 2 min rinse in phosphate buffer. Slides were allowed to air dry completely before coverslipping with Permount. All slides were observed using oil immersion (1000×) with a compound light microscope. The investigator was blind to treatments. White blood cell count differentials were determined according to standard methods (Davis et al., 2008). Briefly, slides were incrementally scanned in a zigzag pattern, and individual cell types (lymphocytes, neutrophils, monocytes, basophils and eosinophils) were tallied until a total of 100 white blood cells had been counted. Slides with poor quality staining (grainy appearance or poorly delineated cells) were excluded from analyses prior to unblinding.

Carcasses were dissected to remove abdominal fat bodies. Both left and right abdominal fat bodies were massed together on a microgram balance. Carcasses were tagged and preserved with 10% neutral buffered formalin for 24 h. Following a 24 h water rinse, carcasses were stored in 70% ethanol. Gonads (testes or ovaries) were later excised and weighed from each preserved carcass. Also at this time, carcass mass (consisting primarily of muscle mass and gut, not including head, gonads and abdominal fat bodies) and length (forelimb to hindlimb) were recorded for calculating body condition (see statistics below).

Immediately following decapitation, the lower jaw was removed and each head (including the brain) was placed in a vial of freshly prepared 4% paraformaldehyde to preserve the tissue. Vials were agitated on a shaker overnight and then heads were rinsed with distilled, deionized water for 24 h. Brains were dissected from the heads and placed in 15% sucrose in phosphate-buffered saline (PBS) for at least 24 h, or until they had sunk to the bottom of the vial. They were transferred into 30% sucrose in PBS for another 24 h, followed by a final 24 h in a 1:1 solution of 30% sucrose and optimal cutting temperature (OCT) compound (Tissue Tek, catalog number 4583). Following cryoprotection, brains were embedded in OCT and snap-frozen using liquid nitrogen. Blocks were stored at −80°C until cryosectioning. Using a cryostat (International Equipment Company Minotome, OM 2488), alternating coronal sections (26 μm) were collected in two series by thaw-mounting onto polylysine-coated Superfrost Plus slides. Slides were allowed to air dry before being stored at −80°C until immunohistochemical analysis for CRH.

To visualize cells that were immunoreactive for CRH, standard immunohistochemistry was performed on one slide series using a primary rabbit anti-human/rat CRH antibody (code no. PBL rC68 from Dr Benedict Kolber, Duquesne University, originally from Dr Paul Sawchenko, Salk Institute). At room temperature, slides were washed in PBS containing 0.3% Triton X-100 (PBST), incubated in PBS with 3.0% H₂O₂, and washed in PBST again. Slides were incubated in a blocking solution of 5% non-fat dry milk in PBS for 1 h. Incubation in primary antibody was conducted overnight at 4°C with slides lying flat in a humidified chamber. Antibody was first reconstituted 1:50 in PBS to make a stock solution, and then diluted an additional 1:1000 in PBS containing 5% non-fat dry milk and 10% normal goat serum and applied directly to slides with a micropipette. Following incubation in primary antibody, slides were washed in PBST before incubation in secondary antibody, which consisted of a 1:200 dilution of biotinylated goat anti-rabbit in PBS containing 5% normal goat serum and 2% bovine serum albumin. Incubation occurred for 2 h at room temperature with slides lying flat in a humidified chamber. Slides were then washed in PBST before incubation in Vectastain Elite ABC (made according to kit instructions, Vector Laboratories, catalog number PK-6100) for 1 h. Following a wash with PBS, slides were incubated in tyramide signal amplification plus cyanine 3 (made according to kit instructions, TSA Cy3, PerkinElmer, NEL744001KT) for 3 min. Following a final wash in PBS, slides were allowed to air dry before coverslipping with Prolong Gold antifade mountant with DAPI (ThermoFisher Scientific, catalog number P36931). Slides were stored at 4°C until viewed and photographed (no more than 1 week post-immunohistochemistry; see below).

We confirmed that the primary antibody was specific for CRH and not for related neuropeptides in Allegheny Mountain dusky salamanders by doing a preabsorption control. Briefly, the primary CRH antibody was incubated overnight in a 50 μg ml⁻¹ solution of its immunogen, human/rat CRH (Sigma-Aldrich, CAS number 86784-80-7). The preincubated primary antibody was then used for immunohistochemistry on test slides containing sections from anatomical regions where CRH immunoreactivity is known to be present. Adjacent slides were treated with CRH antibody that had not been preabsorbed with CRH. Preincubating the primary CRH antibody with its antigen (human/rat CRH) resulted in a loss of fluorescence signal throughout the brain, observed in both fibers/puncta and cell bodies (Fig.1).

Stained slides were viewed and photographed using a TRITC filter on an epifluorescence microscope at 200× magnification. A thorough examination of brain tissue sections was conducted with reference to a brain atlas for plethodontid salamanders (Laberge et al., 2008). Although immunoreactive fibers were broadly distributed throughout the brain, most cells were observed in the ventral preoptic area, the magnocellular preoptic area (mPOA) and the locus coeruleus, with a small number of cells in the subpallial amygdala. While blind to treatment, neurons immunoreactive for CRH were counted. We limited our analyses to the mPOA CRH population because it is the primary population of hypophysiotropic CRH neurons, homologous to the hypophysiotropic paraventricular nucleus in mammals (Fasolo et al., 1984; Tonon et al., 1986). Slides with poor quality tissue and/or staining were excluded from analyses while blind to treatment.

To test the predictions generated by the hypotheses for variation in HPA/I activity (Table 1), we used information theory, specifically the Akaike information criterion adjusted for small samples sizes (IT-AICc) (Burnham and Anderson, 2002, 2014; Symonds and Moussalli, 2011). IT-AIC is a statistical method based on linear regression that allows multiple models to be evaluated and compared while avoiding over-fitting. Using R (https://www.r-project.org/), we modeled the response variables (baseline CORT, stress-induced CORT and the number of CRH neurons in the mPOA) with predictor variables that reflect different physiological functions. We analyzed baseline and stress-induced CORT separately because they were measured in separate animals, and there are separate hypotheses and predictions for baseline versus stress-induced CORT. For CRH neurons, we pooled data from baseline and stressed animals because there was no difference in the number of CRH neurons in baseline versus stress-induced samples (see Results). Note that only those subjects with complete data sets could be included in the modeling. Hence, sample sizes were 38 for baseline CORT, 43 for stress-induced CORT and 66 for mPOA CRH. To satisfy assumptions of the linear regression modeling approach, CORT and mPOA CRH values were square root transformed to achieve normality and homogeneous variance.

To assess the impact of innate immune function on the response variables, we used percentage lymphocyte (for baseline CORT, stress-induced CORT and CRH neurons) to reflect one aspect of immune function. We could not use both percentage lymphocyte and percentage neutrophil in our models because they co-vary, and collinearity should be avoided for IT-AIC. We included percentage lymphocyte in the models rather than percentage neutrophil or the neutrophil:lymphocyte ratio (NL ratio, a composite measure frequently assessed in relation to stressors) because percentage lymphocyte was not altered by the handling stressor whereas neutrophils decreased with the handling stressor (see Results).

To reflect energy stores, we included fat body mass as a predictor variable. To reflect reproductive investment, we included gonad mass. To reflect body condition, we included carcass mass (which excluded abdominal fat bodies and gonads) as a predictor. Because fat body mass, gonad mass and carcass mass are related to overall body size, we corrected these measures for body length using ANCOVA with body length as the covariate. To do so, we first transformed the mass variables to achieve normality and homogeneity of variance using square root (fat body mass, gonad mass) or log (carcass mass) transformations. Next, we checked for linearity using scatterplots. Finally, we confirmed that the relationship between the mass variable and body length was the same for males and females by running an ANCOVA with sex as a factor and body length as a covariate. Because slopes were homogeneous (the sex × length interaction was non-significant) for fat body mass and carcass mass, we saved the unstandardized residuals as the size-corrected variables. For gonad mass, slopes describing the relationship between gonad mass and body length differed between males and females. Therefore, we did separate ANCOVA for gonad mass for males and females to obtain residuals to represent the size-adjusted gonad masses.

Although males and females clearly differ in CORT levels, we did not include sex alone as a predictor variable, because inclusion might mask relationships between the physiological and endocrine variables of interest. Instead, we assessed the impact of sex by including interaction terms of sex and the other predictor variables to assess whether the effects of the predictor variables might differ according to sex.

As shown in the Results, use of IT-AIC did not identify clear ‘best’ models. Instead, many similar models were supported by the analysis. For example, in most cases the weight of the best model was much less than 0.9, meaning that the probability that it was the best model was less than 90%. Also, several models had similar weights with ΔAICc values less than 2. Models with ΔAICc values less than 2 are considered to be as good as the best model, and models with ΔAICc values up to 6 may also be valid (Symonds and Moussalli, 2011). Therefore, to better understand which predictor variables are most likely to be part of the best model, we performed full-model averaging, which allows inferences to be made from the full set of candidate models. Full-model averaging sums the Akaike weights across all the models that contain a particular predictor variable to give a measure of the probability that the predictor variable is in the best model. It also provides an estimate of the mean (β) and error [se(β)] for each predictor variable; 95% confidence intervals were then determined with the equation: β±1.96[se(β)], to assess the magnitude of the estimate and whether it was distinguishable from zero.

The IT-AICc approach described above specifically tested the hypotheses in Table 1, but we also graphed and analyzed how the variables differed according to season, sex and handling to provide context for the IT-AIC results. Hence, we performed 3-way ANOVA using IBM SPSS Statistics 22 to determine the effects of season, sex and handling after transforming data to meet parametric assumptions (normality and homogeneous variance). Analyses of fat body mass, gonad mass and carcass mass included body length (forelimb to hindlimb length) as a covariate to correct for body size.

Using IT-AICc, candidate models were generated from data from 38 animals. The model containing percentage lymphocyte and the sex×percentage lymphocyte interaction was the best approximating model for variation in baseline CORT (Table 2, weight=0.185). The best-fit model demonstrated goodness-of-fit (R²=0.788, adjusted R²=0.776). However, because the weight of the top model indicated only an 18.5% probability that it was the best among the set of candidate models, full-model averaging was calculated for all generated models. This identified percentage lymphocyte and the sex×percentage lymphocyte interaction as the strongest contributors to variation in baseline CORT (Table 3, lymphocyte: w=0.999; sex×lymphocyte interaction: w=0.999). Additionally, 95% confidence intervals for both predictors did not contain zero (percentage lymphocyte: −0.098, −0.044; sex×percentage lymphocyte interaction: 0.03, 0.04). Although it is difficult to visualize multivariate relationships in 2-dimensional space, we present a scatterplot to help understand the sex×percentage lymphocyte interaction in Fig. 2. The scatterplot suggested that the negative relationship between baseline CORT and percentage lymphocyte was particularly evident in females (r=−0.749, P=0.001) but not males (r=0.195, P=0.398).

Candidate models were generated using data from 43 animals. Data from one animal (male, stressed, spring sample) was a highly influential point in model selection (Cook's distance: D≈2.85) because of a very low plasma CORT and small carcass mass. By convention, D>1 is considered to be influential, so we removed data for this animal from further analysis (Cook, 1977). For stress-induced CORT, a model containing percentage lymphocyte and a sex×percentage lymphocyte interaction was identified as the best approximating model for variation (Table 2, weight=0.115). Regression analysis of the best-fit model (percentage lymphocyte, sex×percentage lymphocyte interaction) demonstrated goodness-of-fit (R²=0.673, adjusted R²=0.657). However, the overall probability that this model was the best was 11.5%, and there were three models with ΔAICc<2. Because it was unclear which model was the best, we conducted full-model averaging to determine which predictors were likely to be a part of the best model. This identified the sex×percentage lymphocyte interaction and percentage lymphocyte as the strongest contributors to variation in stress-induced plasma CORT (Table 3, sex×percentage lymphocyte interaction: w=0.996; percentage lymphocyte: w=0.961). Additionally, 95% confidence intervals for both variables did not contain zero (sex×percentage lymphocyte: 0.02, 0.04; percentage lymphocyte: −0.061, −0.011). A scatterplot of percentage lymphocyte and stress-induced CORT did not provide insight for the sex×percentage lymphocyte interaction because bivariate correlations were weak and non-significant in both sexes (males: r=0.299, P=0.2; females: r=0.034, P=0.875).

For CRH neurons in the mPOA, models were generated using data from 66 individuals. As with stress-induced CORT, the influential male observation was removed (see above). A model containing percentage lymphocyte, a sex×percentage lymphocyte interaction, relative fat and relative gonad, and a sex×relative gonad interaction was identified as the best-approximating model for the variation (Table 2, weight=0.1357). Linear regression analysis of the best-fit model (percentage lymphocyte, sex×percentage lymphocyte interaction, sex×relative gonad interaction, relative gonad, relative fat) demonstrated goodness-of-fit (R²=0.301, adjusted R²=0.243). However, the overall probability that this model was the best was 13.6%, and there were five models with ΔAICc<2. Because it was unclear which model was the best, full-model averaging was conducted, and this identified a sex×percentage lymphocyte interaction and percentage lymphocyte as contributors to the variation (Table 3, sex×percentage lymphocyte interaction: w=0.969; percentage lymphocyte: w=0.853). The sign of the coefficient was negative, indicating a negative relationship between percentage lymphocyte and mPOA CRH neurons. Additionally, 95% confidence intervals for both predictors did not contain zero (sex×percentage lymphocyte interaction: 0.005, 0.029; percentage lymphocyte: −0.073, −0.003). A scatterplot of percentage lymphocyte and CRH neurons in the mPOA did not provide insight for the sex×percentage lymphocyte interaction because bivariate correlations were weak and non-significant in both sexes (males: r=−0.093, P=0.578; females: r=−0.052, P=0.794).

Here, we present results on how the variables differed according to season, sex and handling. Although these results do not indicate which hypotheses (Table 1) are supported by the data, they provide context for the IT-AIC results. Plasma CORT levels were higher in males than in females, and handling significantly elevated plasma CORT (Fig. 3A, sex: F_1,105=247.04, P<0.001; treatment: F_1,105=109.69, P<0.001). However, there was no overall effect of season on plasma CORT (season: F_2,105=1.47, P=0.234). The number of CRH neurons in the mPOA did not vary between treatments or across season, but males had more CRH neurons than females (Fig. 3B, treatment: F_1,82=1.23, P=0.271; F_2,82=0.56, P=0.574; sex: F_1,82=14.27, P<0.001). Plasma CORT and the number of CRH neurons in the mPOA were positively correlated (r=0.411, P<0.001).

The percentage lymphocyte of the white blood cell count differential varied according to season but not sex or treatment (Fig. 4C, season: F_2,76=6.1, P=0.003; sex: F_1,76=0.30, P=0.59; treatment: F_1,76=0.32, P=0.86). The percentage neutrophil was similar between males and females and across season (Fig. 3D, sex: F_1,76=0.85, P=0.36; season: F_2,69=2.1, P=0.132), but decreased with treatment (treatment: F_1,76=9.7, P=0.003). Similar to the percentage neutrophil, the NL ratio was similar between males and females (Fig. 4A, sex: F_1,76=0.86, P=0.36) and decreased with treatment (treatment: F_1,76=7.5, P=0.008). There was a trend for NL to vary across season, being lowest in the autumn (season: F_2,76=2.96, P=0.058).

Body length, fat body mass, gonad mass and carcass mass did not differ between baseline and stress-induced treatments, so baseline and stress data were pooled (Fig. 5A, treatment: F_1,111=1.741, P=0.190; Fig. 5B, treatment: F_1,111=0.013, P=0.908; Fig. 5C, treatment: F_1,111=0.848, P=0.359; Fig. 5D, treatment: F_1,111=0.002, P=0.965). Body length was greater in the summer and autumn than in the spring, and males were longer than females (season: F_2,111=3.810, P=0.025; sex: F_1,111=17.375, P<0.001). Overall, fat body mass (corrected for body length by including body length as a covariate) increased from season to season, and females had more fat than males (season: F_2,111=35.95, P<0.001; sex: F_1,111=74.545, P<0.001). Ovaries were heavier than testes (sex: F_1,111=99.003, P<0.001) and ovarian mass was lowest in the summer while testes mass peaked in the summer (season: F_2,111=3.56, P=0.032; sex×season: F_2,111=18.772, P<0.001). Finally, both males and females increased in mass (corrected for body length by including body length as a covariate) over time, and males were heavier than females throughout the year (season: F_2,111=8.87, P<0.001; sex: F_1,111=21.14, P<0.001).

With this study, we amassed a large and variable data set on physiological and endocrine variables from male and female free-living Allegheny Mountain dusky salamanders to distinguish among alternative hypotheses explaining individual variation in patterns of HPA/I activity. Specifically, we used IT-AICc statistical methods to determine which physiological variables were most closely associated with individual variation in baseline CORT, stress-induced CORT and mPOA CRH. The best-fit models explained a substantial amount of the individual variation in baseline CORT (78%) and stress-induced CORT (67%), albeit lower for mPOA CRH (30%). All of the best-fit models, regardless of the measure of HPA/I activity, included percentage lymphocyte and a sex×percentage lymphocyte interaction. The probabilities that percentage lymphocyte and the sex×percentage lymphocyte interaction were part of the best-fit model were greater than 96% in all cases, except for percentage lymphocyte in the model of mPOA CRH neurons (85% probability). Furthermore, the 95% confidence intervals of the coefficients generated by model averaging did not overlap zero, further confirming that percentage lymphocyte and the sex×percentage lymphocyte interaction were associated with patterns of HPA/I activity. In contrast, measures of body condition, energy stores and reproductive investment did not contribute to the variation in plasma CORT or CRH neurons. Below, we discuss our results in more detail.

In the field of conservation physiology, it has been proposed that measurement of plasma CORT gives insight into the health and vigor of an individual or population (Tarlow and Blumstein, 2007; Wikelski and Cooke, 2006), with baseline CORT (and its releasing hormone, CRH) being negatively associated with body condition, and stress-induced CORT being positively associated with body condition. In Allegheny Mountain dusky salamanders, we found no support for the current-condition hypothesis. Neither body condition nor abdominal fat body mass (a measure of energy stores, presumably reflective of condition) explained individual variation in plasma CORT, either baseline or stress induced. These results are consistent with a study in Allegheny Mountain dusky salamanders that also found no relationship between body condition and CORT (Woodley et al., 2014). Thus, in Allegheny Mountain dusky salamanders, measurement of plasma CORT does not provide insight into the health of an individual or population.

Our lack of support for the current-condition hypothesis is consistent with the literature, which has shown mixed support for the CORT-fitness hypothesis, which is conceptually similar to the current-condition hypothesis. In several studies, the correlation between CORT and fitness was highly context dependent (Bonier et al., 2007, 2009b; Dantzer et al., 2014; Lancaster et al., 2008; Magee et al., 2006). Other studies have found unpredicted positive correlations between CORT and indicators of fitness (Cyr and Romero, 2007; Silverin, 1998). Interestingly, a recent meta-analysis found that integrated levels of CORT (such as can be measured in fecal samples or water-borne hormone analyses) were better associated with measures of chronic stress than were more instantaneous, snapshot, measures of CORT like plasma CORT (Dantzer et al., 2014). Thus, future studies should use more integrated measures of GCs in relation to health and fitness.

We found no relationship between individual variation in HPA/I activity and relative gonadal mass, suggesting that baseline CORT (and mPOA CRH) does not facilitate investment in reproductive effort, and that patterns of stress-induced CORT do not reflect trade-offs related to reproduction. In some species, the CORT response to acute stressors like handling is suppressed during the breeding season (reviewed in Wingfield and Sapolsky, 2003). This observation suggests that high stress-induced CORT is incompatible with reproductive functions because the elevated CORT would allocate energy away from reproduction, thus jeopardizing the reproductive effort. In our study of free-living salamanders, gonadal mass was not related to either baseline or stress-induced CORT, despite gonads representing a substantial portion of the total body mass, especially in females. The disconnect between plasma CORT and reproductive effort extends to reproductive behaviors; dusky salamander mating behavior was unaffected by endogenous or exogenous elevations in plasma CORT in laboratory studies (Bliley and Woodley, 2012; Woodley and Lacy, 2010).

The disconnect between CORT and reproductive effort contrasts with studies in other species. In white-crowned sparrows, a species where females do the majority of nest building, incubating and feeding of offspring, reproductive success was positively associated with baseline CORT and was negatively associated with stress-induced CORT (Patterson et al., 2014). Amphibians actively engaged in reproductive behavior often have elevated plasma CORT (Harvey et al., 1997; Leary et al., 2004; Mendonça et al., 1985; Orchinik et al., 1988; Reedy et al., 2014). One explanation for the species differences might relate to the length of the breeding season. In white-crowned sparrows, the breeding season is limited to a few weeks. In the amphibian examples, mating is often limited temporally, with an explosive mating season. In dusky salamanders, the breeding season lasts most of the year, including the autumn–spring mating season and the summer brooding period when females oviposit and guard egg clutches for several weeks. Given that physiological levels of CORT elevate metabolic rate (Wack et al., 2012), dusky salamanders may avoid elevating plasma CORT for such a prolonged period as part of their low-energy lifestyle (Feder, 1983).

Our results indicated that percentage lymphocyte was negatively associated with variation in baseline plasma CORT levels. The best-fit models for baseline CORT included percentage lymphocyte and a sex×percentage lymphocyte interaction and explained 78% of the variation, indicating a good fit. However, interpretation of the decline in blood lymphocytes is unclear (for a review of the interpretation of white blood cell count differentials, see Davis et al., 2008). On the one hand, the lower percentage lymphocyte could result from a redistribution from the blood to the periphery for monitoring potential sites of infection (Viswanathan and Dhabhar, 2005). If so, this result would be consistent with the facilitation hypothesis, which posits that elevated baseline CORT primes the immune systems for challenges. Interestingly, there was a clear sex difference, with females but not males showing a negative relationship between percentage lymphocyte and baseline CORT. The strong negative correlation in females may reflect up-regulation of immunity to protect females from infection during mating, which consists of a prolonged period of male–female physical contact during which the male scratches the female's dorsum with hypertrophied premaxillary teeth (Houck and Arnold, 2003). Indeed, the percentage lymphocyte in the blood is lowest during the spring, which is the peak of the mating season. Infection is associated with lymphopenia in other vertebrates, including salamanders (Davis et al., 2008; Hopkins et al., 2016), suggesting that decreased blood lymphocytes are important for combatting immune challenges.

On the other hand, the decrease in the percentage lymphocyte in the blood could result from a reduced synthesis of lymphocytes or death of lymphocytes (Cidlowski et al., 1996). If so, the negative association with lymphocytes is consistent with the trade-off hypothesis, whereby the immune system is suppressed to allocate energy to other important processes related to survival. Related to this, treatment of Allegheny Mountain dusky salamanders with physiological levels of CORT suppressed wound healing (Thomas and Woodley, 2015), perhaps reflecting a trade-off between healing and survival. To distinguish among these hypotheses, more measures of immune function and determining health outcomes of animals with particular immune profiles are necessary to understand whether the immune system is primed or suppressed in association with elevated plasma CORT. For example, after immobilization restraint, cururu toads had decreased bacterial killing ability along with increased plasma CORT and decreased percentage lymphocytes in the blood, suggesting suppression of immune function (de Assis et al., 2015). Study of additional immune factors would also indicate whether the relationships between HPA/I activity and white blood cell count differentials extend to other measures of innate or cell-mediated immunity.

In addition to measuring baseline CORT, we also counted numbers of hypophysiotropic CRH neurons because in most vertebrates, CRH is the primary releasing hormone for activation of the HPA/I axis. Hence, it was not surprising that mPOA CRH numbers were positively correlated with baseline plasma CORT. Furthermore, we found that percentage lymphocyte was a strong contributor to variation in the number of mPOA CRH, similar to baseline CORT. Similarly, the value of the estimating coefficient was negative, demonstrating a negative relationship between mPOA CRH neurons and percentage lymphocyte. Thus, the CRH results provided additional evidence that HPA/I activity is strongly related to immune factors.

White blood cells also contributed to variation in stress-induced plasma CORT. As with baseline CORT, it is difficult to interpret this result. Our findings could support the facilitation hypothesis or the trade-off hypothesis, depending on the functional consequences of a drop in percentage lymphocyte. Again, more measures of immune function and determining health outcomes of animals with particular immune profiles would aid interpretation.

The relationship between percentage lymphocyte and the different measures of HPA/I activity (baseline CORT, stress-induced CORT and CRH) differed between males and females. The sex difference was most evident for baseline CORT, where baseline CORT was negatively associated with percentage lymphocyte in females but not males. However, the sex differences were more difficult to assess for stress-induced CORT and CRH. Scatterplots and simple correlations did not show strong sex differences, suggesting that the sex differences were relatively subtle, perhaps reflecting slight differences in slopes. Nonetheless, the consistent sex×percentage lymphocyte interaction results indicate that additional studies of sex differences in immune function are warranted.

It is important to note that the associations revealed by the IT-AICc here do not indicate causality. In fact, interactions between HPA/I activity and immune function are bi-directional (Silverman and Sternberg, 2012). Baseline circulating GCs stimulate early increases in immune responses associated with a challenge, but they inhibit immune function long-term, perhaps to prevent the development of morbidities associated with autoimmune diseases (Besedovsky and Sorkin, 1977; Hardy et al., 2012; Munck and Náray-Fejes-Tóth, 1994; Sapolsky, 2002). Also, white blood cells and associated cytokines can influence GC variability via CRH-mediated pathways (Bethin et al., 2000; Chrousos, 1995; Dunn, 2000).

It was not our goal to determine whether endocrine and physiological variables changed on a seasonal basis, but to model how different physiological variables were associated with patterns of HPA/I activity. However, we graphed and analyzed the data according to season to illustrate the broad range of variation that we captured for our physiological variables. Most variables differed seasonally, such as body and fat mass, body length, gonad mass, and baseline percentage lymphocyte and percentage neutrophil. Although plasma CORT and the number of mPOA CRH neurons did not vary statistically across season, mean baseline CORT was slightly higher in the spring than in summer and autumn, as was found in a previous study (Ricciardella et al., 2010). Also, males and females differed for many of the variables, justifying our rationale to include interactions with sex in the IT-AICc modeling.

The simultaneous assessment of multiple physiological and endocrine variables in a free-living amphibian, including novel data for CRH neurons and white blood cells, represents an important contribution to the field of environmental endocrinology. With these data, we tested predictions about the meaning of patterns of HPA/I activity. Although it is often hypothesized that much of the variation in CORT reflects animal or population health, we found no support for the current-condition hypothesis. We also found a surprising disconnect between plasma CORT and reproduction. We identified the percentage lymphocyte as the strongest predictor of HPA/I axis variation compared with other physiological variables. We also found evidence for sex differences in the relationship between percentage lymphocyte and measures of HPA/I activity. The next step is to examine more measures of immune function to determine whether the relationships we found are reflective of other aspects of immune function or are limited to white blood cell count differentials.

We thank Woodley lab members for assistance with field work.