Hair cortisol concentration reflects the life cycle and management of grey wolves across four European populations

The grey wolf (Canis lupus) persists in a variety of human-dominated landscapes and is subjected to various legal management regimes throughout Europe. Our aim was to assess the effects of intrinsic and methodological determinants on the hair cortisol concentration (HCC) of wolves from four European populations under different legal management. We determined HCC by an enzyme-linked immune assay in 259 hair samples of 133 wolves from the Iberian, Alpine, Dinaric-Balkan, and Scandinavian populations. The HCC showed significant differences between body regions. Mean HCC in lumbar guard hair was 11.6 ± 9.7 pg/mg (range 1.6–108.8 pg/mg). Wolves from the Dinaric-Balkan and Scandinavian populations showed significantly higher HCC than Iberian wolves, suggesting that harvest policies could reflected in the level of chronic stress. A significant negative relationship with body size was found. The seasonal, sex and age patterns are consistent with other studies, supporting HCC as a biomarker of chronic stress in wolves for a retrospective time frame of several weeks. Our results highlight the need for standardization of sampling and analytical techniques to ensure the value of HCC in informing management at a continental scale.

, thus addressing the pressing issues of animal welfare in wildlife research 12 . Additionally, hair collection is simple and inexpensive when performed on specimens in archives and on recently dead animals 13,14 . This facilitates achieving larger sample sizes and wider spatial-temporal coverage than with invasive sampling.
Cortisol in blood or feces reflects the activation of the HPA axis from minutes to a few days before sample collection 6,9 . In contrast, cortisol is hypothesized to be integrated into the hair over a longer time frame, providing a picture of the activation of the HPA axis over a time frame of weeks to months 4,6,11 . Hence, this method offers an opportunity to study the physiological responses of wildlife to natural processes and potential long-term stressors such as social interactions, infectious diseases, or human perturbation 15 .
The grey wolf (Canis lupus; hereafter, wolf) inhabits different human-dominated landscapes across its range 16,17 . While the historical range of the wolf covered most of the northern hemisphere, by the end of the nineteenth century it was exterminated from most of western, central, and northern Europe 16,18 . In the last decades, the wolf has recolonized large parts of its former European range, including human-dominated landscapes, with an average human density of 36.7 inhabitants/km 216 . Wolf management regimes vary across Europe, with the species being fully protected in parts of the Iberian and Alpine population ranges, while hunted or subjected to legal population control in much of the Scandinavian and Dinaric-Balkan ranges 19,20 . This makes European wolf populations a relevant model for studying conservation physiology.
The aim of this study was to assess the effect of intrinsic determinants (sex, age, body condition, body structural size, month, and cause of death/capture, and wolf population) on HCC, as determined by enzyme-linked immune assay (ELISA), of European wolves from the Iberian, Alpine, Dinaric-Balkan, and Scandinavian populations while controlling for the confounding effect of variables related to the methods employed (weight and length of hair used for cortisol extraction, sample storage time, body region and methanol evaporation protocol).
The tail and ventral thorax regions, but not dorsal thorax, showed significantly lower HCC than the lumbar region (Table 1 and Fig. 2A,B), after the exclusion of 4 outliers. All subsequent analyses were performed for lumbar hair samples only.
The mean HCC of lumbar guard hair, excluding 2 outliers (see "Methods"), was 11.64 ± 9.65 pg/mg, ranging 1.58-108.81 pg/mg. The most supported model (Model 1, Table S4), with a weight 0.122, included the fixed effects sex, age, population, body size (total length), and cause of death/capture ( Table 2). The cumulative weight of the 4 models with ΔAICc < 2 was 0.320. The sum of the variable weights in the 4 models with ΔAICc < 2 (Table S4) was 1 for age, population, body size, and cause death, 0.85 for sex, 0.32 for evaporation protocol, and 0.15 for hair length. Model 1 conditional R 2 = 0.232, and the null model's ΔAICc = 24.5, supporting the ability of the former to explain the variation in HCC.
The most supported model revealed significantly higher HCC in the Dinaric-Balkans and Scandinavian populations (Fig. 3A). Hair cortisol concentration showed a significant negative relationship with total length (Fig. 3B), and higher values in the winter and lowest in May and November (Fig. 3C).

