Evidence that stress-induced changes in surface temperature serve a thermoregulatory function

Changes in body temperature following perception of a stressor were first reported almost 2000 years ago. In the 2nd century CE, for example, Galen described the presence of an ‘ephemeral fever’ in humans that was thought to be evoked by an abundance of humour (Yeo, 2005); a concept later expanded upon by Ibn Sina (11th century CE) who described ‘fevers’ driven by emotions, including grief, anger and dread (reviewed in Parviz et al., 2013). To date, changes in body temperature following exposure to stressors have been reported across numerous non-hominid species (i.e. lizards: Cabanac and Gosselin, 1993; turtles: Cabanac and Bernieri, 2000; birds: Greenacre and Lusby, 2004, and fish: Rey et al., 2015) and the area has attracted significant research attention. Despite such long-standing recognition in philosophical and scientific literature, and near vertebrate-wide conservation (but see Cabanac and Laberge, 1998; Jones et al., 2019), the functional significance of stress-induced changes in body temperature remain elusive. Indeed, while many studies have done well to uncover proximate mechanisms underpinning this phenomenon (i.e. mediation by peripheral vasoconstriction, pyrogenic cytokines, prostaglandins or glucocorticoids; Yokoi, 1966; reviewed in Oka et al., 2001; and Jerem et al., 2018, respectively), remarkably few have empirically interrogated a functional role of changing one's body temperature, per se, in response to stress exposure (Cabanac and Gosselin, 1993).

Thermal responses to stress exposure can be broadly categorized according to their position of occurrence within an organism (i.e. at core tissues or surface tissues). Of these, core temperature responses have, perhaps, received the greatest theoretical attention with respect to functional significance (reviewed in Oka et al., 2001; Oka, 2018). For example, some immunological studies have posited that changes in core body temperature following stress exposure represent a true ‘fever’ (Singer et al., 1986; Sanches et al., 2002) and endow individuals with a defensive advantage if they experience injury and pathogen exposure during stress exposure (discussed in Oka et al., 2001). Others, however, have contested this hypothesis by failing to substantiate immune mediation of core temperature responses to stress (Long et al., 1990b; Soszynski et al., 1998; Hiramoto et al., 2009; Vinkers et al., 2009). At the level of surface tissues, however, a functional role of stress-induced changes in temperature, per se, appears yet to be raised (though briefly discussed in Herborn et al., 2018). Indeed, dominant theory explaining stress-induced changes in surface temperature posits that this phenomenon is merely a consequence of haemetic redistribution driven by peripheral vasoconstriction (Jerem et al., 2015; Jerem et al., 2018; Nord and Folkow, 2019) and holds no direct functional role; rather, it is haemetic redistribution, but not thermal modulation, that carries functional significance by attenuating blood loss in the event of injury (as long shown following haemorrhage by McGuigan and Atkinson, 1921; Freeman, 1932; Darlington et al., 1986).

Despite an absence of functionally guided research, studies describing stress-induced changes in surface temperature do allude to a functional value of this phenomenon. Nord and Folkow (2019), for example, reported that while the skin temperature of Svalbard rock ptarmigans (Lagopus muta hyperborea) typically falls after handling, the magnitude of the skin temperature response varies according to ambient temperature (T_a); specifically, ptarmigans handled at low temperature (−20°C) display a larger change in skin temperature than those handled at warmer temperatures (0°C). An effect of T_a on stress-induced changes in skin temperature suggests that the function of this phenomenon extends beyond haemetic redistribution and may provide thermoregulatory advantages, where heat conservation at low T_a (sub-thermoneutral) is enhanced during perception of environmental challenges. Supporting a thermoregulatory function to stress-induced changes in surface temperature, Herborn et al. (2018) described an increase in skin temperature of domestic chickens (Gallus gallus) after exposure to chronic stress treatments at constant, thermoneutral temperatures, with respect to controls. As exposure to stressors is thought to elevate metabolically generated heat (i.e. as a consequence of tachycardia, tachypnea and avoidance behaviour; Cabanac and Aizawa, 2000; Cabanac and Guillemette, 2001; Greenacre and Lusby, 2004; Long et al., 1990a), elevation of skin temperature under chronically challenging environments may facilitate dissipation of excess metabolically generated heat, thereby reducing thermal load associated with activation of a stress response. Neither Nord and Folkow (2019), nor Herborn et al. (2018), however, tested patterns of heat conservation or dissipation following stress exposure, rendering such thermoregulatory consequences of surface temperature responses speculative.

In this study, we propose that stress-induced changes in surface temperature, per se, serve a functional role that may be understood when contextualized according to energetic and thermal load. Specifically, we argue that changes in surface temperature following stress exposure reduce energetic costs that are incurred during activation of a stress response, by promoting heat conservation at low temperatures (conservatively, below thermoneutrality), and heat dissipation at high temperature (conservatively, above thermoneutrality). This hypothesis (henceforth defined as the ‘thermoprotective hypothesis’) contrasts the dominant hypothesis stating that surface temperature responses are a functionally neutral corollary of haemetic redistribution (discussed above, and in Jerem et al., 2015; henceforth defined as the ‘haemoprotective hypothesis’). Importantly, the thermoprotective hypothesis does not preclude a selective advantage to blood-loss avoidance on stress-induced peripheral vascular motion; rather, it dictates that such stress-induced vascular motion is further shaped by thermoregulatory requirements experienced during stress exposure.

