The farm visit was planned so that the undisturbed animals could be observed for 20 min around morning feeding time as this is a time point of high activity (Bokkers and Koene, 2001). If multiple pens were included on a farm, the observer rotated between the pens every five min until each pen had been observed for a total of 20 min. Two out of 42 pens contained more than 20 animals and were too large to observe the whole pen from one location. In these cases, the pen was observed for 2 × 20 min from two opposite locations to account for the possibility of missing behaviors performed by animals standing far away from the assessor. The protocol, which consisted of clinical scoring, undisturbed behavioral observations, and a behavioral test, was developed based on previous publications on the assessment of calf welfare (Lidfors, 1993; Jensen and Kyhn, 2000; Bokkers and Koene, 2001; Whay et al., 2003; Windschnurer et al., 2008; Bokkers et al., 2009; Leach et al., 2009; Welfare Quality®, 2009). The full protocol is presented in Table 1. The measures consisted of a combination of long-term indicators, such as body condition score and coat condition, and measures reflect a shorter time period, such as ruminal fill. When applicable for the behavior being observed, such as mounting and licking, both performer and receiver were noted. After the undisturbed behavioral observations, fear of unknown humans was tested using a distance avoidance test. Finally, calves were restrained and individually assessed e.g., for nasal discharge, tear staining, and sampled for VTEC and hair. At this time, the reactivity of the calf to restraint was also scored. All animal-based assessments were carried out by the first author.

Hair samples were collected from the tail of each calf (just above the tail switch) with scissors. Hair was cut as close to the skin as possible. Samples were wrapped in aluminum foil and stored dark at room temperature for 3–7 months before analysis. From 318 sampled animals, 196 samples including both individuals colonized by VTEC O157:H7 and randomly selected negative controls (1–3 controls per colonized individual) were selected for analysis of hair cortisol. The protocol for preparing and analyzing hair was based on previous studies on cattle and other species (Meyer et al., 2014; Tallo-Parra et al., 2015; Roth et al., 2016) but adapted and modified by authors. Prior to analysis, hair from the centimeter closest to the cut end of the hair (closest to the skin) was cut with scissors and placed in a 15-ml tube. Samples were washed by adding 4 ml of isopropanol, and samples were vortexed for 4 min before the isopropanol was decanted. The procedure was repeated three times. Samples were left to dry at room temperature under a fume hood for 2-3 days until all the isopropanol had evaporated. For samples where fecal contamination remained after the washing procedure, the dirtiness of the hair was scored from 1 to 4 (1= no dirt, shiny, 2 = <5 dirt particles, shiny, 3 = > 5 dirt particles, shiny, 4 = > 5 particles of dirt, dry, or dull). Approximately 50 mg of each hair sample was added to a 2-ml microtube with a screwcap. If less than 50 mg of washed hair was available, as much as possible was used. Exact duplicates were prepared for 30 samples (12%). In addition, a duplicate sample with a different degree of remaining dirt or of a different color was prepared from 21 of the individual samples. Samples were prepared farm by farm in random order. For exact duplicates, every 10th sample was selected during preparation (if the sample did not have enough hair, the 11th sample was taken). The 21 samples of different colors and dirtiness were selected based on availability during sample preparation.

Three chrome steel balls (3.2 mm in diameter, BioSpec Products, Cat. No. 11079132c) per tube were added, and the tubes were frozen in liquid nitrogen for 2 min before grinding the hair for 3 × 60 s at 6,500 RPM in a bead beater (Precellys Evolution). After grinding, 1.2 ml of methanol was added and samples placed on a rocker. After 18 h, the samples were centrifuged. The first 41 samples were centrifuged for 2 min at 7,000 g, and 0.6 ml of methanol was transferred to new 1.5-ml tubes. As the first samples had occasional hair particles in the final product, the centrifugation protocol was changed and the remaining samples were centrifuged two times After the initial centrifugation, 0.8 ml of methanol was transferred to new 1.5-ml tubes, and the samples were recentrifuged at 10,000 g for 5 min, and 0.6 ml of methanol was transferred to new 1.5-ml tubes. Methanol was evaporated at 38°C under a fume hood after the centrifugation. Phosphate-buffered saline (0.01 M PBS, pH 7.4) was added to the dried samples. The first 92 samples were dissolved in 200 μl PBS, the last 159 samples in 150 μl as the dried sample appeared to dissolve better in this volume. The concentration of cortisol was determined using an ELISA kit designed for salivary cortisol (Salimetrics Europe Ltd, Art 1-3002) according to the manufacturer's instructions (validated by Moya et al., 2013). Each sample was analyzed in duplicate, and high and low controls were run on each assay. Inter- and intraassay coefficients of variation (CV) values were calculated according to the instructions from the manufacturer (Salimetrics, 2018). Final cortisol content (pg/mg) in hair was calculated using the formula presented by Meyer et al. (2014).

