Inter- and intra-specific variation in hair cortisol concentrations of Neotropical bats

Abstract Quantifying hair cortisol has become popular in wildlife ecology for its practical advantages for evaluating stress. Before hair cortisol levels can be reliably interpreted, however, it is key to first understand the intrinsic factors explaining intra- and inter-specific variation. Bats are an ecologically diverse group of mammals that allow studying such variation. Given that many bat species are threatened or have declining populations in parts of their range, minimally invasive tools for monitoring colony health and identifying cryptic stressors are needed to efficiently direct conservation efforts. Here we describe intra- and inter-specific sources of variation in hair cortisol levels in 18 Neotropical bat species from Belize and Mexico. We found that fecundity is an important ecological trait explaining inter-specific variation in bat hair cortisol. Other ecological variables such as colony size, roost durability and basal metabolic rate did not explain hair cortisol variation among species. At the individual level, females exhibited higher hair cortisol levels than males and the effect of body mass varied among species. Overall, our findings help validate and accurately apply hair cortisol as a monitoring tool in free-ranging bats.


Introduction
Free-living animals face multiple natural and anthropogenic challenges that threaten their survival and thus are of considerable interest to ecophysiologists concerned with the study of effects of stress on vertebrates. One of the most extensively studied processes associated with response to stressors (biotic or abiotic environmental factors that disrupt homeostasis; Schulte, 2014) is the release of glucocorticoid (GC) hormones (Creagh and Brendan Delehanty, 2013;MacDougall-Shackle-ton et al., 2019). GCs are known to facilitate the mobilization of energy required to cope with stressors and, during normal conditions, play a key role in regulating growth, circadian activity and energy metabolism (review in Landys et al., 2006). Levels of GCs are commonly employed as a biomarker of allostatic load or stress (indirect indicators of health) (Sapolsky et al., 2000;Wikelski and Cooke, 2006;Pearson Murphy, 2007;Busch and Hayward, 2009). GC secretion is a well-conserved process across vertebrates and involves activation of the hypothalamic-pituitary-adrenal (HPA) axis Although many of the environmental challenges that wild populations experience are chronic (e.g. prolonged food deprivation, climate change, habitat disturbance, pollution), studies of stress physiology have focused on detecting acute stress by looking at GC levels in blood, urine and faeces (Sheriff et al., 2011;Creagh and Brendan Delehanty, 2013). The rapid turnover of these tissues, however, only gives shortterm information of HPA activity over periods of hours or days (Sheriff et al., 2011), which may not be an appropriate time scale. Assessment of cortisol in tissues with slower turnover rates, such as hair, may reflect circulating cortisol levels over longer periods of several weeks or even months, which is the time scale over which chronic environmentally induced stress would be expected to occur (Davenport et al., 2006;Macbeth et al., 2010;Ashley et al., 2011;Mastromonaco et al., 2014). Cortisol is incorporated into developing hairs from the blood stream during periods of active hair growth, allowing researchers to retrospectively examine cortisol production at the time that a stressor or stressors were faced (Davenport et al., 2006;Pragst and Balikova, 2006). Hair can be collected in a minimally invasive manner, is usually easily accessible in relatively large amounts and is easy to store and transport, all of which make it particularly useful for wildlife studies, especially those involving threatened or endangered species (Koren et al., 2002;Macbeth et al., 2010Macbeth et al., , 2012. Hair cortisol levels are not likely affected by stress induced by capture and/or handling, which is one of the main limitations of blood GC analysis (Russell et al., 2012). A single sample of hair can also provide complementary and valuable information about ecology and behaviour, including diet and movement (e.g. using stable isotope analyses; Fraser et al., 2010;Sullivan et al., 2012;Voigt et al., 2012;Oelbaum et al., 2019), condition (e.g. nutrition; Montillo et al., 2019), toxicant exposure (Hernout et al., 2016;Becker et al., 2018) and molecular identification (Magioli et al., 2019), opening possibilities for more integrative studies. However, analyses of hair samples can be challenging. Despite being a very promising tool for assessing wildlife health, quantifying hair cortisol is a method that has limitations; although these are largely based on lack of detailed knowledge of patterns of hair growth (Meyer and Novak, 2012;Russell et al., 2012;Sharpley et al., 2012). For example, the exact time scale reflected in any given sample will depend on the rate of hair growth and moulting patterns; this information is unknown for most species, which makes the time window being evaluated unclear (Koren et al., 2002;Fourie et al., 2016). Moreover, rates of cortisol incorporation to the hair shaft are known to differ across body regions and among species (Sharpley et al., 2012;Acker et al., 2018;Lavergne et al., 2020). Nevertheless, hair cortisol levels offer a potentially powerful tool for assessing relatively long-term stress levels in mammals.
Hair cortisol and its correlation with natural and anthropogenic stressors has been explored for different wild mammals, including rhesus monkeys (Macaca mulatta; Dettmer et al., 2014), grizzly bears (Ursus arctos; Macbeth et al., 2010), reindeer/caribou (Rangifer tarandus; Ashley et al., 2011), lynx (Lynx canadensis; Terwissen et al., 2013;Azevedo et al., 2020), mongoose (Herpestes ichneumon; Azevedo et al., 2019) and snowshoe hares (Lepus americanus; Lavergne et al., 2020); other examples reviewed by Kalliokoski et al., 2019). Although most of these studies support hair cortisol as an informative measure of central HPA activity, they also identified intrinsic factors such as age, sex, reproductive stage and social status that modulate GC levels in different contexts (Wingfield and Romero, 2011;Crespi et al., 2013;Hau et al., 2016). Not accounting for these intrinsic sources of variation in GC levels may lead to incorrect or misleading estimates of the effects of stressors on individual fitness and population health (Sapolsky et al., 2000;Reeder and Kramer, 2005;Busch and Hayward, 2009;Wingfield and Romero, 2011;Kalliokoski et al., 2019).
Ecological traits such as diet, fecundity and lifespan, as well as phylogenetic relatedness, have been proposed to explain differences in baseline cortisol levels in wild species (Wingfield and Romero, 2011;Patterson et al., 2014). Evolution of different life-history strategies are also thought to have led to different adaptations in HPA activity modulation so as to maximize individual fitness within species (Bonier et al., 2009;Bonier and Martin, 2016). Bats are a very ecologically diverse group comprising over 1400 species that live in most terrestrial ecosystems and have a wide variety of diets, use many different roost types and have many different social systems (Kunz and Fenton, 2005;Dumont et al., 2012;Gunnell and Simmons, 2012;Simmons and Cirranello, 2020). This diversity provides the opportunity to study the ecological correlates of cortisol levels among phylogenetically related species with different life history traits. Few ecological correlates of GCs have been evaluated simultaneously in mammalian groups in the context of cortisol studies, and fewer studies have further related cortisol levels to life history traits across multiple species from a single mammalian clade. Among bats, variation in plasma cortisol levels associated with seasonal food availability has been studied in two species with contrasting diets, Carollia perspicillata and Desmodus rotundus (Lewanzik et al., 2012), but no other comparative studies have been conducted within this order. Furthermore, little is known about the modulation of the stress response in bats, despite Chiroptera being the second-most speciose order of mammals. Bat populations are declining worldwide due to ongoing habitat destruction and land use changes, increased interaction with human environments and associated threats including wind turbine fatalities, hunting and targeted killing, pesticide exposure and emerging infectious diseases such as whitenose syndrome (O'Donnell, 2000;Mickleburgh et al., 2002;Kunz et al., 2007;Frick et al., 2010;Racey, 2013;Voigt and Kingston, 2015). Because many bat species are threatened or have declining populations in parts of their range (IUCN Red List of Threatened Species, 2020), minimally invasive tools to monitor colony health and identify cryptic stressors are critically needed to efficiently direct conservation efforts. It is essential to investigate the factors influencing baseline GCs to properly detect elevated cortisol levels due to long-term stressors.
In this study, we describe intra-and inter-specific sources of variation in baseline hair cortisol levels in bats, which contributes to better understanding the potential for hair cortisol to be an indicator of HPA activity in this taxon. We hypothesize that interspecific variation in hair cortisol of bats will be greater than intra-specific variation and that such heterogeneity will be best explained by ecological traits directly related to energy expenditure, such as basal metabolic rate (BMR), dietary guild, foraging style and roost durability. We expect that species with high energetic demands or less predictable energy acquisition (e.g. less reliable food sources) will have higher hair cortisol. Specifically, we predict the following: (i) a positive relationship between BMR and hair cortisol;(ii) bats that feed on fruit and nectar, which are energy rich and readily available, will have lower hair cortisol; (iii) bats that actively hunt prey during flight, such aerial hawkers, will have higher GC levels owing to greater energetic demands compared to gleaners that can hunt from perches (Norberg and Rayner, 1987;Fenton, 1990); and (iv) species using more ephemeral day roosts (e.g. foliage or crevices under exfoliating bark) will have higher hair cortisol than species using more stable structures (Kunz and Fenton, 2005).

Study sites
We sampled bats from northern Belize (Orange Walk District) and two locations in Mexico (Colima and Chihuahua States). In each region, we sampled sites with different levels of habitat fragmentation and agricultural intensity. We used the global Human Modification Index (HMI; Kennedy et al., 2019) as a standardized measure of disturbance, using a 5-km buffer around each collection site. The HMI is a cumulative measurement with possible values between 0 (no disturbance) and 1 (highest disturbance) that includes transportation, human settlement, agriculture, extractive activities and electric infrastructure (Kennedy et al., 2019). Sites were classified as low (0 ≤ HMI ≤ 0.10), moderate (0.10 < HMI ≤ 0.40), high (0.40 < HMI ≤ 0.70) and very high (0.70 < median HMI ≤ 1.00). At all sites, bats were captured from 18:00 to 22:00 using mist nets and from 18:00 to 5:00 using harp traps (only in Belize) set along flight paths. Bats sampled during the day were captured in their roosts, mainly caves, using hand nets. We recorded sex, size (body mass, g) and reproductive stage (females: pregnant, lactating; males: active, inactive; Kunz and Parsons, 2009). Reproductive stage was assessed by checking for the presence of scrotal testes in males (indicating that the individual was reproductively active at the time of capture) and by the evidence of pregnancy or lactation (enlarge nipples) in females (Racey, 2009 . We refer here to these locations collectively as central Mexico, and we took samples from three species in this region: Leptonycteris yerbabuenae, Macrotus waterhosii, and Pteronotus mexicanus. We also sampled bats foraging close to pecan nut croplands near the town of Jimenez, Chihuahua (northern Mexico). This region is entirely dedicated to the production of pecan nuts with thousands of squared kilometres of cultivated land (Orona Castillo et al., 2018). We visited one that farms using organic practices and another that farms with intensive use of pesticides. However, the estimated HMI was the same for the two sites (HMI = 0.49, high disturbance). We collected hair samples from three bat species (Antrozous pallidus, Tadarida brasilensis and Myotis velifer) at both northern Mexico sites.
Our field sites in Belize consisted of two forest patches of very different size located ∼10 km apart and separated by a heterogeneous, largely agricultural landscape. Lamanai Archaeological Reserve (LAR) is a protected secondary semideciduous forest of 450 ha with a high canopy and with relatively low disturbance (HMI = 0.17) (Herrera et al., 2018). In contrast, the Ka'Kabish archaeological site (KK) is a small remnant forest patch of ∼45 ha surrounded by cattle pastures and local croplands (Fig. 1). Although the landscape in Belize is apparently disturbed and highly fragmented, agricultural activity and urban development is not as intense as the field sites in Mexico, which is reflected in their moderate HMI scores (LAR: 0.17; KK: 0.18). We collected hair samples from 12 different species (Table 1)  Mammalogists (Sikes and Bryan, 2016) and were approved by the Institutional Animal Care and Use Committees of the University of Georgia (A2014 04-016-Y3-A5), University of Toronto (20012113) and American Museum of Natural History (AMNHIACUC-20180123). Fieldwork was authorized by the Belize Forest Department under permits WL/2/1/18(16) and WL/1/19(06). Sample collection in Mexico was approved under the permit #FAUT-0069.

Sample collection
We manually trimmed a hair sample (3-10 mg) using round tip curved dissection scissors from the scapular region on the back of each bat. The resulting samples were placed individually in 1-2-ml plastic tubes using flat tweezers. The dissection tools were cleaned with 70% ethanol between sampling different individuals to avoid cross contamination. The amount of hair removed from each bat depended on the hair density of each species. The hair shaft was carefully cut close to the root avoiding removing skin or follicle tissue. From pilot analyses, we determined a minimum amount of 3 mg of hair was necessary to obtain values around 50% binding on the standard curve thereby accurately estimating cortisol concentration in the sample.

Extraction and quantification of cortisol
Hair samples were processed and analysed at the endocrinology laboratory at the Toronto Zoo following methods described by Acker et al., 2018. Each hair sample was spread apart and weighed in a 7-ml glass scintillation vial. To avoid contamination with other biological fluids, all hair samples were washed with 100% methanol by vortexing in a tube for 10 s and immediately removing the methanol using a pipettor. Immediately thereafter, 100% methanol was added to each sample, at a ratio of 0.005 g/ml. Samples were then mixed for 24 hrs on a plate shaker (MBI Orbital Shaker; Montreal Biotechnologies Inc., Montreal, Quebec City, Canada). After 24 hrs the vials were centrifuged for 10 min at 2400 g. The supernatants were pipetted off into clean glass vials and dried down under air in a fume hood. The dried extracts were stored at −20 • C until analysis.
Samples were brought to room temperature prior to analysis. Reconstitution of the desiccated extracts was done by adding phosphate buffer and vortexing for 10 s. Belize samples were reconstituted neat (i.e. evaporated 150 ul and reconstituted with 150 ul), and Mexico samples were reconstituted as follows: three species were neat, two species diluted 1:5 and one species diluted 1:50 in phosphate buffer (Andreasson et al., 2015). Cortisol concentrations were determined using an enzyme immunoassay (EIA) previously described  and horseradish peroxidase dilutions were 1:10200 and 1:33400, respectively. Cortisol, rather than corticosterone, was targeted because it has been found to be the primary circulating GC in bats, with concentrations four times higher than corticosterone (reviewed in Kwiecinski and Damassa, 2000). Biochemical validation (parallelism and recovery) of the cortisol EIA was done using pooled hair extracts (see Supplementary data). We used pooled hair samples from Big Brown bats (Eptesicus fuscus) from a captive colony at McMaster University as a model for species with low cortisol concentrations (i.e. extracts reconstituted neat). For species with diluted extracts (Tadarida brasiliensis, P. mexicanus and L. yerbabuenae) there was sufficient hair for separate pools. The inter-assay coefficient of variation (CV) for high control (24% binding) and low control (60% binding) were 9.3% and 9.8%, respectively. The intra-assay CV was 8.5%. The limits of detection and quantitation were 56 pg/ml and 153 pg/ml, respectively. Results are presented as nanograms of cortisol per gram of hair.

Validation of immunoassay
Biochemical validations showed that the cortisol assay was suitable for hair.

Species ecological traits
We compiled data on ecological traits considered relevant to cortisol mobilization from previously published literature and databases. Values for traits are species-level averages and may not reflect specific values at these sites (Table 1). Data on BMR was extracted from the literature (Cruz-Neto et al., 2001;Genoud et al., 2018; see Supplementary data) and when not available (n = 2) the following formula was used for the estimation: ln BMR = 0.744 × ln mass(in g) + 1.0895 (Speakman and Thomas, 2003). Information on diet, foraging style, percentage of invertebrates in the diet and fecundity was extracted from the Elton Traits, PanTHERIA and Amniote Life History databases (Jones et al., 2009b;Wilman et al., 2014;Myhrvold et al., 2015). We collapsed variation in diet into two dietary guilds: phytophagy (including nectarivores and frugivores) and animalivory (insectivores and carnivores) because many bat species in our study have diets that combine more than one food source within these categories (Fenton et al., 2001;Kunz and Fenton, 2005;Reid, 2009;Oelbaum et al., 2019). We also considered the percentage of invertebrates in the diet of the animalivorous bats, which can vary significantly among species. Because foraging behaviour is a complex and plastic trait, we simplified this variable into two categories: aerial foragers (i.e. hawkers) and gleaners (including species that glean plant products like fruit as well as insects) since these behaviours may reflect differences in energetic demands associated with foraging (Herrera et al., 2018). Because wing morphology can strongly influence the energetic costs of flight, we also included the mean wing aspect ratio for each species (Norberg and Rayner, 1987;Bullen et al., 2014). Fecundity was defined as the annual average fecundity (litter size × number of litters per year). We estimated roost durability following the methods of Patterson et al. (2007), where 1 indicates the most ephemeral and least protected roost types (e.g. rolled leaves and foliage) and 6 indicates the most permanent and protected roost types (e.g. caves). For species known to multiple use different kinds of roost, intermediate ranks were calculated, weighting roost categories according to the relative frequency of use reported in the literature (Schneeberger et al., 2013). Lifespan was drawn from the Animal Ageing and Longevity database (AnAge: The Animal Ageing and Longevity Database, 2020) and DATLife (DATLife Database, 2020). Lifespan was grouped in five categories: 0-5, 5-10, 10-15, 15-20 and >20 years. For many of the species in these databases, longevity estimates are based on captive animals, which likely overestimates life expectancy in the wild. Because bats of a single species may live in colonies of varying sizes, and most values of colony size are reported in ranges in the literature, we classified maximum colony sizes reported for each species as small (1-50), medium (50-500) or large (>500) sensu Santana et al. (2011).

Data analysis
We first used phylogenetic generalized least squares (PGLS) models to evaluate the effect of species-level ecological variables on hair cortisol concentrations while accounting for bat phylogenetic relatedness. We used the rotl and ape packages in R to extract the bat phylogeny from the Open Tree of Life and calculate branch lengths with Grafen's method (Paradis et al., 2004;Michonneau et al., 2016). We first fit a null PGLS model (intercept only) using the nlme package to estimate phylogenetic signal as Pagel's λ (Pagel, 1999 sample sizes (AICc) and assessed fit with an adjusted R 2 (Burnham and Anderson, 2002). All PGLS models included weighting by sampling variance to account for variable sample sizes per species (Pennell, 2015).
We used generalized linear models (GLMs) to determine which individual-and habitat-level factors influence hair cortisol for each bat species. We first evaluated the relationship between body mass and hair cortisol separately for each species. Next, we ran species-specific GLMs including sex, reproductive stage (by sex) and site disturbance as predictors. Not all covariates were tested for all species due to sample size restrictions. Total sample size and balanced sample sizes among levels were considered to select the number of covariates to include in the model for each species. We included disturbance in GLMs only for species present in more than one site (Pteronotus mesoamericanus, P. mexicanus, Mactotus waterhousii, T. brasiliensis, Glossophaga soricina, D. rotundus) since disturbance was treated as constant within sites. The only genus sampled in both Belize and Mexico was Pteronotus. The two species P. mesoamericanus (Belize) and P. mexicanus (Mexico) represent lineages considered conspecific until a few years ago, but are now thought to represent distinct species that diverged very recently based on molecular and morphometric evidence (Pavan and Marroig, 2016). Because their phenotypes and ecology are still very similar, we treated these as conspecific to test if there were differences in hair cortisol between representatives from the two regions (Mexico and Belize). Tukey post hoc tests were conducted for significant covariates. We compared effect sizes across bat species by evaluating the degree of overlap in 95% confidence interval for each GLM coefficient. All analyses used the natural logarithm of hair cortisol as the response variable and assumed Gaussian errors. We confirmed that all models fulfilled assumptions of normality, homoscedasticity and non-multicollinearity (variance inflation factors < 3). We report data as mean ± SD, unless otherwise noted.

Ecological and evolutionary predictors of hair cortisol
We analysed 259 hair samples from 18 different bat species representing 5 families in Belize and Mexico (Table 1). Hair cortisol concentration across species varied by four orders of magnitude, ranging from 36.6 ± 40.5 ng/g in Eptesicus furinalis to 24 614 ± 14 780 ng/g in L. yerbabuenae (Table 1). Even though mean hair cortisol apparently differed among families (F 4,252 = 18.89; P < 0.01; R 2 = 0.23; Fig. 2), this effect did not hold after accounting for phylogenetic relatedness (F 4,13 = 1.84; P = 0.18; R 2 = 0.16). Accordingly, we did not find strong phylogenetic signal in cortisol (Pagel λ = 0). Ecological and life history traits were instead better predictors of species-level cortisol levels (Table 2). Mean hair cortisol was best predicted by a model including both dietary guild and fecundity; however, only fecundity had a significant effect (F 2,15 = 5.51; P = 0.01; R 2 = 0.34; Table 2). Annual fecundity explained 24% of the variance in Neotropical bat mean hair cortisol. Species reported to have more than one pup per year had significantly lower cortisol than bats having only one pup per year (F 1,16 = 6.22; P = 0.02; Fig. 3A). While phytophagous bats seem to have higher levels of cortisol in hair than animalivorous bats, this difference is not significant when considering the phylogenetic relatedness (Fig. 3B). Other ecological traits including roost durability, foraging style, and colony size were uninformative ( Table 2). As cortisol levels varied by sex for some species (see below), we reran our model comparison after calculating species-level means for males and females separately. However, the above aggregate species-level results held when analysing the sexes separately (Supplementary Tables 2-3).
The lesser long-nosed bat (L. yerbabuenae) showed particularly high hair cortisol (24 614 ± 14 780 ng/g). Because the high values of this species could bias inter-species comparisons, we assessed the sensitivity of our top models by excluding L. yerbabuenae. In Fig. 3C, we show the coefficients from the top PGLS models with and without this species. In both cases, fecundity was the best species-level predictor of hair cortisol regardless of including L. yerbabuenae.

Discussion
Analysis of hair cortisol has become a popular method to study long-term stress in wild animals, offering several practical advantages (e.g. minimally invasive collection, easy sample storage and transport). An accurate interpretation of cortisol levels attributed to stress, however, requires a good understanding of the intrinsic and extrinsic drivers of baseline variation. Factors influencing hair cortisol in bats must be identified before hair cortisol can be used as a conservation tool to assess effects of environmental conditions on bat population health. In this study, we present the first quantification of hair cortisol in bats and describe relationships between hair cortisol levels and both intrinsic and ecological traits. Cortisol in blood, faeces and hair are known to be highly correlated in various mammals (e.g. chimpanzees, chipmunks and mice; Kalliokoski et al., 2019). Therefore, although the concentration values are not directly comparable across different matrices, the effects of covariates can still be compared with our results.
One of the advantages of the use of hair as a physiological biomarker of stress is the longer window of time that it provides compared to other samples like blood and faeces. This feature, however, is only useful if the rate of hair growth and moulting patterns are known for the species of concern. Moulting cycles in bats are understudied, especially in tropical bats (Fraser et al., 2013). Therefore, we cannot be certain of the precise temporal window that this biomarker could represent for our study species. While moult patterns for most of our tropical study species are unknown, we know that many temperate species have a one annual summer/fall moult cycle (Fraser et al., 2013). If we assume our study species have a similar cycle, then the hair cortisol levels detected in our study would be reflecting primarily the circulating cortisol levels during one to two months of rapid hair growth between July and September. Moult in tropical bats, however, might not be confined to a distinct or single annual event, given the relative stable environmental conditions the tropics offer. Moult patterns in tropical bats therefore should be investigated further to allow inferences about the window of time over which cortisol levels in hair integrate and reflect circulating levels. Other aspects associated with differences in the timing of new fur growth, and therefore cortisol deposition in hair, could be related to local environmental variables (e.g. precipitation, temperature and food seasonality; Tiunov and Makarikova, 2007), as well as species-specific life history traits (e.g. reproduction and migration; Heydon et al., 1995).
Overall, we found particularly high hair cortisol in most bat species compared to levels reported in hair for other small mammals (e.g. chipmunk [66-110 g] = 40.27-260.22 ng/g of hair; Mastromonaco et al., 2014). Previous studies in bats, examining plasma and faeces, have also reported higher cortisol levels relative to similar samples from other mammal species (Lewanzik et al., 2012;Kelm et al., 2016;Hald, 2019). Some of the exceptional life history traits of bats, such as long lifespan and low fecundity, could explain why bats exhibit higher levels of GCs compared to other mammals (Austad and Fischer, 1991). According to life history theory, long-lived species with low reproduction rates are expected to prioritize their adult survival (i.e. future offspring) over current reproduction (Stearns, 1992), which could in turn favour higher investment in self-maintenance that might be facilitated by high baseline levels of GCs (Ricklefs and Wikelski, 2002).

Ecological factors among species
Despite the similarity in the HPA hormonal pathways across vertebrates, baseline and stress-induced GC levels are context and species specific (Romero, 2004;French et al., 2008;Crespi et al., 2013;Kalliokoski et al., 2019). In light of this, it is not surprising that hair cortisol levels in bats showed broad interspecific variation. The differences found could not be explained solely by taxonomic family or phylogenetic relatedness (λ = 0), which suggests that other environmental and ecological factors are influencing hair cortisol in Neotropical bats.
Among all the ecological traits evaluated, annual fecundity was the best predictor of hair cortisol. Species with lower fecundity showed higher concentrations of cortisol in hair. This relationship can be supported from a physiological perspective, considering that GCs play an important role modulating the production of reproductive hormones upstream in the hypothalamic-pituitary-gonadal axis, both under homeostatic and challenging conditions. This interaction, however, has only been clearly demonstrated within populations (e.g. individuals with high levels of cortisol reduce reproduction, reallocating energy current needs). It is unknown if interspecific variation in GCs is driven by the same mechanisms. Looking to other taxa, studies on birds have found similar results to those here in Neotropical bats, where avian species with low clutch size and few breeding events showed higher circulating GCs (Bókony et al., 2009;Ouyang et al., 2011). This pattern is supported by predictions from life history theory; for species with lower fecundity, the value of each offspring is higher than in species with relatively high fecundity (Lendvai et al., 2007). Parents of more valuable broods would be predicted to be more 'willing' to invest in offspring survival, which might be facilitated by high baseline GC levels (Bókony et al., 2009). Confirming if this mechanism explains the GCs levels found in our bat species would require more information about reproductive strategies (e.g. monoestry and polyoestry) and seasonality (Wingfield and Sapolsky, 2003), which is limited for most Neotropical species. A more practical rationale for our results could be related to differences in hair cortisol deposition among species with different reproductive cycles. For many species, pelage moulting occurs after breeding (Constantine, 1957;Dwyer, 1963;Cryan et al., 2004;Measor et al., 2017), given that both reproduction and hair growth are energetically demanding processes (Ling, 1970). Based on this pattern, we could expect that cortisol deposition in hair may reflect differences in breeding events among bat species.
Although diet explained additional variation in Neotropical bat hair cortisol, this variable was uninformative when considering phylogenetic relationships; hair cortisol did not vary significantly between our simplified dietary guilds. Further studies using a more accurate classification of diet (e.g. using stable isotopes or metabarcoding; Oelbaum et al., 2019;Ingala et al., 2021) could give more conclusive insights into the links between feeding strategies and cortisol levels in bats.
GCs play a key role in metabolic function, facilitating fuel mobilization (e.g. glucose, fatty acids) under normal and challenging conditions (Kuo et al., 2015). A positive relationship between resting metabolic rate (RMR) and plasma cortisol levels has been reported for various mammalian species (including four species of bats), and this relationship has been suggested as a general pattern for mammals (Haase et al., 2016). Due to the limited data on RMR for our study species, we used BMR as an indicator of energy expenditure. Different to what we expected, BMR was not an informative factor for cortisol variation in hair among the bats in our study. The positive relationship between cortisol levels and metabolic rate previously found in plasma might be obscured in studies of hair cortisol like ours, due to confounding factors such as moulting cycles and cortisol deposition rate. In addition, obtaining accurate BMRs in wildlife species (particularly freeranging animals) is challenging, which raises questions about the quality of BMR data, especially in comparative studies (Genoud et al., 2018). For future studies, a more realistic and informative indicator of energy turnover in free-ranging animals is the Daily Energy Expenditure (Speakman, 1997), which integrates the energy allocated in different activities such as foraging, commuting and thermoregulation (Butler et al., 2004).
The relationship between body condition and GC release has been widely evaluated, because weight loss is one of the early responses to long-term stress in many species ( Angelier et al., 2009;Dickens and Romero, 2013). However, the direction of the effect of body condition on cortisol is context and species dependent (Crespi et al., 2013). We used body mass as an indicator of body condition because it is a more informative metric than other indices in bats (McGuire et al., 2018). We found different directions of the effect of body mass on hair cortisol. For two of the studied species (M. nigricans and P. mesoamericanus), heavier individuals showed higher concentrations of hair cortisol. In contrast, P. mexicanus showed a negative relationship between body mass and cortisol. These divergent results suggest that the relationship between cortisol and body condition of bats is not generally predictable and might be species specific. Hair cortisol only reflects the time window of fur growth. Therefore, it is possible that body condition might have changed since the individual's last moult.
One species that stood out for its particularly high levels of cortisol was L. yerbabuenae. This species is highly mobile and migratory (Horner et al., 1998;Buecher and Sidner, 2013;Medellin et al., 2018). Migration itself was not considered in our analyses, because the degree to which bats may migrate seasonally is unclear for many of the species in our sample. Migratory behaviour, however, could explain such high cortisol concentrations in L. yerbabuenae. La Fábrica caves in Colima, one of our field sites, is known to be one of the starting points of the annual migration of L. yerbabuenae (Medellin et al., 2018). The role of GCs during migration has been widely studied in birds, fish and some large mammals, but not in bats (Holberton, 1999;Romero, 2002;Wada, 2008). We hypothesize that premigratory fattening could explain the high hair cortisol levels observed in L. yerbabuenae, and we encourage future studies to address this question.
Consistent with other studies, we found differences in hair cortisol levels between sexes, albeit for only 4 of our 18 studied species: D. rotundus, M. nigricans, P. mexicanus and P. mesoamericanus. For these species, females showed higher cortisol than males, a trend that appears to hold for many mammalian species (Bechshøft et al., 2011;Hau et al., 2016;Rakotoniaina et al., 2017;Dettmer et al., 2018). Higher levels of GCs in females can be attributed to sex differences in the HPA axis activity, which are mainly mediated by gonadal steroid hormones (i.e. androgens and estrogens). For example, estradiol, which is more abundant in females than males, can enhance cortisol release, while androgens tend to reduce its production (Handa et al., 1994). Other mechanisms underlying sex differences in HPA axis regulation and stressrelated behaviours in mammals are reviewed by Zuloaga et al. (2020). Females have also shown differences in HPA axis activity depending on their life history stage, with GCs being higher during the late stages of pregnancy (Reeder et al., 2004). Studies in a fruit-eating bat (Artibeus jamaicensis) and little brown myotis (Myotis lucifugus) have reported higher levels of plasma GCs in pregnant females (Reeder and Kramer, 2005;Klose et al., 2006). Contrary to those findings, we did not find reproductive state to influence hair cortisol in our female-only model (i.e. for P. mesoamericanus in Belize). However, it may have been difficult to detect an effect, given the low number of pregnant females in our sample (n = 5, 24%) and the fact that moulting might not occur in conjunction with mating. Although cortisol has been proposed as the primary GC in bats (reviewed by Kwiecinski and Damassa, 2000), corticosterone is also detectable in circulation and has been identified to play an important role in reproduction . Therefore, some effects of ecological traits could go unnoticed and the complexity of the stress response in bats could be oversimplified by quantifying only cortisol.
Bats have been proposed as good indicators of habitat quality due to their ecological diversity, wide distribution and potential sensitivity to disturbance (Jones et al., 2009a;Cunto and Bernard, 2012;Stahlschmidt and Brühl, 2012). However, a clear correlation between environmental disturbance and cortisol levels, in faeces and blood, has not been reported in bats (Wada et al., 2010;Allen et al., 2011;Kelm et al., 2016). Cortisol in hair could reflect better the effects of chronic stressors such as human settlements, and it is not sensitive to capture stress. We compared hair cortisol in three of our study species found in sites with varying fragmentation and agricultural activities in Mexico. We found an effect of disturbance in only one of these species, P. mexicanus, for which bats roosting in Don Pancho Cave island, a site with moderate disturbance, showed the highest concentrations of cortisol (Fig. 1). We speculate that the high levels found in this population could reflect differences in the cave microhabitat compared to the other caves in our Mexican sample. Don Pancho Cave is a narrow crevice estimated to have a higher colony size (100 000 individuals from 6 species; Téllez et al., 2018) than the other sampled caves El Salitre (∼10 000 individuals from 10 species; Torres-Flores et al., 2012) and La Fábrica (>5000 individuals from four species). The high density of bats in the Don Pancho cave may increase agonistic social interactions (Creel et al., 2013) and parasite transmission (Langwig et al., 2012;Postawa and Szubert-Kruszyńska, 2014), factors that have been shown to increase cortisol levels in other mammals.
Physiological responses to chronic stress in wildlife are difficult to unravel and predict unless multiple responses at different levels of biological organization are evaluated simultaneously (Dickens and Romero, 2013). Hair cortisol offers great potential as a tool to monitor health in wild populations, particularly those already identified at risk (Kalliokoski et al., 2019). For instance, chronically elevated cortisol levels have been linked to greater susceptibility to infection and disease severity (Davy et al., 2017). Periodic surveys of hair cortisol could therefore help identify periods when bats might be more vulnerable to infection (e.g. white nose syndrome). Further, such surveys might also inform when individuals are more likely to shed zoonotic pathogens (e.g. henipaviruses and filoviruses; Plowright et al., 2008;Davy et al., 2017;McMichael et al., 2017;Kessler et al., 2018).

Conclusions
The current study reports cortisol levels in hair of 18 Neotropical bat species from two countries and serves as a reference for future research using this method in wild bat populations. We found that fecundity and potentially diet are important ecological traits explaining interspecific variation in bat hair cortisol. Within species, female bats exhibited higher cortisol than males and the effect of body mass varied among species. Other factors that may be important at the individual level, such as parasite load and colony size, should be considered in future studies to have a more complete understanding of sources of variation on baseline GC levels within species. Importantly, studies looking at hair growth rate and moulting cycles in Neotropical bat species are imperative to give an accurate interpretation of hair cortisol as a biomarker of stress response. Applied properly, hair cortisol quantification is a powerful minimally invasive technique with multiple potential applications in bat ecology, physiology and conservation. Our findings and ongoing work will help to validate and apply hair cortisol as a monitoring tool in wild bat populations.