To test whether stress-induced changes in surface temperature themselves serve a functional value, we generated predications for both the thermoprotective and the haemoprotective hypotheses, with reference to a temperate endotherm, the black-capped chickadee (Poecile atricapilus Linnaeus 1766). According to the thermoprotective hypothesis, we predicted that surface temperature and dry heat-loss (here, heat lost by radiation and convection, but not evaporative cooling or conduction) of chickadees would fall under stress exposure when T_a is low (below thermoneutrality), and rise when T_a is high (above thermoneutrality), thereby reducing energetic and thermal load during activation of a stress response. Alternatively, under the haemoprotective hypothesis, we predicted that surface temperature and dry heat-loss of chickadees would also fall under stress exposure, however, only until T_a and surface temperature are approximately equal (i.e. when acquisition of heat at surface tissues matches loss of heat from local ischemia), whereafter, surface temperature would remain equivalent to T_a. Because female parids are thought to differ in both metabolic costs of thermoregulation (Nilsson et al., 2011) and heat dissipation capacity compared with male conspecifics (i.e. as a consequence of gynolateral brood-patch development and sexual size dimorphism; Foote et al., 2010; Alonso et al., 2016; Cooper and Voss, 2013; discussed in Nilsson and Nord, 2018), we further predicted that thermal responses to stress exposure would differ between females and males under the thermoprotective hypothesis, but not the haemoprotective hypothesis. Specifically, if stress-induced changes in surface temperature follow the thermoprotective hypothesis, and therefore occur to offset energetic costs incurred during activation of the stress response, we predicted that female chickadees would exhibit a larger fall in surface temperature and dry heat loss following stress exposure than male chickadees when T_a is low (i.e. a lower intercept of the relationship between T_a and surface temperature following stress exposure). Additionally, when T_a is high, we predicted that under the thermoprotective hypothesis, male chickadees would display a larger increase in surface temperature and dry heat loss following stress exposure than females (that is, a steeper slope of the relationship between T_a and surface temperature following stress exposure), owing to their lower heat dissipation capacity.

To our knowledge, this is the first study to empirically test a direct functional role of stress-induced changes in surface temperature. It is also the first study to investigate stress-induced changes in surface temperature across an ecologically relevant T_a range that extends both below and above the thermoneutral zone of the study organism. In light of a changing global climate, elucidating a thermoregulatory function to the vertebrate stress response may provide insight into the adaptive capacity of organisms that face both physiological and thermal stress from anthropogenic, environmental change.

All methods used for animal capture, handling and experimentation were approved by the Trent University Animal Care Committee (AUP no. 24614) and Environment and Climate Change Canada (permit no. 10756E).

During the months of March and April in 2018, we captured 20 free-living black-capped chickadees within a 100 km² area of south-central Ontario (Canada) for captive experimentation. Because conspecific songbirds from urban and rural populations have been shown to exhibit differences in the magnitude of the stress response (e.g. Abolins-Abols et al., 2016), 10 birds (n_females=5; n_males=5) were captured from across known urban populations (Cambridge, Ontario, 43.3789°N, 80.3525°W; Guelph, Ontario, 43.3300°N, 80.1500°W; Brantford, Ontario, 43.1345°N, 80.3439°W) and 10 birds (n_females=5; n_males=5) were captured across three known rural populations (Erin, Ontario, 43.7617°N, 80.1529°W; Corwhin, Ontario, 43.5090°N, 80.0899°W; Ruthven Park National Historic Site, Ontario, 42.9797°N, 79.8745°W). All chickadees were trapped using remotely operated potter-traps (90×70×70 cm; l×w×h), baited with sunflower seeds and suet on the day of capture. To further attract individuals to trapping locations, black-capped chickadee breeding and mobbing songs were broadcasted alternately from a remote call-box (FoxPro™ Patriot; Lewisville, PA, USA) at approximately 80 decibels, which was stopped when at least one individual had approached within a 4 m radius of a baited trap.

Once captured, individuals were immediately blood sampled (50 µl) by brachial venipuncture and capillary tube collection (<5 min post-capture), then assigned a unique combination of one government issued, stainless steel leg band and two coloured darvic leg bands for identification. Individuals were then weighed (nearest 0.1 g using a digital platform scale), measured (tarsus and flattened wing-chord to the nearest 0.1 mm, using analogue callipers) and secured in covered carrier cages (30×30×15 cm; l×w×h) for transport to our long-term holding facilities at the Ruthven Park National Historic Site, Cayuga, Ontario, Canada (maximum of 90 km; 2 h travel by vehicle). Finally, erythrocytes were isolated from whole blood samples on site by centrifugation (12,000 rpm), then lysed in 500 µl Queen's lysis buffer for long-term preservation of genetic material (Seutin et al., 1991). Plasma was isolated and stored for another study (not described here). All lysed blood samples were held on ice until transfer to permanent holding at 4°C.

Once at our long-term holding facility, all chickadees were haphazardly assigned to one of four identical and visually isolated outdoor flight enclosures (n=5 per flight enclosure; 183×122×244 cm; l×w×h). Each flight enclosure was supplied with an insulated roosting box (60×20×20 cm; l×w×h) mounted at 1.2 m in height, one roosting tree (white cedar, Thuja occidentalis; 1.0 m), and two perching branches (80 cm in length), mounted at approximately 1.5 and 1.8 m in height. Chickadees were provided water, and a mixture of meal worms (Tenebrio molitor), house crickets (Acheta domesticus), shelled peanuts, apple pieces, boiled egg, sunflower seeds, safflower seeds and Mazuri (St Louis, MO, USA) Small Bird Maintenance diet ad libitum. Both food and water were distributed twice to three times daily, with food being exclusively dispensed on a 400 cm² platform, raised to 1.2 m in height. To minimize disturbance during feeding, all food and water were distributed through opaque, hinged doors (15 cm×15 cm), such that chickadees were blind to experimenter presence. All individuals were acclimated for a minimum of two weeks prior to experimentation.