Data were entered in Excel (Microsoft Corp., Redmond, WA) and exported to R Statistical Software (R Core Team, 2018) for statistical analysis. Hair cortisol was log10-transformed to achieve normal distribution before calculation of Pearson's correlation coefficient for exact duplicates and duplicates of different dirtiness or colors. Plots of the correlation were generated using the package ggplot (Wickham, 2016). The impact of remaining dirt, age, hair color, month, and methodological changes during the process was analyzed using a generalized additive model (GAM) with a log link function in the package gamm4 (Wood, 2004; Wood and Scheipl, 2017). In this analysis, one sample from each pair of exact duplicates was excluded from the analysis, but all other remaining samples were included. A random effect of pen nested within-herd was included to account for the clustering of the observations. The model was reduced based on the change in R-squared after the removal of one variable at a time, and the final model was evaluated by residual plots.

To explore and visualize the observations of the calves and HCC, a cluster analysis of mixed data was performed using the CluMix package in which distance between individual calves was based on Gower's general similarity coefficient and variables clustered based on similarity using measures of association (Hummel et al., 2017, 2019).

The associations between hair cortisol (log10-transformed) and welfare indicators were also analyzed using regularized regression using package glmnet (Friedman et al., 2010). For animals with multiple samples, only one sample (with the lowest CV value) was included. Associations between welfare parameters and HCC were first explored excluding the samples with fecal contamination (n = 139) and second with all samples, but then with the dirtiness score included in the analysis (n = 196). There were small amounts of missing data in the welfare indicators randomly scattered across observations. The proportion of missing data per indicator was low, for the majority of indicators <1%. The indicator with the highest proportion of missing values was cleanliness body (8 missing) and distance avoidance test (7 missing at random, 23 missing by design). These were followed by cleanliness upper hind legs (6), cleanliness below hocks and coat condition (5), and ruminal fill (4). Another 10 variables were missing 1 observation. Missing values were assumed to be missing at random except for the distance avoidance test where 23 observations were missing due to the pens being too small or crowded for carrying out the test. These were treated as a separate category in the analysis. The remaining missing values were imputed using non-parametric random forest imputation in the package missForest before the elastic net regression (Stekhoven and Bühlmann, 2012; Stekhoven, 2013). Two variables (lameness and conjunctivitis) were excluded due to little variance between the calves, and qualitative indicators with several levels were reduced by grouping if some levels contained very few observations and grouping was biologically plausible (Supplementary Table S2). The counts of hairless patches and wounds or swellings were also condensed into categories (Supplementary Table S2). For many behavioral observations, a high number of calves did not perform the behavior leading to a distribution with the high frequency of zero observations. This was handled by adding a categorical variable describing if the behavior was observed and a quantitative variable with the number of times the behavior was observed (Robertson et al., 1994). The models were tuned and fitted using leave one out cross-validation using the package caret (Kuhn, 2020), and the models that generated the smallest mean squared errors were selected (James et al., 2013). The herd was included as a fixed effect to account for clustering.

The GAM showed that the dirtiness of hair was associated with a significant increase in cortisol extracted and that the increase was larger for the dirtiest samples (Table 2). There was a significant peak in hair cortisol levels in calves around 50 days of age and a tendency for darker hair to be associated with higher HCC (Figure 1). Other parameters (storage time, number of centrifugations, methanol) did not influence the cortisol concentrations (Table 2). A preliminary analysis indicated that lower HCC was seen in the samples where 200 μl of PBS was used to dissolve the cortisol, which supported our suspicion that the dried samples did not dissolve completely. However, in a stratified analysis of samples from farms where both preparation volumes were applied, no difference was observed. It was, therefore, concluded that the effect observed was a result of farm differences and not the amount of PBS. Also, removing the PBS variable did not impact the explanatory abilities of the model, further supporting that the effect was related to farm variations. Thus, this variable was excluded from the final model.

Detailed results of the welfare and behavior assessments (including all samples regardless of contamination level) and a visualization of the cluster analysis can be found in the (Supplementary Tables S1, S2; Supplementary Figure S1). HCC clustered (as expected from the previous analysis) closely to the dirtiness score. Interestingly, hair cortisol also clustered closely to score in the distance avoidance test (which was in the same cluster as the dirtiness score) although the pattern of association between HCC and the score in the avoidance distance test was not distinguishable. Most notably, a group from Herd 1, where the distance avoidance test was missing by design (indicating a small and crowded pen), appeared to have higher HCC levels. In addition, HCC clustered close to age, herd, sex, being fed milk, and stretching (although the latter had relatively few observations). Young animals being fed milk formed a cluster that appeared to have higher HCC values compared with the older animals. Animals with poor body conditions (and poor coat conditions) also showed high similarity and clustered closely together. However, in this group, cortisol levels did not appear to increase but decreased. There was also a cluster where calves with a high frequency of coughing during the behavioral observations and with nasal discharge and coughing observed during the clinical observations appeared to have a slightly lower HCC. There was some clustering associated with the herd, but the degree of similarity within-herd differed between herds.