Discussion
The mean HCC in lumbar guard hairs of wolves from the four European populations studied (11.6 ± 9.7 pg/ mg) was lower than that reported for wolves from the Canadian tundra-taiga (16.5 pg/mg) and northern boreal forest (13.5 pg/mg) 15 . It was suggested that the degree of legal lethal removal by humans (licensed hunting or population control) was reflected on the HCC of Canadian wolf populations, possibly mediated by the increased social instability of wolf populations subjected to intense harvest 15 . Therefore, the mean HCC in our sample, being lower than in the less heavily harvested Canadian boreal forest population, could imply a lower intensity of legal lethal removal of the European wolf populations studied. Bryan et al. 15 do not quantify the proportion of their study populations that is harvested but mention it can be up to 50% annually. In the same Canadian province of Alberta, 34% of the wolf population was estimated to be legally harvested every year 21 . In Europe, the disappearance rate of Scandinavian wolves was estimated at 23%, with legal harvest amounting to 5-15% annually 22 . In the Alpine wolf population, annual mortality was estimated at 10-25% 23,24 . The apparently lower human-related harvest of wolves in Europe, compared to Canada, could explain the lower HCC in our study.
Wolves were fully protected in Portugal and Italy but subjected to varying degrees of legal harvest in Spain, Slovenia, Sweden as well as in other countries of the Alpine, Dinaric-Balkan, and Scandinavian population ranges that were not sampled in this study 19 . The populations showing significantly higher predicted HCC (Dinaric-Balkan and Scandinavian) were subjected to legal harvest, while the partially protected Iberian and Alpine populations did not significantly differ between them. Even fully protected wolf populations in Europe experience high levels of poaching 22,25,26 , which could dampen the differences in the legal wolf removal and in  www.nature.com/scientificreports/ HCC across the four European populations studied. Furthermore, extrinsic factors not addressed in this study, such as wolf density, prey availability and non-lethal human disturbance can also influence the HCC 11,13 . The potential of HCC as biomarker of the intensity of legal and illegal human-related harvest of wolves across Europe needs to be further established.
On the other hand, the differences in HCC between European and Canadian wolves could also be caused by slightly different methodologies employed by Bryan et al. 15 , such as a non-standardized body region for the collection of hair samples, hair grinding by ball mill or a different commercial ELISA. The grinding method was shown to have a significant effect on the amount of cortisol extracted from hair, possibly mediated by matrix and surface area mechanisms 27,28 . In our study, it was decided to use a closed-tube method, with the advantage of preventing loss of hair powder and sample contamination that can occur when using a ball mill 27 . Moreover, large differences in HCC according to the ELISA protocol utilized have also been described 29,30 and could contribute to explain the differences in measured HCC to Canadian wolves 15 . Despite the unplanned change in the evaporation protocol (nitrogen gas stream vs suction hood evaporation) accounting for a non-significant effect on HCC, this was the most relevant methodological variable (weight 0.32). Full standardization of the analytical protocol is strongly recommended.
Our results show difference in HCC across the body regions from which hair was collected. Results from other studies are inconsistent, as differences in HCC according to body region have been reported in some species but not in others 11 . The main mechanism hypothesized for the incorporation of glucocorticoids into growing hair is via the blood vessels and eccrine and sebaceous glands that surround the hair follicle 4 . Guard hair from different body regions may present disparate periods and rates of growth and therefore cortisol will not be incorporated equally throughout the body 13 . Wolves seem to have one annual molt, beginning in late spring when the old coat is shed, and the new summer coat grows continuously through fall and winter 18 . Relating hair samples to a specific period is more challenging in species with poorly known molt patterns 15 . Our results support that hair should be collected from standardized body regions to improve inference from HCC.
The small number of wolves in this study showing viral, bacterial, or ectoparasitic infections or neoplasia (n = 5, including 1 outlier excluded from Model 1-see "Methods"), some of which caused their deaths, tended to present extreme HCC compared to apparently healthy wolves. Chronic diseases have been shown to cause major increases in HCC in many species 11 and higher levels were found in chronic compared to acute diseases in cattle 31 . The highest HCC in our sample were associated with sarcoptic mange, which can course for many weeks to months in individual wolves 32 . In contrast, death by events likely lasting minutes to a few days or live www.nature.com/scientificreports/ trapping did not translate into significant alterations to HCC. Together, these results support the hypothesis that HCC reflects cortisol levels in the organism in a retrospective time frame of several weeks 4,11 . No evidence was found for lower HCC in samples stored for up to 4 years, suggesting that archived specimens could be used to provide insights into the long-term dynamics of chronic stress in wolf populations. Nevertheless, this study was not designed to evaluate the effect of storage time on HCC, which would be better accomplished  Other studies report inconsistent results on this potentially confounding variable, from a significant negative effect of storage time in the Egyptian mongoose Herpestes ichneumon 33 to higher HCC in polar bear Ursus maritimus hair stored for > 100 years than in recent samples 14 .
Males showed non-significantly higher HCC than females. Other studies consistently report non-significant differences between sexes, with higher HCC in females 15,35 or males 34 . Younger age classes, particularly subadults, showed a tendency for lower HCC than the adults. The same age pattern was reported for HCC and fecal cortisol metabolites in North American and European wolves 34,35 .
The monthly variation in HCC predicted by the model is consistent with the annual life cycle of wolves and with the results of fecal glucocorticoid metabolites [34][35][36] . The highest HCC was found in the winter, at the start of the mating season 37 and could be a consequence of increased social instability related to sexual and territorial behavior 11,34 . Wolves are strict seasonal monoestrous breeders 18 . It has been shown in other species that cortisol can play a role in reproduction by inhibiting the action of oestradiol and luteinizing hormone and, consequently, inhibiting female receptivity to mating and ovulation 38 . The surge in cortisol during the mating season could thus be a contributing mechanism to regulate reproduction within wolf packs, where usually only one pair breed, even when other potential breeders are present 18 . The lowest HCC was found in May, coinciding with parturition, during which intraspecific competition and social instability are lessened 37 . Between May and October the HCC progressively increases which could be related to the growing effort to supply pups. The peak in HCC in October coincides with juveniles leaving their homesites and start travelling with the pack 36,39 .
Our results support the suitability of hair cortisol concentration as biomarker of chronic stress in wolves, suggesting that it allows the assessment of cortisol levels in a retrospective time frame of weeks to months, as postulated for this biological matrix 4,6 . This is supported by the observation that causes of death with an expected duration of days to weeks do not translate in higher HCC, while the small number of wolves with chronic infections, particularly sarcoptic mange, showed the highest values. The reliability of this approach is further supported as the observed annual cycle of HCC is consistent with the levels of chronic stress inferred from the wolves' life cycle and other studies measuring fecal cortisol metabolites [34][35][36] . The sex and age patterns are also similar to those reported in other studies with wolves, measuring hair cortisol or fecal cortisol metabolites 15,[34][35][36] . This study supports that sampling (e.g. body region, cause of death) and analytical (e.g. methanol evaporation protocol and length of hair) techniques should be standardized or accounted for in the statistical analysis to improve inference from HCC data.