From April until late June of 2018, we tested a functional role of stress-induced changes in surface temperature by using a repeated sampling design to maximize statistical power. Here, each individual was exposed to one control (untouched; n=30 days) and one stress exposure treatment (n=30 days), separated by a rest period of 2 days (total experimental duration=62 days). All chickadees were maintained in outdoor flight enclosures for the duration of the study (see above) and were therefore exposed to seasonal changes in temperature and day length that could alter patterns of stress responsiveness (Wingfield et al., 1992; Astheimer et al., 1995). To account for these seasonal changes, individuals within two flight enclosures (n=10; Group A) were exposed to a control treatment followed by a stress exposure treatment, while the remaining individuals (n=10 within two flight pens; Group B) were exposed to a reversed treatment order, such that stress exposure treatments in Group B co-occurred with control treatments in Group A, and control treatments in Group B co-occurred with stress exposure treatments in Group A.

Stress exposure treatments followed a protocol of rotational stressors similar to Rich and Romero (2005) and Cyr and Romero (2007); however, no auditory stressors were used to ensure that the application of a stressor to individuals within one flight enclosure did not elicit a stress response in individuals held within remaining and nearby flight enclosures. Because our flight enclosures were not auditorily separated, however, we could not control for the exposure of control individuals to alarm calls elicited by stress-exposed individuals; consequently, our analyses comparing surface temperature profiles between control and stress-exposed individual are expected to be conservative.