Regularized regression using clean samples and all samples yielded comparable models with similar R-squared (0.59 and 0.58) and mean squared error (0.42 and 0.46). In the model including all samples, the dirtiness score of hair (Scores 3 and 4) had the largest coefficients (Score 2 was also included but with a smaller coefficient). The coefficients of variables describing the association between HCC and the welfare indicators in the two models are presented in Figure 2. Variables with positive coefficients (i.e., associated with increased hair cortisol) in both models were a mixture of potential indicators for both poor and good welfare (Boissy et al., 2007; Botreau et al., 2007). Being fed milk had large coefficients in both models. The behaviors with the largest positive coefficients, i.e., affecting cortisol the most, present in both models were the qualitative and quantitative variables for being displaced and being rubbed on by others. That both the qualitative and the quantitative indicators for displacement were included in the final models which implies that the effect of being displaced on cortisol concentration increases the more calves are being displaced. Together with poor ruminal fill (model with all samples), these variables are potential indicators for competition and/or low dominance. Moderate tear staining was also included in both models. There were also variables with smaller positive coefficients not usually related to poor welfare (e.g., locomotor play, which is considered an indicator of good welfare, and grooming or comfort behaviors such as self-licking and rubbing). The coefficients of the variables associated with low cortisol levels (i.e., variables with negative coefficients) in both models were related to poor health (nasal discharge, coughing), dirtiness on the upper hind legs (indicating diarrhea), and high reactivity during sampling. Performing stereotypic behaviors, such as cross-sucking and tongue movements, also had small negative coefficients in both models.

There were some differences between the two models. In the model including only clean samples, mounting behaviors had relatively large coefficients while poor ruminal fill and dirtiness below hocks had relatively large coefficients in the model with all samples. For the negative coefficients, poor cleanliness of the body had a larger negative coefficient in the model with only clean samples compared with the model with all samples. These differences are likely associated with correlations between variables (leading to different variables being picked in the different models) and in some cases (such as for mounting) a result of few observations.

The mean cortisol content in the samples was 15.4 pg/mg, which is in the higher range of what previous studies using an ELISA from the same manufacturer have reported (Moya et al., 2013; Burnett et al., 2014, 2015; Ghassemi Nejad et al., 2017). However, considering the peak in hair cortisol values observed in young animals this is not surprising as these other studies included older animals. As visualized in Figure 1, a significant increase in hair cortisol levels of dairy calves with a peak around 50 days of age was observed in the GAM. In the regularized regression, the variable “Fed milk,” i.e., belonging to the age group where the peak was observed, was associated with increased hair cortisol levels. The effect of young age on HCC is consistent with the increased HCC in 6-week-old calves compared with 24-week-old calves observed by Wood and Scheipl (2017).

Calves were weaned around 2 months of age on the farms included in this study, and the observed peak may be related to increased stress due to reduced milk allowance before weaning. There are also additional changes in the life of a young calf occurring during the first 2 months. For example, most of the farms kept animals in single pens for about 2 weeks before they were group housed. This transition, including social challenges and increased activity, may have influenced hair cortisol levels. Calves are also disbudded during this time period and although this is done with local analgesics in Sweden, the procedure may induce stress and pain when the analgesics wear off. Experimental studies where hair is clipped and the newly grown hair regularly sampled should provide additional information about the role of different challenges in a young calf's life.

Correlation between duplicates was high for the exact duplicates while duplicates with different levels of contamination and color showed more variations. Considering this, in combination with the results of the GAM, emphasizes that fecal contamination had a large influence on HCC and the association is supported by a dose-response curve, i.e., increasing levels of contamination are associated with a larger effect (Hill, 1965). As the risk of dirtiness on the tail is high, we suggest avoiding the use of hair from this location when studying dairy calves, contradictory to the recommendations in previous studies on cattle (Moya et al., 2013; Burnett et al., 2014; Vesel et al., 2020). Dark hair color was also associated with increased HCC while mixed samples did not differ from white hair samples. This could be a result of the relatively low number of samples with dark and mixed hair in this study, that it was more difficult to observe contamination in dark samples or be related to other unmeasured parameters (such as the thickness of hair).

