The experiment was carried out in accordance with the ARRIVE guidelines, the EU Directive 2010/63/EU for animal experiments and the protocol (19-2181) was approved by the ethical commission of Liège University. The sample size was calculated in order to use only the necessary number of cows. Among the biomarkers listed in the literature, the hair cortisol is the most frequently mentioned and was then used to calculate the sample size. Effect of stress on hair cortisol (d=(µ-µ0)/σ) was found 0.5 24 or to be higher than 1 30,31. Then an intermediate level of 0.75 was selected, and combined with an α risk of 5% and a test power (1 – β) of 80%, to reach a minimum number of 12 cows per group. To prevent from potential removing of cows during the trial for ethical reasons 3 cows were added to each group to reach a number of 15 cows per group. The total number of animals being limited to 30, the experiment was not replicated.
Animals and induction of stress
The data in this study were collected in the experimental herd of the Walloon Agricultural Research Centre (CRA-W, Gembloux, Belgium), from February to March 2020. A total of 30 cows were involved in the experiment: 25 Holstein and 5 Holstein x Simmental crossbreed F1 cows. To avoid pathological or metabolic biases, the cows were selected regarding absence of diseases and with a lactation stage greater than 30 days in milk (DIM). The cows were divided into control and stress groups of 15 cows each. Groups were constituted to have similar mean and standard deviation regarding parity, milk yield, lactation stage, and equivalent proportion of pregnant, dominant and crossbred cows in both group. Parities were comprised between 1 and 6, with an average of 2, and DIM ranged from 43 to 400, with an average of 168. All cows were originally housed in a common straw-bedded free stall barn pen, with more than 10 m² per cow and more individual places at the feedbunk than cow number. Then, for a period of 4 weeks, the 15 cows of the stress group were housed in overstocked condition, with 4.6 m² per cow, including resting and feeding areas, by moving them into a smaller straw-bedded free stall pen of the same building. In this area, only 7 feed bunks were available to generate competition for feed access. Once weekly, for 2 hours, and additional stress was induced by moving cows to an unfamiliar barn and diffusion of stressing noises (dog barking). During the 4 weeks, the 15 cows of the control group stayed in the original barn, with more than 10 m² per cow and more individual places at the feedbunk than cow number. The two pens for the stress and the control group were in the same barn, facing each other and only separated by the feeding area of 4m, with identical environment regarding exposure, temperature, materials, design or feeding times. Thus the pen effect was considered limited or non-existent. All cows received after the morning milking, approximately at 09h00, the same TMR diet composed by maize silage, grass silage and concentrate. The stress period finished at the end of the week 4, and all 30 cows were gathered into the original barn with more than 10 m² per cow and more feedbunk places than cow number. Schematic representation of the experiment is presented in Fig. 1.
Milk yield was measured daily during the experiment. The milk yield dynamic evolution of individual cows was calculated as the daily percentage change compared to the average of week 0. Body condition score (BCS) were recorded by two trained observers using a five-point scale with quarters 32 and cows were weighted weekly for stress group, and only at the beginning and at the end of the trial for the control group in order to avoid induction of stress. Clinical disease and oestrus were observed as well.
Recording of behavior
Activity and rumination were recorded continuously during all the experiment, with a time resolution of 2 hours, using the system SCR Heatime ® Pro (Allflex, Palmerston North, New Zealand). The avoidance distance test were realized following the Welfare Quality® protocol, weekly for the stress group, and only at the beginning and the end of the trial for the control group in order to avoid induction of stress. Once weekly, at the 6th day of each week as described in Fig. 1, all the 30 cows were observed during one hour, by two observers, at the rate 15 minutes per group of 15 cows, repeated four times by alternating observers. Observations took into account the interactions between animals such as given and received chasing, head butts and grooming. For further analysis, the social position of cows (i.e. dominant, neutral or dominated) was determined from the experience of the herd managers and updated after separation of the 2 groups using the difference between given and received aversive interactions observations (i.e. chasing and head-butts). The 30% cows with higher and positive differences were considered as dominant, the 30% cows with lower and negative differences were considered as dominated, and the remaining 40% cows were considered as neutral.
Heart rate variability
Heart rate and heart rate variability were measured weekly for the stress group, and only at the beginning and the end of the trial for the control group in order to avoid induction of stress. Measurements were done at the 6th day within the corresponding weeks, as described in Fig. 1. Heart rate recording were obtained using mobile Equine Polar H10 transmitters (Polar Electro Oy, Kempele, Finland) and Polar Equine belts equipped with electrodes. Signals were collected by the Polar Equine App (Polar Electro Oy, Kempele, Finland) installed on smartphones (Wiko Y50, Tinno, Shenzhen, China). The electrodes were positioned on the left side of the chest with one electrode placed close to the sternum and the other over the right scapula. The coat was first cleaned and water dampened, and electrode gel was applied to optimize electrode-skin contact. In weeks 0 and 4, belts were placed on 15 animals, half from the stress group and half from the control group, and recording were obtained from 10.00 to 12.00. From 13.00 to 15.00, recording were obtained on the remaining 15 animals. The first hour of measurement was considered as acclimatization period and associated recordings were discarded as suggested by von Borell et al. (2007). Data were first cleaned manually after visual detection of time periods with artefacts or loss of signal. Data was treated as described by Kovács et al. (2015). For analysis, 5-min time windows were selected. A total of 925 valid 5-min time windows were used for HRV analysis, 698 from stressed cows [10.0 ± 3.7 observations per date per cow] and 227 from control ones [7.9 ± 2.5 observations per date per cow]. The Kubios HRV software (version 2.1, Biomedical Signal Analysis Group, Department of Applied Physics, University of Kuopio, Finland) was used for HRV analysis. Means of heart rate, in beat per minute (bpm), and interbeat intervals (IBIs) were calculated. The root mean square of successive differences (RMSSD) between consecutive IBIs were calculated to assess the regularity of the heart rate. The correlation between successive IBIs, where each interval in the time series (IBIi + 1) is plotted against its successor (IBIi), was evaluated through Poincaré plot analysis. Standard deviation 1 (SD1) and the ratio between standard deviation 2 (SD2) and SD1 (SD2/SD1) were calculated to analyze the discontinuity and the continuity between successive IBIs, respectively.
Saliva and hair cortisol
Saliva and hair samples were collected at the 7th day within the corresponding weeks, as described in Fig. 1, right after the morning milking and before the diet distribution. Sampling were done weekly for stress group in order to follow the dynamic of the cortisol concentration, and only at the beginning and the end of the trial for the control group in order to avoid induction of stress. Therefore, the hairs collected in the week 4 for the control group corresponded to a period and length of 4 weeks growth. To compare on a similar period and length of growth, the hair cortisol concentration of the week 1 to 4 were averaged for the stress group and compared to the hairs of the week 4 of the control group. To take into account the lag time, of one to two weeks, for cortisol deposition in the hair shaft due to its initial deposition in the hair root which is beneath the skin surface (Vesel et al., 2020), a last hair sample was collected in week 5, one week after the end of the stress period.
Saliva samples were collected using a sponge held with a string and placed inside the mouth of the cows until saturated (approx. 10 ml, 1 to 2 min). The sponges were manually pressed to gather saliva, and collected samples were stored on ice, centrifuged at 2000 × g for 10 min to separate from feed particles (within 2 h after collection) and immediately frozen at -20°C until analysis for cortisol. For analysis, samples were thawed at room temperature, vortexed and centrifuged at 1500 × g for 15 minutes at 4°C. Cortisol content was determined using Salimetrics extended range salivary cortisol ELISA kit (1-3002, Salimetrics, State College, PA, USA) and following manufacturer protocol as described in Schwinn et al. (2016). Repeatability of the Elisa was 4%CV. Hair sample were collected at the extremity of the tail switch as the hair on the tail switch grows more rapidly than other sites, and is sensitive enough to capture changes in cortisol over intervals as short as 3 weeks 34. For the first hair sampling, in week 0, the hairs were first cut at a distance of 2 cm from the skin, which correspond approximately to 40 days of hair growth. Hairs were collected by shaving as close as possible to the skin with an electric clipper. For the following sampling, only the re-grown hairs of the same area were collected to avoid contamination with old hairs and observe only the cortisol deposit in hair due to the current experiment. The tail switch extremity was always entirely shaved after each sampling to maximize the re-grown hair surface to collect at the following sampling and reach 250 mg of hairs. For cows in the stress group, sampled weekly, collected regrown hairs one week after shaving were approximately 3–4 mm long, for a total weight of 486 mg per cow in average (from 220 to 670mg). The clipper was cleaned with a brush between each animal. Hair were collected in large metallic tub, dried at room temperature for 1 week, and store at -20°C until treatment. Before analysis, hairs were separated from skin follicles, dirt and faeces by mechanic sieving for 5 minutes, using 3 sieves of 400, 250 and 200 µm. The sieves were cleaned with a paintbrush between each individual samples. Then hair sample were washed and cortisol extracted using a protocol adapted from Tallo-Parra et al. (2014). From each sample, 250 mg of hair were weighed and placed into a 15-ml conical tube. In order to evaluate chronic stress, hairs were washed by adding 2.5 ml of isopropanol (2-propanol 99.5%) and vortexed at 1800 r.p.m. for 2.5 min to remove saliva, sweat, and sebum as diffusion of cortisol to these fluids is influenced by acute stress 36. The supernatant was separated by decantation and the process was repeated three times in total. The hair samples were left to dry completely for 5 days at room temperature. Washed hair samples were grounded using a ball mill, during 5 minutes at 22Hz with a 12mm metallic ball. For cortisol extraction, 50 mg of grounded hair were weighted and placed into a 2-ml eppendorf tube with 1.5 ml pure methanol and the samples were shaken at 100 r.p.m. for 18 h at 30°C. Samples were centrifuged at 7000 × g for 2 min and 0.750 ml of supernatant were transferred into a new 2-ml eppendorf tube and then placed in an oven at 38°C for 24H to evaporate methanol. The dried extracts were reconstituted with 0.25 ml buffer provided in the ELISA kit and stored at -20°C. Cortisol content was determined using Salimetrics extended range salivary cortisol ELISA kit (1-3002, Salimetrics, State College, PA, USA) following manufacturer protocol. Repeatability of the Elisa was 15.7 %CV, whereas when taking two different sub-sample after grounding, reproducibility was 28%CV, suggesting a high cortisol variation due to extraction steps or high heterogeneity within hair powder.
Blood sampling and analysis
Samples were collected to analyse β-endorphin, glycaemia, fructosamine, thyroxine (T4) and leucocytes. Blood samples were collected weekly for stress group, and only at the beginning and the end of the trial for the control group in order to avoid induction of stress. Sampling were done at the 7th day within the corresponding weeks, as described in Fig. 1, right after the saliva sampling and before the diet distribution and the hair sampling. Sample were collected at the tail vein (vena caudalis), in yellow tubes with serum separating gel for fructosamine and T4 analysis, in green heparinized tubes to harvest plasma for β-endorphin analysis, in purples tubes with EDTA for leucocytes count and in grey tubes with antiglycolytic agent for glucose analysis, and stored on ice until treatment or analysis. Fructosamine and glucose contents were analyzed with calorimetric methods (Alinity C, Abbott®), T4 was analyzed by chimiluminescence (Immulite 2000, Siemens®) and leucocytes count was realized by flow cytometry (Advia, Siemens®) at Synlab (Liège, Belgium). Treatment and analysis for β-endorphin were realized at CRA-W. Within the 30 minutes after sampling, tubes were centrifuged at 4°C, 1000 × g for 15min, and 300µl plasma were pipetted in 2ml tubes and preserved at -80°C until analysis. Analysis of β-endorphin content were realized with Mybiosource Bovine beta-endorphin ELISA kits, (MBS2000120-96, Mybiosource Inc, CA, USA) following the manufacturer protocol and repeatability of the Elisa was 18%CV.
The collected variables had different time resolution, with majority having one observation per cow per week, milk production data having one observation per cow per day, and activity and rumination having one observation per 2H per cow. The different time resolutions were harmonized by performing weekly averages for milk yield, activity and rumination. To take into account the intra-week variability of variables with high time resolutions (i.e. activity and rumination), the standard deviation (SD) per cow per week were calculated.
The main objective was to highlight biomarkers having equivalent distribution for stress and control group in week 0, and having a different level in week 4, showing thus a level modification due to stress induction. Because of the long-term effect of chronic stress, it was not possible to implement a cross-over design because the stressed cows could not be moved into a non-stressed group afterward. A common practice when cross-over is not possible, such as for studies on heat stress or diseases, is to consider the individual cows effects with repeated data in time within a mixed model 37–39. The duration of the experiment being limited to 4 weeks, the evolution of biomarkers was considered parallel among cows, with each cows having an individual intercept. For this, linear mixed repeated models were performed using the PROC MIXED procedure of SAS (SAS Institute, Cary, USA), with the random effect of cow being REPEATED along the weeks:
Yijklmn = µ + groupi + weekj + groupi*weekj + cowk + eijklmn
where Yijklmn is the observation for the potential biomarker, µ = overall mean; groupj = the fixed effect of group i (control or stress); weekj = the fixed effect of week j; cowk = the random effect of cow k and eijklmn = the experimental error. Different covariance structures were tested for each biomarker: AR(1), ARH(1), ANTE(1), CSH, TOEP, TOEPH and UN, and the lowest AIC was selected. All the records were included in the analysis. Plots of residuals were used to ensure an approximate normal distribution and no log transformations were necessary. For each week and biomarker (dependent variable), difference of least square means were used to assess the difference between the stress and control groups. In the objective to perform a large-scale assessment of chronic stress, it would be important to know what other factors are influencing the potential biomarkers. For this purpose, more complex linear mixed models were also used to evaluate the effect of DIM, parity, breed and social position:
Yijklmn = µ + groupi + weekj + groupi*weekj + cowk + b1*DIM + b2*DIM2 + parityl + breedm + socialpositionn + eijklmn
where b1 and b2 are the regression coefficients for DIM and squared DIM (DIM²); parityl= the fixed effect for parity l (1, 2 or 3+); breedm= the fixed effect for breed m (Holstein or cross-bred); socialpositionn=the fixed effect for the social position n (dominant, neutral or dominated) and eijklmn = the experimental error. The Type III sums of squares were used to determine whether these effects were significant.