Collection of wolf hair samples. Hair samples were collected by researchers from opportunistically
found-dead wolves upon standard necropsy (all the Alpine and part of the Iberian samples) or in the field (all the Dinaric-Balkan and most of the Iberian samples), or from legally harvested wolves (only in the Scandinavian population). At the time of sample collection, wolves were legally harvested in Sweden, Slovenia, and Spain, and under total protection in Portugal and Italy. Hair samples were collected from four body regions, when possible: lumbar (n = 133), dorsal cervical (n = 66), tail (n = 33) and ventral thorax (n = 27) (Tables S1 and S2). The hair was cut as close as possible to the skin with scissors to avoid collecting hair follicles, but in some samples, hairs were pulled from the carcass. Samples were stored at room temperature in paper envelopes. Age, sex, date, and cause of death/capture, geographical location, body mass, and total length were obtained for most of the wolves.
Age was estimated by the dental eruption and wear or cementum age analysis and classified as 'juveniles' (< 1 year old), 'subadults' (1-2 years), 'adults' (> 2 years) 40 , or 'unknown' . Sex was assessed by inspection of genitalia. Causes of death were classified as 'acute' , likely lasting minutes to hours (vehicle accident and legal or illegal shooting); 'subacute' , likely lasting hours to days (drowning, poisoning, trapping and intraspecific aggression); 'chronic' , likely lasting several weeks (infectious diseases-canine distemper, canine parvovirosis, leptospirosis; sarcoptic mange; or neoplastic diseases) or 'unknown' . Total length was obtained by measuring with metric tape (1 mm precision) the distance from snout to the distal end of the last tail vertebrae. The body mass was measured with 100 g precision with scales.
The detailed protocol for the handling of wolves live trapped in the scope of ecological and conservation studies (n = 7, all from the Iberian population) has been previously described 5 . Traps were monitored twice every day, in the early morning and late afternoon, hence the duration of restraint after capture was unknown Cortisol extraction. The protocol for the extraction of cortisol from the hair was adapted from previously described procedures 15,27 . Forty mg of guard hairs were separated from the undercoat and placed in 15 ml falcon tubes. Hair follicles were cut whenever found in the sample. For each sample, the length of three intact hairs was recorded. The samples were washed twice with 40 µl of distilled water/mg hair and three times with the same amount of isopropanol. In each washing step, the samples and washing solution were vortexed, the supernatant discarded, and the hair dried using clean paper towels. After the final wash, samples were dried overnight at www.nature.com/scientificreports/ room temperature and 30 mg of hair cut into a 2 ml polypropylene screw cap plastic tube with five 4 mm steel beads added to each tube. The hair was ground to a fine powder in a FastPrep sample homogenizer (MP Biomedicals, USA) for four times 1 min at 6.0 m/s. 50 µl methanol/mg hair were added to each sample and sonicated for 30 min at 50 Hz at 50 °C. The samples were incubated for 18 h at 50 °C in an orbital shaker at 160 rpm, centrifuged for 15 min at 14,000g at 20 °C, and 1000 µl of supernatant was collected to a screw cap glass chromatography vial and dried at room temperature in a gentle stream of nitrogen gas. Due to restrictions on laboratory use during the SARS-Cov-2 pandemic, some batches of samples were instead evaporated overnight on a suction hood. This unexpected change in the methanol evaporation protocol was recorded and accounted for in the statistical analysis.
Cortisol quantification. A commercial competitive ELISA kit (Cortisol free in Saliva ELISA, Demeditec, Germany) was used to quantify the concentration of cortisol, following the manufacturer's instructions. The kit plate wells are provided coated with polyclonal rabbit antibody against cortisol, and cortisol-horseradish peroxidase was used as conjugate. According to the manufacturer, the cross-reactivity of the test to selected steroids is low (Table S3), the intra-assay variation is 3.8-5.8% and the inter-assay variation is 6.2-6.4%. Samples, standards, and controls were tested in duplicate.
The 4-parameter standard curve was calculated from the log-transformed cortisol concentration of the standard solutions and their measured OD 450 44 . Standard curves were estimated using the software GraphPad Prism 6.04 (GraphPad Software, La Jolla, California USA), and yielded an average R 2 adjusted = 0.991 (range 0.968-0.999). The cortisol concentration of the reconstituted samples was estimated from the standard curve and converted to cortisol concentration as picograms (pg) of cortisol/mg of guard hair.
Intra and inter-assay coefficients of variation were estimated for six ELISA assays of 37-40 samples each. The low and high controls included in the kit were used to estimate the inter-assay coefficient of variation and the duplicate runs of each sample were used to estimate the intra-assay coefficient of variation. Linearity was assessed by two-fold dilutions (1:1, 1:2, 1:4 and 1:8) of 4 extracted samples, comparing the expected and observed concentrations. Recovery was assessed by spiking 6 ground hair samples with known concentrations of cortisol (50, 25, 12.5, 6.25 pg/mg, and no spiking), comparing the expected and observed concentrations.