Fecal contamination of the hair sample was negatively correlated with the cleanliness of the calf in the welfare assessment (i.e., the samples from clean calves were less likely to be contaminated), and the variables clustered relatively close in the cluster analysis. Due to this correlation, it is possible that some aspects of the high HCC observed in dirty samples are actually associated with calves living in a suboptimal, dirty environment. There is some support for this in cluster analysis where it appears that high HCC was generally associated with poor cleanliness, except among the younger individuals that stand out with very high cortisol levels and good cleanliness. This discrepancy between the group of young animals and others may explain the relatively large change in the coefficient for poor cleanliness of body when more samples were added in the regularized regression. With the larger sample size, the influence of young animals, and thereby the usefulness of the variable cleanliness of body, decreased. In addition, some variables with few observations, such as mounting others, also decreased in importance when additional samples were added. However, we note that the majority of the larger coefficients not related to cleanliness remained comparable between models. Thus, we argue that using a contamination score can be a way to handle samples exposed to contamination in the analysis, as long as the potential for confounding between environment and contamination is kept in mind.

It is important to note that in this study, welfare was only assessed at the one-time point and the behaviors observed during a relatively short time period. This means that calves are compared based on a snapshot during an active time period. Hair cortisol on the other hand represents an average of cortisol levels during the time of hair growth, and there are likely past events influencing hair cortisol levels that are unaccounted for in this study. Still, it is interesting that, despite the limited observation time, we find relatively short-term measures, such as poor ruminal fill and being displaced, related to increased hair cortisol levels. However, as the personality and feeding behavior of calves have been shown to be relatively constant over time (Lecorps et al., 2018; Neave et al., 2018), the short-term measures may still be related to retrospective cortisol levels. BCS is also a measure of feed intake, although more long-term, as opposed to ruminal fill, this variable was not associated with high HCC levels. In the cluster analysis, this measure clustered with another long-term measure, coat condition, and these individuals appeared to have low HCC. A possible explanation could be that this represents a group that has been exposed to stress for a longer time and exhibit a habituated cortisol response. Another possible explanation is that these long-term characteristics are associated with a poor, non-stimulating environment which causes hypostimulation (Korte et al., 2007). Longitudinal studies where confounders can be limited or controlled and HCC followed over time are needed to confirm these associations.

The usefulness of cortisol as an indicator of welfare lies in the potential to identify chronically stressed individuals struggling with coping. In this study, indicators of poor welfare were found among variables associated with high and low hair cortisol levels. The variables associated with the largest negative coefficients in the regularized regression were related to poor health while the largest positive coefficients appeared to be related to competition and dominance. However, some of the variables with smaller positive coefficients were not related to stress and poor welfare. That locomotor play was related to higher cortisol levels and may be related to activity influencing cortisol, for example, human endurance athletes have higher HCC compared with controls (Skoluda et al., 2012) or a result of the correlation between play and young age. However, the association between high reactivity and low cortisol was not anticipated as previous studies on cattle have not identified associations between HCC and temperament (Cooke et al., 2017; Lockwood et al., 2017). In many other species, fearfulness and reactivity have been associated with different coping styles and used to separate proactive and reactive individuals (Koolhaas et al., 1999). Proactive coping is characterized by a low HPA axis activity but a high sympathetic reactivity, whereas reactive coping is associated with a more active HPA axis activity and low sympathetic reactivity. The association between high reactivity and low cortisol, observed in our study, may be attributable to calves having different coping styles. The animals that were fearful and reacted violently to restraining may be the examples of proactive individuals (with a low HPA axis activity) exposed to an unfamiliar situation, whereas the reactive individuals (with high HPA axis activity) are generally more prone to handle new situations in a more passive manner (Koolhaas et al., 2007). This could explain why high cortisol levels were observed in dynamic calves interacting with other calves and the environment.

While high HCC was associated with indicators of poor and good welfare, low HCC was mainly associated with indicators related to poor welfare and reduced health. This is in agreement with a study that found lower hair cortisol levels in horses with clearly compromised welfare (Pawluski et al., 2017). It also supports the idea that the most severe forms of chronic stress may be associated with lower levels of cortisol, due to a decreased responsiveness to adrenocorticotropic hormone or habituation (Korte et al., 2007; Mormède et al., 2007). However, the association between low cortisol levels and indicators of poor clinical health (such as coughing and nasal discharge) contradicts that high cortisol levels can be useful to identify stress associated with clinical disease in calves, as previous studies on cows have suggested (Comin et al., 2013; Burnett et al., 2015). This discrepancy may also be a result of confounding due to age or environmental factors as mentioned above. A restrictive, non-stimulating environment can lead to a lower baseline of blood cortisol in healthy calves (Fisher et al., 1997). Such non-stimulating environments could be connected to environments where calves are more likely to develop coughing and other health-related problems (such as poor coat condition and BCS mentioned above). If so, low hair cortisol levels could be used as a herd level indicator of a suboptimal environment and/or poor management of calves where calves are also more likely to become sick. However, additional confirmation of these associations and increased understanding of how environment, personality, and coping styles influence hair cortisol levels are needed.