Experimental individuals were administered five, randomly selected passive stressors each day, with each stressor persisting for 20 min, and being separated from subsequent stressors by 1 h. Daily randomization of stressors was used to circumvent habituation to each individual stressor throughout the course of experimentation. Stressors included; capture and restraint, presence of a mock predator (a taxidermied adult Cooper's hawk Accipiter cooperii), presence of a novel object placed in the centre of holding enclosures (garden gnome), presence of a human within holding enclosures, absence of light (by wrapping enclosures with opaque fabric) and presence of a mounted conspecific placed in the centre of a feeding platform, mimicking an unfamiliar and dominant individual. Hormonal responses to stress exposure treatments were not measured since blood collection is likely to interfere with peripheral thermal profiles; however, behavioural responses to each stressor were visually and statistically confirmed (alarm calling, panting and elicitation of avoidance behaviour during stress exposure: J.K.R., personal obs.; reduction in feeding: see Appendix and Fig. A1). Individuals within control treatments were maintained according to acclimation conditions and were not subject to handling or disturbance by experimenters.

At the onset and completion of each treatment (control and stress exposure), all individuals were re-weighed to monitor changes in body condition across treatments (see Appendix), and upon completion of the experiment, individuals were released to their site of capture.

We were interested in testing the effects of repeated stress induction on surface temperature profiles, within the context of naturally cycling temperatures. We therefore sampled T_a (°C) at the location of experimentation using a ThermoChron iButton™ (Maxim Integrated, DS1922L-F5, San Jose, CA, USA) placed in the shade, for the duration of experimentation. To capture rapid and subtle changes in T_a across time, T_a readings were sampled at a frequency of 20 samples h⁻¹, and at a resolution of 0.5°C. Because relative humidity can influence the transmission of infrared radiation, and therefore estimates of object temperature by thermography (reviewed in Tattersall, 2016), we collected local relative humidity in addition to T_a (one sample h⁻¹) from Environment Canada climate repositories (https://climate.weather.gc.ca/; Hamilton A, 22 km from the experimental holding location).

Surface temperature responses to stress exposure and control treatments were measured by capturing maximum eye region temperature of chickadees during feeding, using time-lapse infrared thermography. Temperature of the eye region was assessed because it contains exposed integument that may be readily imaged in birds while stationary (e.g. during feeding) and is capable of heat exchange unfettered by insulatory keratinous tissues (i.e. feathers or leg scale; discussed in Jerem et al., 2018). Furthermore, in domestic chickens, exposure to hyperthemic conditions has been shown to increase blood flow in capillary beds of the head (Wolfenson et al., 1981), suggesting that cephalic vasculature, including the vessels located near the eye region (i.e. the opthalmotemporal, ethmoid and facial veins), may serve as a location for heat exchange (discussed in Midtgård, 1983). Here, we chose to use the maximum temperature of the eye region as a metric of eye region temperature, rather than mean temperature, because it is thought to be less susceptible to measurement error (Jerem et al., 2015, 2018) and is less likely to fluctuate according to the angle at which an individual was imaged.

Thermographic images were captured using a remotely operated thermographic camera (VueProR™, FLIR, Wilsonville, OR, USA; 13 mm lens, 336×256 resolution; image frequency=1 Hz) that was rotated among flight enclosures each day, according to cardinal direction. Specifically, beginning at 08:00 h each day, time-lapse thermographic imaging was conducted for approximately 1 h at one flight enclosure, after which the thermographic camera was transferred cardinally clockwise to a second flight enclosure, and imaging was repeated for 1 h. This rotational process was repeated until 16:00 h each day, to ensure that each flight enclosure was subjected to at least 1 h of thermographic filming per day. To control for possible effects of circadian rhythms on surface temperature profiles, the first flight enclosure to be filmed each day was also rotated cardinally clockwise direction. Because individuals within a flight enclosure could not be identified by thermography (according to band colour combinations), a remotely operated digital camera was rotated alongside the thermographic camera, allowing for post hoc individual identification.

To ensure that chickadees were blind to the presence of both the thermographic and digital camera during experimentation, flight enclosures were equipped with water-tight camera boxes that were mounted to an exterior wall adjacent to the feeding platforms, and were perforated with two 30 mm diameter holes through which thermographic imaging and digital filming were conducted (distance of 0.5 m from feeding platforms). Camera boxes were solely accessed from the exterior of flight enclosures, where an experimenter could not be seen by chickadees within the flight enclosure. When thermographic imaging was not being conducted, both 30 mm holes were covered.

Throughout experimentation, 1,035,512 thermographic images were captured across 60 days. Raw radiance values (kW m⁻²) per pixel were extracted from all thermographic images in R (v.3.6.1; https://www.r-project.org/), then converted to temperature values (°C) using Planck's law, and following equations described by Minkina and Dudzik (2009) and Tattersall (2016). Here, temperature values were calibrated according to T_a and relative humidity (determined as described above), and calibration and atmospheric constants for our thermographic camera were identified using Exiftool (https://exiftool.org/). Emissivity of the eye region was assumed to be 0.95, according to estimates for avian integument by Best and Fowler (1981).

Following conversion of radiance to temperature in thermographic images, we used FIJI (https://imagej.net/Fiji) to determine maximum eye region temperature (°C) for chickadees that were present in thermographic images by drawing a region of interest (ROI) around the eye region (∼1.1 cm in diameter) and extracting the maximum temperature from within the ROI (∼230 pixels; Fig. A2). Because motion of an object can lead to underestimation of its surface temperature in thermographic imaging (discussed in Jerem et al., 2018), only images where the feeding individual was stationary (i.e. not in flight, landing from flight or departing) were included in our analyses. Furthermore, maximum eye region temperature was only selected from images where the identity of the individual imaged could be determined from parallel digital video (n=6431; n_stress=3397, n_control=3034).

The thermoprotective hypothesis predicts that stress-induced changes in surface temperature alleviate thermal burdens incurred during a stress response, with stressed individuals conserving more heat at low T_a and dissipating more heat at high T_a than control individuals. To test this prediction, we sought to quantify and compare total heat transfer (q_total) from the eye region of captive chickadee, in stress-induced and control treatments. Here, heat transfer was calculated using maximum eye temperature measurements that were extracted from thermographic images, and by following methodology described by Ward et al. (1999), McCafferty et al. (2011) and Nord and Nilsson (2019), with slight modifications (see Appendix). Because individuals were sheltered from wind in flight enclosures, and were unlikely to transfer heat by direct contact between the eye region and a medium other than air, q_total was assumed to be the sum of radiative heat transfer (q_rad.) and free convective heat transfer (q_conv.). Here, the expected values for the kinematic viscosity of air (m² s⁻¹; at an assumed atmospheric pressure of 101.325 kPa), thermal conductivity of air (W m⁻¹ °C⁻¹) and thermal expansion coefficient of air (1 K⁻¹) used in the calculation of q_rad. and q_conv. were calculated according to T_a at the time of image capture (see Appendix). Eye region surfaces were treated as planar structures, similar to the ventral surfaces of blue tits (Parus caeruleus) described in Nord and Nilsson (2019), and the surface area of the imaged eye region was estimated to be an oval of 1.0 cm vertical diameter and 1.1 cm horizontal diameter. Final q_total values were multiplied by two to represent total heat transfer across both eye regions of an individual.

Sex of black-capped chickadees was determined genetically according to Griffiths et al. (1996) and Fridolfsson and Ellegren (1999). Briefly, whole DNA was isolated from lysed erythrocyte samples by phenol:chloroform:isoamyl alcohol (25:24:1; Thermo Fisher Scientific, Waltham, MA, USA) extraction and precipitation in 100% 2-propanol, then stored at −20°C. Following DNA isolation, sex of individuals was determined by PCR amplification of chromohelicase DNA binding protein intron 16 (Griffiths et al., 1996; Fridolfsson and Ellegren, 1999), and size separation of amplicons on 3% agarose gels.

All statistical analyses were conducted using R and α levels were set to 0.05. All generalized additive mixed-effects models (GAMMs) were fitted with restricted maximum likelihood (REML), and model assumptions (i.e. normality, homoskedasticity and temporal independence of residuals) were assessed visually. Marginal means used in plots were calculated using the R package ‘emmeans’ (https://cran.r-project.org/web/packages/emmeans/index.html).

We were interested in testing whether changes in surface temperature following stress exposure serve a thermoregulatory function. We therefore asked whether individuals exposed to repeated stressors conserved more heat below thermoneutrality and dissipated more heat above thermoneutrality, than those exposed to control treatments. Because measurements of dry heat transfer are likely to be strongly correlated with T_a in endotherms (e.g. Simmons et al., 1997), we first controlled for the global effects of T_a on q_total by calculating how different an individual's q_total was from that expected at a given T_a. This was accomplished by quantifying residual distance between observed q_total values (here, mean q_total per hour and per day, for each individual; n_obs=1006 across n_days=60 and n_birds=20), and those predicted by T_a (mean per respective hour) using a GAMM with q_total as the response variable and mean T_a per hour as a cubic regression spline (restricted to four knots) in mgcv (Wood, 2011) [mean T_a: F=1.456×10⁴, estimated degrees of freedom (edf)=3.000, P<0.0001].

According to our thermoprotective hypothesis, q_total (heat lost to the environment) from an individual under stress should decrease when T_a is below thermoneutrality and increase above thermoneutrality. We therefore divided T_a measurements into temperature zones, according to their position with respect to thermoneutrality; specifically, below thermoneutrality (‘low’), thermoneutral (‘mid’), or above thermoneutrality (‘high’). Because our experiment spanned late winter and early summer (early April to late June; T_a range: 2.5–38.5°C), limits of thermoneutrality in black-capped chickadees were generously set to 14°C and 30°C, according to Grossman and West (1977) and Rising and Hudson (1974). The effect of stress exposure on q_total across T_a zones was then tested using a linear mixed effect model (LMM) in the R package ‘lme4’ (https://cran.r-project.org/web/packages/lme4/index.html), with residual q_total as the response variable. Treatment, temperature zone (factorial: low, mid or high; total observations: n_low=110, n_mid=710, n_high=185), and an interaction between treatment and temperature zone were included in our model as fixed effects. To test whether sex influenced the direction and magnitude of q_total across T_a and stress treatments (as predicted by the thermoprotective hypothesis), sex and a three-way interaction between sex, temperature zone and treatment were initially included as fixed predictors in our model. Neither sex nor the three-way interaction between sex, temperature zone and treatment were significantly correlated with q_total (sex: P=0.068; interaction between treatment, sex and temperature zone: P=0.146), and were subsequently removed to test the effects of temperature zone and treatment on q_total alone, with increased statistical power (Aiken and West, 1991). Finally, individual identity, flight enclosure and date were included as random intercepts to account for residual variance explained by each.

Eye region temperature significantly increased with T_a (P<0.0001; Table 1) and this relationship was significantly influenced by treatment (T_a×Treatment: P=0.0195; Table 1; Fig. 1), as predicted by both the thermoprotective and haemoprotective hypotheses. The directionality of surface temperature responses to stress, however, supported predictions of the thermoprotective hypothesis alone, with eye region temperature decreasing at low T_a and increasing at high T_a with respect to control treatments (Fig. 1). Interestingly, treatment alone (regardless of T_a) did not have a significant effect on eye region temperature (P=0.995; Table 1).

Although we did not detect circadian changes in eye region temperature within our observed time period (Time of day; P=0.840; Table 1), treatment significantly influenced the relationship between time of day and eye region temperature in our birds (P=0.026; Table 1). Both date of image capture and flight enclosure explained considerable variance in eye region temperature (s.e.m._date=0.425; s.e.m._enclosure=1.819; Table 1), probably owing to differing degrees of exposure to solar radiation per day in each enclosure.

Mean±s.d. q_total from the eye region was 25.275±11.903 mW and residual q_total from the eye region differed significantly between temperature categories (Table 2; Fig. 2), as predicted by the thermoprotective hypothesis. Specifically, individuals experiencing stress exposure conserved significantly more heat from the eye region than controls, when held below thermoneutral temperatures (P=0.010; Table 2, Fig. 2) and lost significantly more heat from the eye region than control birds, when held above thermoneutral temperatures (P=0.003; Table 2; Fig. 2). At thermoneutral environments, however, stress-exposed and control individuals did not differ with respect to heat exchange at the eye region (P=0.064; Table 2; Fig. 2).

Females and males significantly differed in their responses to stress exposure across T_a, as predicted by the thermoprotective hypothesis but not the haemoprotective hypothesis. Specifically, eye region temperature in males displayed a more robust response to stress exposure than that of females across T_a (n_females=10, n_males=10; β_males=0.038, s.e.m.=0.017, t=2.269, d.f.=17, P=0.037; β values are means±s.e.m.), with the slope between T_a and eye region temperature being steeper following stress exposure in males than females (Fig. 3). Contrary to predictions of the thermoprotective hypothesis, however, females and males did not significantly differ in the magnitude of their response to stress exposure at low T_a (0°C; as measured by the intercept of their eye region temperature by T_a curve, under stress exposure; β_males=0.219±0.257, t=0.850, d.f.=17, P=0.407; Fig. 3). Mass was not significantly correlated with the slope or intercept of the relationship between eye region temperature and T_a under stress exposure (slope: β_mass=−0.009±0.010, t=−0.916, d.f.=17, P=0.373; intercept: β_mass=−0.201±0.151, t=−1.327, d.f.=17, P=0.202). Interestingly, q_total from the eye region did not reflect patterns that were shown for eye region temperature. A global three-way interaction between sex, treatment and temperature zone was not significant (F=1.928, d.f.=799, P=0.146), demonstrating that we could not detect sex-specific changes in heat transfer across temperature zones and treatment, probably owing to our small sample size (n_females=10, n_males=10).

Our results provide evidence for the thermoprotective hypothesis over the haemoprotective hypothesis and lend support to direct functional significance of stress-induced changes in surface temperature. Indeed, as predicted by the thermoprotective hypothesis, we report a significant interaction between T_a and treatment on eye region temperature in our experimental population, where eye region temperature in stressed individuals was lower at low T_a, and higher at high T_a than that of control individuals (Fig. 1). Furthermore, individuals exposed to repeated stressors conserved more heat at temperatures below thermoneutrality and dissipated more heat at temperatures above thermoneutrality with respect to controls (Fig. 2). These results demonstrate that changes in surface temperature experienced during stress exposure provide a thermoregulatory benefit across ecologically relevant temperature gradients, and likely contribute to lessening thermal burdens experienced during activation of a stress response. In our study, the effect of treatment on eye region temperature across T_a could not be explained by feeding behaviour or individual condition alone (Figs A3–A5), supporting a direct physiological modification of surface temperature in response to stress exposure.

Although our reported patterns of q_total under stress exposure are suggestive of a thermoregulatory adaptation, we recognize that our observed differences in q_total between stress-exposed and control individuals are small (mean Δq_total=2.732 mW or 0.937 J g⁻¹ h⁻¹; Fig. 2). Indeed, for a stress-exposed individual of average mass within our study population (10.5 g), energetic savings by such changes in q_total alone total ∼1% of daily energetic expenditure, as estimated according to Karasov et al. (1992). Given we assessed q_total across a small area of surface tissue alone (the eye region), however, our estimates of energy savings likely underestimate true energy savings associated with whole-body changes in surface temperature following stress exposure. In European starlings (Sturnus vulgaris), for example, most dry heat transfer occurs at the ventral brachial area and legs (Ward et al., 1999), while in the toco toucan (Ramphastos toco), dry heat transfer is predominantly sanctioned to the bill (Tattersall et al., 2009). Combined heat transfer experienced across the legs, bill and eye region in our birds almost certainly exceeded that experienced across the eye region alone, if thermal responses to stress exposure at the legs and bill mimic those observed at the eye region (but see Bech and Midtgård, 1981). Indeed, Herborn et al. (2015) reported a concomitant change in wattle, comb and eye temperature following stress exposure in the domestic chicken, demonstrating that surface tissues are unlikely to respond in isolation to stress exposure. In this study, however, we were unable to estimate q_total across the legs and bill of individuals because each structure could not be readily observed within our thermographic images (see Fig. A2). Furthermore, without calorimetric measurements of our chickadees, the precise energetic merits of our observed changes in q_total following stress exposure remain in question.

Beyond providing direct thermoregulatory advantages, stress-induced changes in surface temperature may also yield indirect advantages with respect to water retention. In small birds, active heat dissipation is predominantly achieved by evaporative cooling at cutaneous tissues in the respiratory tract (via panting; Tieleman and Williams, 2002). Cooling by this mechanism, however, can provide significant water loss that may be challenging to sustain. Black-capped chickadees, for example, have been shown to lose up to 400 mg of water h⁻¹ by evaporative cooling at T_a above thermoneutrality (4% of their body mass; Rising and Hudson, 1974). Estimates for other temperate species are comparable, with red breasted nuthatches (Sitta canadensis) losing up to 260 mg of water h⁻¹ by evaporative cooling at high T_a (2.6% of their body mass; Mugaas and Templeton, 1970). As evaporative heat loss at surface tissues is thought to be negligible in contrast to respiratory tissues (Bernstein, 1971), heat dissipation at surface tissues by radiation or conduction therefore provides significant water-retentive advantages when compared with cooling by panting. In desert-dwelling birds, for example, total evaporative water loss may be reduced by 50% when employing methods of dry heat exchange (i.e. peripheral hyperthermia) rather than evaporative cooling (Tieleman and Williams, 1999). Under stressful conditions, water stores may already be compromised (i.e. by stress-induced defecation; Jones et al., 1995; Haas et al., 2010) and the ability to seek and acquire water may be over-ruled by combative or avoidance behaviours. For this reason, elevations in heat transfer from surface tissues of stress-exposed chickadees could also reflect the outcome of a thermoregulatory trade-off at high T_a, where dissipation of heat generated during a stress response, albeit small, is balanced by a pressure to retain water stores.

Our results regarding sex-specific thermal responses to stress exposure provide mixed support for the thermoprotective hypothesis but little to no support for the haemoprotective hypothesis. In our study, the effect of stress exposure on the relationship between eye region temperature and T_a was significantly greater in males than in females (Fig. 3; slope of the interaction between treatment and T_a). This trend was consistent with predictions of the thermoprotective hypothesis, as male chickadees are thought to have lower heat dissipation capacity than females, owing to their larger body size (Foote et al., 2010; Aschoff, 1981; Porter and Kearney, 2009) and lack of brood-patch development (Odum, 1941; Cooper and Voss, 2013). Contrasting the thermoprotective hypothesis and supporting the haemoprotective hypothesis, however, we did not detect a significant effect of sex on the magnitude of stress-induced changes in eye region temperature at low T_a (i.e. the intercept of the relationship between T_a and eye region temperature; Fig. 3), although females were trending toward a larger stress-induced fall in eye region temperature than males at low T_a (Fig. 3). While such similarity between sexes suggests that stress-induced changes in surface temperature do not occur in proportion to metabolic costs in cold environments, these results may more simply reflect an inability of females to further increase heat conservation by peripheral vasoconstriction at such low observed temperatures (as low as ∼2°C in our study). In domestic ducks (Ana boscas), for example, local peripheral vasoconstriction is used to minimize heat loss in the feet during cold exposure, but only until subcutaneous temperatures reach a critically low level (∼8°C; Reite et al., 1977). Similar thermal limits to vasoconstriction are likely to exist for vasculature surrounding the eye region (i.e. the opthalmotemporal, ethmoid and facial veins), and indeed probably exceed that of vasculature in the feet given their close proximity to the brain (Midtgård, 1983). To our surprise, our analyses did not detect differences in the effect of stress exposure on q_total between sexes, however, a low sample size (n_females=10, n_males=10) probably contributed to a lack of statistical significance (P=0.146). To our knowledge, no other studies have tested the effect of sex differences in surface temperature responses to stress exposure.

Interestingly, mass did not explain individual variation in the intercept (P=0.373) or slope (P=0.202) of the T_a by treatment interaction in our sample population. While it is possible that variance in mass within our sample population was too low to bestow meaningful influences on heat dissipation capacity (s.d.=0.841), it is more probable that differences in thermal windows between sexes overwhelm mass effects in our experiment; indeed, given our experiment was initiated at the onset of the breeding season for black-capped chickadees (April; Foote et al., 2010), our observed differences in surface temperature under stress exposure are more likely a consequence of females alone displaying brood-patch development. Female brood patch development begins in April at a comparable latitude to our field location (Odum, 1941; Cooper and Voss, 2013); approximately the same time-point at which our experiments were initiated. Although we did not observe abdominal defeathering in our birds during experimentation, females may have displayed greater abdominal vascularization inherent to brood patch development than males (Hinde, 1961; Etche et al., 1979; but see Bailey, 1952), allowing them enhanced heat dissipation capacity at high temperatures (Hill et al., 2014; Nilsson and Nord, 2018), and levelling the slope between T_a and eye region temperature in stress treatments, particularly above thermoneutrality.

In addition to differences in thermal windows between sexes, a more robust change in surface temperature following stress exposure in males could reflect a stronger selective pressure for thermoregulation in this sex, particularly when both T_a and metabolic heat production are high (i.e. during activation of a stress response in warm environments). Differences in thermal fertility limits between sexes have recently been argued as an emerging trend across species, with males typically displaying higher thermal sensitivity than females (Greenwood and Wheeler, 1985; Iossa, 2019). Such sex differences may be exacerbated in taxa where males have internal testes (i.e. birds). Sex-specific thermal limits to fertility could exert a stronger selection pressure for thermoregulatory capacity, or heat dissipation capacity, in male chickadees at high T_a when compared with female chickadees. Higher eye region temperature in male chickadees exposed to both high T_a and perceived stressors may therefore be explained by adaptations to dissipate more heat (that is generated by heat exchange with the environment, and as a byproduct of mounting a physiological and behavioural stress response) at the level of surface tissues when compared with females.

Beyond differences in thermal windows and heat dissipation capacity, the disparity in stress-induced changes in surface temperature between sexes may also be explained by differences in autonomic and steroidal pathways underpinning the thermal response to stress. In birds, males typically display a more robust activation of the hypothalamic–pituitary–adrenal (HPA) axis following stress exposure (Silverin, 1998; Marin et al., 2002; Madison et al., 2008; Wada et al., 2008). Across vertebrates, HPA axis activation is thought to both enhance sympathetic function (Fisher et al., 1982) and sensitize the body to catecholamine availability (Fisher et al., 1982; Krakoff, 1988), thereby exacerbating the tachycardic, tachypneic and haemodynamic response to a perceived stressors (reviewed in Sapolsky et al., 2000). In our study, differences in thermal response curves between sexes (across T_a and treatments) may simply be explained by exacerbated vasoconstriction at the periphery in males below thermoneutrality (Fig. 3), and elevated total heat production in males above thermoneutrality (Fig. 3), each as a consequence of enhanced physiological stress-responsiveness (i.e. tachycardia, tachypnea and avoidance behaviour; Cabanac and Aizawa, 2000; Cabanac and Guillemette, 2001; Greenacre and Lusby, 2004; Long et al., 1990a).

Here, we proposed a new hypothesis stating that changes in surface temperature following exposure to stressors hold direct functional significance by offsetting thermoregulatory costs experienced during activation of a stress response (the thermoprotective hypothesis). This hypothesis contrasts with the dominant theory, which argues that changes in surface temperature following stress exposure are a mere consequence of haemetic redistribution and hold no direct functional significance in themselves. Although our results do not preclude a selective influence of haemorrhage avoidance on stress-induced changes in surface temperature, they do provide support for a functional significance of stress-induced changes in surface temperature, per se, with respect to thermoregulation. Indeed, we show that in a temperate endotherm, stress-induced changes in surface temperature influence dry heat transfer to the environment and allow individuals to conserve more heat at low T_a, and dissipate more heat at high T_a, when experiencing prolonged stress exposure. While such changes in heat conservation and dissipation are suggestive of energetic merits associated with stress-induced changes in surface temperature, direct metabolic consequences of this response remain in question.

Overall, these findings emphasize a seldom-recognized influence of the thermal environment in shaping physiological responses to stress exposure in vertebrates and offer new explanations about the ultimate mechanisms driving a long-recognized physiological phenomenon. In doing so, however, our study raises questions about the suitability of stress-induced thermal flux within the context of a changing global climate, whereby occurrence of extreme heat-waves may negate the thermoregulatory benefits of these stress-induced responses.

The average diameter of the eye region in our chickadees was estimated as 1.1×10⁻² m, and the area of the eye region was assumed to be an ovoid with a vertical diameter of 1.0×10⁻² m.

Exposure to psychogenic stressors is capable of influencing feeding behaviour (Favreau-Peigné et al., 2014), with psychologically stressed individuals feeding either less or more frequently than unstressed individuals (Oliver and Wardle, 1999). To confirm whether repeated stressors applied during stress-induced treatments were sufficient to elicit a behavioural response in our chickadees, we tested whether feeding behaviour differed between control and stress-induced treatments. Because feeding behaviour can be altered by either shifting feeding rate and/or shifting total time spent feeding, we combined these two metrics using a scaled and centred principal component analysis (PCA). As PC1 explained 93.168% of variance in our data and was positively correlated with both loading variables, we used PC1 as our metric for feeding behaviour (termed ‘feeding score’) in subsequent analyses. To aid in visual interpretation, feeding score was normalized between 0 and 100, such that individuals with a high feeding score were observed to feed more frequently and remain feeding for longer periods of time at each visit when compared with those with a low feeding score.

The effect of stress exposure on feeding behaviour was modelled using a generalized additive mixed model (GAMM) in the R package ‘mgcv’ (https://cran.r-project.org/web/packages/mgcv/index.html), with feeding score as the response variable, treatment as a binomial parametric predictor, T_a (°C, mean h⁻¹) and time of day (h) as cubic regression splines, and an interaction term between treatment and each cubic spline (T_a and time of day). Splines were restricted to five knots to circumvent model over-fit while capturing the predicted curvilinear relationship between smooth predictors and the response. Individual identification, flight enclosure and date were included as random intercepts to account for variance in feeding behaviour explained by each parameter. A negative-binomial distribution with log-link and a θ of 2.6 was assumed to account for over-dispersion and Poisson distributed residuals.

Results of our additive model show that feeding score (PC1) significantly differed between treatments, with individuals exposed to rotational stressors feeding less than those experiencing control treatments (n_obs=1014, β=−0.068±0.033, z=−2.089, P=0.037; Fig. A1). This effect of treatment on feeding score was not dependent on T_a (n_obs=1014, edf=1.864, s.e.m.=0.166, χ²=125.012, P=0.069), nor did T_a alone correlate with feeding score (n_obs=1014, edf=3.435, s.e.m.=0.192, χ²=62.187, P=0.403). Time of day, however, was significantly correlated with feeding score (n_obs=1014, edf=3.778, s.e.m.=0.109, χ²=36.678, P=0.012), and this relationship significantly differed between treatments, with individuals in control treatments feeding more in morning than those in stress-induced treatments (n_obs=1014, edf=0.951, s.e.m.=0.072, χ²=7.811, P=0.017; Fig. A1). Our GAMM explained 37.7% of the deviance in our data.

Exposure to repeated stressors significantly influenced feeding behaviour of black-capped chickadees in our experiment (see Table 2). Because both feeding behaviour and the metabolic consequences of movement implicit in foraging are sufficient to modulate body temperature (see Zhou and Yamamoto, 1997; Sarmiento-Franco et al., 2000; Nord and Nilsson, 2019), we tested whether the effect of repeated stressors on feeding behaviour was sufficient to explain our observed differences in eye temperature across T_a between treatment groups (see Table 1). To do so, we asked whether eye temperature was correlated with feeding score, when controlling for the influence of T_a on both parameters. We addressed this question by averaging eye temperature for each individual per hour, per day, then calculating average T_a during these respective time period. We then quantified the relative difference between observed eye temperature means and those expected at any given T_a (regardless of treatment) by constructing a simple GAM in the R package ‘mgcv’ with mean eye temperature as the response variable, and mean T_a as a smoothed predictor (cubic regression spline with four knots), then extracting residual eye temperature from this GAM. Our previous results (see Table 1) show that the direction of the effect of treatment on eye temperature varies according to T_a. We therefore calculated the absolute values of the eye temperature residuals to create a directionless metric of distance between observed eye temperature means, and expected eye temperature means for a given T_a.

To test whether differences in feeding rate, but not treatment, could explain discrepancies between mean eye temperature for a given individual and that expected for a respective T_a, we constructed a generalized linear mixed effects model [GLMM; R package ‘cplm’ (https://cran.r-project.org/web/packages/cplm/index.html); Zhang, 2013] with residual eye temperature (absolute values, as described above) as the response variable, and feeding score and treatment (binomial factor) as fixed effect predictors. Time (h) was also included is a fixed effect to account for circadian rhythms in body temperature, as well as an interaction between time and treatment to account for significant interaction between each variable observed in our previous model (see Table 1). Residual effects of individual identity, flight enclosure identity, and date were initially accounted for by including each parameter as a random intercept, however, individual identity explained zero variance in our model and was therefore removed. A Tweedie distribution was assumed for our model as our response variable was right-skewed, non-negative, and a non-integer (Tweedie, 1984). Our final model included 581 observations across 50 days, with all four flight enclosures considered.

Results of our GLMM show that feeding score was not significantly correlated with residual eye temperature (β=−0.001, t=−0.491, d.f.=581, P=0.312; Fig. A3) but was significantly influenced by treatment, with stressed individuals displaying high residual eye temperature value (β=1.444, t=2.811, d.f.=581, P=0.003; Fig. A4). Similar to our previous results (Table 1), time alone was not significantly correlated with eye temperature (β=0.041, t=1.105, d.f.=581, P=0.135), however, a significant interaction between time and treatment was observed (β=−0.137, t=−2.794, d.f.=581, P=0.003).

A positive relationship between body condition and surface temperature at the eye has been reported in blue tits (Jerem et al., 2018). We therefore tested whether repeated stress exposure was sufficient to decrease body condition, and therefore modulate eye temperature in our study system. To do so, we quantified individual body condition by regressing mass (g) against wing chord (mm) in a linear model and extracting model residuals. Mass measurement at each time intervals were not available for all individuals (n=16 of 20), and individuals without mass measurements at each time interval were removed from analyses (n=4). Residuals were then normalized between 0 and 1 to aid in visual interpretation and are hereafter referred to as ‘body condition’. First, we tested whether body condition of each individual changed between treatments using two paired Student's t-tests. Next, we quantified the total change in body condition (‘Δbody condition’) across each treatment by subtracting each individual's body condition after a given treatment by their respective body condition before the same treatment. The effect of treatment on condition was then tested using an LMM with Δbody condition as the response variable, treatment as a factorial fixed effect and individual identification and flight enclosure as random intercepts. Individual identity and flight enclosure ID, however, were later removed from our model as they explained zero residual variance. Sex and the interaction Treatment×Sex were included to test for sex-specific responses in condition to treatment regimes. The effect of condition on eye temperature alone, however, was not assessed because condition metrics were only available for two time points per treatment (beginning and end), thereby limiting statistical power (π=0.56 at an estimated Cohen's d of 0.50).

Interestingly, we did not detect statistical support for a change in condition across stress treatments (t=−0.959, d.f.=15, P=0.352, n=16) or control treatments (t=−1.008, d.f.=15, P=0.330, n=16; Fig. A5). Similarly, the total change in condition scores across a treatment did not statistically differ between treatment types (n_obs=32, β=0.035, t=0.611, P=0.546; Fig. A5), sexes (n_obs=32, β=0.069, t=1.040, P=0.307; Fig. A4) or an interaction between both factors (n_obs=32, β=−0.058, t=−0.623, P=0.538).

We thank Lianne Ralph for exceptional patience and field assistance, Glenn Tattersall for thermographic advice, Simon Tapper and Robby Marrotte for statistical input, Jan-Åke Nilsson and Andreas Nord for theoretical conversations, James Quinn for his generous contribution of laboratory space, and Kimberley Tasker and Tyler Maksymiw for their assistance in aviary construction. Additionally, we thank the staff and avian banding team at the Ruthven Park National Historic Site for their immense support in experimental procedures.