Determinants of hair cortisol concentration. The potential determinants of HCC investigated
included wolf intrinsic variables: sex, age, body condition, body structural size, month of death/capture, and wolf population. The scaled mass index was selected as a measure of body condition 45 and estimated from the log-transformed body weight (g) and total length (mm). Log-transformed total length was used as an indicator of body structural size 46 . Samples were assigned to the Iberian, Alpine, Dinaric-Balkan, or Scandinavian wolf populations 16 from the geographical location of the death or live-trapping sites (Fig. 1).
The relationship between HCC and additional variables related to the sampling procedure or to the work conducted in the laboratory (length of hair used for cortisol extraction, sample storage time, body region, cause of death/capture, and methanol evaporation protocol), herein referred to as methodological variables, was also investigated as potential confounding variables. Sample storage time was the period in months between death/ capture and cortisol extraction. In those samples for which only the year of death was available, 30 June was assigned as the date of death, solely to estimate storage time. All continuous variables were standardized to their z-scores.
Statistical analysis. First, the effect of body region was investigated by a linear mixed model with HCC as the dependent variable, and the independent variables body region, as a categorical fixed effect, and individual wolf, as a random effect. The lumbar region was set as the reference class as it was the most represented in our sample (Table S1). Data from 27 wolves for which samples were available from all 4 body regions were used in this analysis. Four outliers in the dataset violated the assumption of normality in the residuals of the model comparing HCC across body regions (Fig. S1A) and were excluded from this model's dataset (Fig. S1B).
Second, the effect of intrinsic and methodological variables on HCC from the lumbar body region was investigated by another linear mixed model with sex, age, body condition, body structural size (standardized log-transformed total length), cause of death/capture, wolf population, hair length, sample storage time, and methanol evaporation protocol as fixed effect independent variables. The month of death/capture was included as a random effect. Reference classes for the categorical variables were set as female, adult, acute death, Iberian population, and methanol evaporation by nitrogen gas stream. Two outliers in the dataset violated the assumption of normality in the residuals of the model (Full model, Table S4) and were excluded from this analysis (Fig. S1C,D).
The goal of this analysis was to assess the relationship between HCC and wolf intrinsic variables, controlling for the potential confounding effect of the methodological variables. Starting from the full model (Table S4), models including all possible combinations of variables were ranked by their AICc using the package "MuMIn" 47 in R 3.6.1 48 . The most supported model was selected for inference and models with ΔAICc < 2 are reported in Table S4.
The HCC predicted by the most supported model was estimated using the package "merTools" 49 . The function "predictInterval" was used, which fits multivariate normal distributions to the random and fixed effects. www.nature.com/scientificreports/ 1000 values were sampled from these distributions for each category of the random (month) and fixed effects, capturing the full uncertainty in predictions as 95% confidence intervals. Linear mixed models were fitted using the package "lme4" 50 in R 3.6.1 48 . The correlation between fixed effects was estimated with a threshold for acceptance of 0.700. The conditional R 2 of the model was estimated according to Nakagawa and Schielzeth 51 implemented in the package "MuMIn" 47 . The assumption of normality of the model residuals was checked by inspection of quantile-quantile plots. Graphics were produced using the R package "ggplot2" 52 .

Data availability
The datasets generated during the current study are available from the corresponding authors on reasonable request.

Scientific Reports
| (2022) 12:5697 | https://doi.org/10.1038/s41598-022-09711-x www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.